clear all
kL = 1;  Gamma = 1; Delta = -3*Gamma;Omega = 0.5*Gamma;
Je = 3; Jg = 2;
vlist = [linspace(-1,-0.2,31), linspace(-0.2,0.2,801), linspace(0.2,1,31)]/kL *Gamma ;
nv = length(vlist);
flist = zeros(nv,1);

trapdepthlist = [0,0.1,1];
for d=1:3
    trapdepth=trapdepthlist(d); % mk
    for pol=0:1
    %pol = 0; % 0=circular , 1=linear
        for direction=0:1
        %direction = 0; % 0=perpendicular to dipole trap axis, 1=parallel
            if pol==0
                %circular trap, 851nm
                level_shift_F3 = [3.36025e6,2.91374e6,3.95920e6,6.49664e6,1.05260e7,1.60474e7,2.30608e7]/6.07e6*trapdepth; % in units of Gamma
                level_shift_F2 = [-2.25999e7,-2.17135e7,-2.0827e7,-1.99406e7,-1.90542e7]/6.07e6*trapdepth; % in units of Gamma
                if direction==0
                    % beams perpendicular to optical dipole trap axis
                    Omega_1m=Omega/sqrt(4); Omega_10=Omega/sqrt(2); Omega_1p=Omega/sqrt(4); 
                    Omega_2m=-Omega/sqrt(4); Omega_20=Omega/sqrt(2); Omega_2p=-Omega/sqrt(4);
                    filename=strcat(num2str(trapdepth),'mk_circulartrap_per.dat');
                else
                    % beams parallel to optical dipole trap axis
                    Omega_1m=Omega; Omega_10=0; Omega_1p=0; 
                    Omega_2m=0; Omega_20=0; Omega_2p=Omega;
                    filename=strcat(num2str(trapdepth),'mk_circulartrap_par.dat');
                end
            else
                %linear trap, 851nm
                level_shift_F3 = [2.02072e6,9.48058e6,1.39565e7,1.54485e7,1.39565e7,9.48058e6,2.02072e6]/6.07e6*trapdepth; % in units of Gamma
                level_shift_F2 = [-2.0827e7,-2.0827e7,-2.0827e7,-2.0827e7,-2.0827e7]/6.07e6*trapdepth; % in units of Gamma   
                if direction==0
                    Omega_1m=Omega; Omega_10=0; Omega_1p=0; 
                    Omega_2m=0; Omega_20=0; Omega_2p=Omega;
                    filename=strcat(num2str(trapdepth),'mk_lineartrap_per.dat');
                else
                    Omega_1m=Omega/sqrt(4); Omega_10=Omega/sqrt(2); Omega_1p=Omega/sqrt(4); 
                    Omega_2m=-Omega/sqrt(4); Omega_20=Omega/sqrt(2); Omega_2p=-Omega/sqrt(4);
                    filename=strcat(num2str(trapdepth),'mk_lineartrap_par.dat');
                end
            end

            for kv = 1:nv
               v = vlist(kv);
               force = force_pgc_FORT(v,kL,Jg,Je,Gamma,Delta, Omega_1m, Omega_10, Omega_1p,  Omega_2m, Omega_20, Omega_2p,level_shift_F2,level_shift_F3);
               flist(kv) = force(0).ampl;
               echo off
               if rem(kv,5)==0, fprintf('%d ',kv); end
            end
            fprintf('\n');
            dlmwrite(filename,[vlist',real(flist)],'delimiter','\t');
        end
    end
end

