function c0=hbar_sq_state_generator(x1,p1,x2,p2,lambda,epsilon1,epsilon2,N,hbar)
%lambda=0.0075;
alpha1=(x1+1i*p1)/sqrt(2*hbar);
%x2=sqrt(5);
%p2=sqrt(2*(Ener-x1^2/2-p1^2/2-x2^2/2-lambda*x1^2*x2^2));
alpha2=(x2+1i*p2)/sqrt(2*hbar);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
a=(1/sqrt(hbar))*spdiags(sqrt(0:N-1)',1,N,N);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
D1 = sparse(expm(alpha1*a'-alpha1'*a)); %Displacement operator
if epsilon1==0.0
psi1 = D1(:,1); %easier than mutliplying with vacuum vector
else
S1 = sparse(expm(0.5*epsilon1'*(a^2)-0.5*epsilon1*(a')^2)); %Squeezing operator
psi1=D1*S1(:,1);
end
%%%%Now, Do the same for the second mode%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
D2 = sparse(expm(alpha2*a'-alpha2'*a));
if epsilon2==0.0
psi2=D2(:,1); %*sparse(qbasis(M,1));
else
S2 = sparse(expm(0.5*epsilon2'*a^2-0.5*epsilon2*(a')^2));
psi2=D2*S2(:,1);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
c0=psi1*psi2.'; %tensor product is taken in the matrix form for the computational purpose
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end