Anyway implicitplot should be used only when explicit or parametric are out of the question. So:

plot([seq(arctan(surd(3*cos(x),3))+k*Pi,k=-1..1)], x=-Pi..Pi,y = -(1/2)*Pi-1 .. (1/2)*Pi+1, color=red) ;

@hitstudent 

Maple has it but the user must formulate mathematically the problem.
In this problem Maple will by really helpful if e.g. the direction of the projection is not the normal of the plane or if the sides of the square are not paralel/orthogonal to the base of the cylinder.

@Markiyan Hirnyk 

Thank you. Not a bad job done by Mathematica.
However the two samples N(1,1),  N(3,1/2) seem to be too separated, so the task is easier.

@Carl Love 

I don't think that the command works this way. It could return a convex combination of a normal and some other distribution. It would be interesting to merge two equal samples N(m1,s1), N(m2,s2) and see whether m1,m2,s1,s2 are retrived.

Maybe Markiyan could do this.

@Markiyan Hirnyk 

It was just an opinion, I am not a statistician. However, it does not seem to be  a very complicated command. In your example it returned a convex combination of 2 standard distributions. Having a finite set F of such distributions, it is not difficult to look for a distribution c*f + (1-c)*g,  f,g in F, c in [0,1]  and determine the pair (including parameters) which best fits for a given sample.

@TheFixer 

For cones and cylinder "h" is not the same. Look at my vol.

@TheFixer 

In my formula the max volume is 3.175419717

while X=Y=sqrt(2)/2 corresponds to t=Pi/4 and the volume is 2.528158531.

Yes, but probably it would be used only by curiosity, because AFAIK the user usually knows the theoretical distribution of his data and has to determine only the parameters.
Of course it would be interesting to know a statistician's opinion about this.

@nika 

To plot in the same window the histogram and the density function you need relative frequency or scale the density function.

Also, it seems that "a" does not follow very closely a Gamma distribution. (See a true sample in the attached worksheet).

gamma.mw

@nika 

Your data "a" vanished.
It would be a good idea to upload your worksheet.

@J4James 

Actually you don't need, because the fact that f is undefined at those points does not help much:

f:=x->x-floor(x/4)*4:
A:=plot(f, -12..12-0.0001, scaling=constrained, discont):
B:=plot([seq([-12+4*k,0],k=0..5), seq([-8+4*k,4],k=0..5) ],style=point, symbol=circle):
plots[display](A, B, size=[1200,400]);

@torabi 

It seems that you want to obtain some polynomial approximations for an ODE system. You impose e.g. K(0)=0 but in your system K appears at the denominator; it cannot work.

If you change
s2:=K(eta)(0) = 0.001: s3:=Omega(eta)(0) = 0.001:
then fsolve finds a solution, but I don't know if it is related to the solution of the ODE (provided it exists, which is problematic).

@Markiyan Hirnyk 

I don't think that you are in a good position to talk about seriouseness and also politeness here.

@Markiyan Hirnyk 

The post was edited. The Gamma stuff was in the second question of the OP.

@Markiyan Hirnyk 

This is not a scientific paper. The OP assumed that the sample comes from a Gamma distribution and I have shown how the parameters a,b can be obtained. That.s all.