function [m_Order] = calculateOrderParameter(~,psi,Transf,Params,VDk,V,T,muchem)

    NumRealiz = 100;

    Mx = numel(Transf.x);
    My = numel(Transf.y);
    Mz = numel(Transf.z);

    r = normrnd(0,1,size(psi));
    theta = rand(size(psi));
    noise = r.*exp(2*pi*1i*theta);
    
    KEop= 0.5*(Transf.KX.^2+Transf.KY.^2+Transf.KZ.^2);
    Gamma = 1-1i*Params.gamma_S;
    dt = Params.dt;

    avgpsi = 0;
    avgpsi2 = 0;

    for jj = 1:NumRealiz
        %generate initial state
        xi = sqrt(2*Params.gamma_S*Params.kbol*T*10^(-9)*dt/(Params.hbar*Params.w0*Transf.dx*Transf.dy*Transf.dz));
        swapx = randi(length(Transf.x),1,length(Transf.x));
        swapy = randi(length(Transf.y),1,length(Transf.y));
        swapz = randi(length(Transf.z),1,length(Transf.z));
        psi_j = psi + xi * noise(swapx,swapy,swapz);
        
        % --- % propagate forward in time 1 time step:
        %kin
        psi_j = fftn(psi_j);
        psi_j = psi_j.*exp(-0.5*1i*Gamma*dt*KEop);
        psi_j = ifftn(psi_j);
    
        %DDI
        frho = fftn(abs(psi_j).^2);
        Phi = real(ifftn(frho.*VDk));
    
        %Real-space
        psi_j = psi_j.*exp(-1i*Gamma*dt*(V + Params.gs*abs(psi_j).^2 + Params.gammaQF*abs(psi_j).^3 + Params.gdd*Phi - muchem));

        %kin
        psi_j = fftn(psi_j);
        psi_j = psi_j.*exp(-0.5*1i*Gamma*dt*KEop);
        psi_j = ifftn(psi_j);
    
        %Projection
        kcut = sqrt(2*Params.e_cut);
        K = (Transf.KX.^2+Transf.KY.^2+Transf.KZ.^2)<kcut.^2;
        psi_j = ifftn(K.*fftn(psi_j));
    
        % --- %
        avgpsi = avgpsi + abs(sum(psi_j(:)))/NumRealiz;
        avgpsi2 = avgpsi2 + sum(abs(psi_j(:)).^2)/NumRealiz;

    end

    m_Order = 1/sqrt(Mx*My*Mz)*avgpsi/sqrt(avgpsi2);
end