function 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)
% Calculate the force on an atom with light fields in a 
%  sigma+ sigma- configuration.
Ne = 2*Je+1; Ng = 2*Jg+1; Nat = Ne+Ng;
% Form the lowering operators for various polarizations of light
[am,a0,ap] = murelj(Jg,Je);
Am = sparse(Nat,Nat); Am(Ne+1:Ne+Ng,1:Ne)=am; Am = qo(Am);
A0 = sparse(Nat,Nat); A0(Ne+1:Ne+Ng,1:Ne)=a0; A0 = qo(A0);
Ap = sparse(Nat,Nat); Ap(Ne+1:Ne+Ng,1:Ne)=ap; Ap = qo(Ap);
%
gp  = eseries(0.5*Omega_1p,i*kL*v) +  eseries(0.5*Omega_2p,-i*kL*v);
gm  = eseries(0.5*Omega_2m,-i*kL*v) +  eseries(0.5*Omega_1m,i*kL*v);
g0  = eseries(0.5*Omega_10,i*kL*v) +  eseries(0.5*Omega_20,-i*kL*v);
ggp = eseries(0.5*i*kL*Omega_1p,i*kL*v) + eseries(-0.5*i*kL*Omega_2p,-i*kL*v);
ggm = eseries(-0.5*i*kL*Omega_2m,-i*kL*v) +  eseries(0.5*i*kL*Omega_1m,i*kL*v);
gg0 = eseries(0.5*i*kL*Omega_10,i*kL*v) +  eseries(-0.5*i*kL*Omega_20,-i*kL*v);
% Form interaction picture Hamiltonian
% Excited state terms
excited = basis(Nat,1:Ne);
H0 = -Delta*sum(excited*excited');
% Zeeman shit
f3m3 = basis(Nat,1);
f3m2 = basis(Nat,2);
f3m1 = basis(Nat,3);
f30 = basis(Nat,4);
f3p1= basis(Nat,5);
f3p2 = basis(Nat,6);
f3p3 = basis(Nat,7);

f2m2 = basis(Nat,8);
f2m1 = basis(Nat,9);
f20 = basis(Nat,10);
f2p1 = basis(Nat,11);
f2p2= basis(Nat,12);
Hz =  level_shift_F3(1)*(f3m3*f3m3') + level_shift_F3(2)*(f3m2*f3m2') + level_shift_F3(3)*(f3m1*f3m1')+level_shift_F3(4)*(f30*f30') + level_shift_F3(5)*(f3p1*f3p1') +level_shift_F3(6)*(f3p2*f3p2') +level_shift_F3(7)*(f3p3*f3p3') +  level_shift_F2(1)*(f2m2*f2m2') + level_shift_F2(2)*(f2m1*f2m1') + level_shift_F2(3)*(f20*f20') + level_shift_F2(4)*(f2p1*f2p1')+ level_shift_F2(5)*(f2p2*f2p2') ; 

% Interaction with sigma+ light
Hp = -(gp'*Ap + gp*Ap');
% Interaction with sigma- light
Hm = -(gm'*Am + gm*Am');
% Interaction with pi light
Hpi = -(g0'*A0 + g0*A0');

H  = H0 + Hp + Hm + Hpi + Hz;
% Force operator
Fp = (ggp'*Ap + ggp*Ap');
Fm = (ggm'*Am + ggm*Am');
F0 = (gg0'*A0 + gg0*A0');
F = Fp + Fm + F0;
% Collapse operators as an array
C = sqrt(Gamma)*[Am,A0,Ap];
% Loss Liouvillians
Lloss = spre(C)*spost(C')-0.5*spre(C'*C)-0.5*spost(C'*C);
% Calculate total Liouvillian 
L = -1i*(spre(H)-spost(H)) + sum(Lloss);
% Matrix continued fraction calculation
rho = mcfrac(L,20,1);
force = expect(F,rho);
