1;
useful_constants;

% note all frequency values are in MHz

%  87Rb D1 line
%
%     m=-2   m=-1   m=0    m=1    m=2
%     ----   ----   ----   ----   ----   |P,F=2>
%
%            ----   ----   ----          |P,F=1>
%            m=-1   m=0    m=1        
%
%
%
%
%
%     m=-2   m=-1   m=0    m=1    m=2
%     ----   ----   ----   ----   ----   |S,F=2>
%
%            ----   ----   ----          |S,F=1>
%            m=-1   m=0    m=1        


w_hpf_ground=6834;
w_hpf_exited=817;
w_sf2 = w_hpf_ground; %            Distance from |S,F=1> to |S,F=2>
w_pf1       =1e9;  % something big Distance from |S,F=1> to |P,F=1>
w_pf2       = w_pf1+w_hpf_exited; %Distance from |S,F=1> to |P,F=2> 
gmg=.7; % gyro magnetic ration for ground level
gme=.23; % gyro magnetic ration for exited level % CHECKME

zeeman_splitting=.025;
B_field=zeeman_splitting/gmg;

%bottom level   |F=1>
levels( 1)=struct( "ang_momentum",  0, "total_momentum",   1, "m",  -1,  "energy",       0, "gm",  -gmg);
levels( 2)=struct( "ang_momentum",  0, "total_momentum",   1, "m",   0,  "energy",       0, "gm",  -gmg);
levels( 3)=struct( "ang_momentum",  0, "total_momentum",   1, "m",   1,  "energy",       0, "gm",  -gmg);

%second bottom level |F=2>
levels( 4)=struct( "ang_momentum",  0, "total_momentum",   2, "m",  -2,  "energy",   w_sf2, "gm",   gmg);
levels( 5)=struct( "ang_momentum",  0, "total_momentum",   2, "m",  -1,  "energy",   w_sf2, "gm",   gmg);
levels( 6)=struct( "ang_momentum",  0, "total_momentum",   2, "m",   0,  "energy",   w_sf2, "gm",   gmg);
levels( 7)=struct( "ang_momentum",  0, "total_momentum",   2, "m",   1,  "energy",   w_sf2, "gm",   gmg);
levels( 8)=struct( "ang_momentum",  0, "total_momentum",   2, "m",   2,  "energy",   w_sf2, "gm",   gmg);

% first exited level |F=1>
levels( 9)=struct( "ang_momentum",  1, "total_momentum",   1, "m",  -1,  "energy",   w_pf1, "gm",  -gme);
levels(10)=struct( "ang_momentum",  1, "total_momentum",   1, "m",   0,  "energy",   w_pf1, "gm",  -gme);
levels(11)=struct( "ang_momentum",  1, "total_momentum",   1, "m",   1,  "energy",   w_pf1, "gm",  -gme);

% second exited level |F=2>
%levels(12)=struct( "ang_momentum",  1, "total_momentum",   2, "m",  -2,  "energy",   w_pf2, "gm",   gme);
%levels(13)=struct( "ang_momentum",  1, "total_momentum",   2, "m",  -1,  "energy",   w_pf2, "gm",   gme);
%levels(14)=struct( "ang_momentum",  1, "total_momentum",   2, "m",   0,  "energy",   w_pf2, "gm",   gme);
%levels(15)=struct( "ang_momentum",  1, "total_momentum",   2, "m",   1,  "energy",   w_pf2, "gm",   gme);
%levels(16)=struct( "ang_momentum",  1, "total_momentum",   2, "m",   2,  "energy",   w_pf2, "gm",   gme);

Nlevels=size(levels)(2);


H0=zeros(Nlevels);
energy            = [1:Nlevels].*0;
ang_momentum      = [1:Nlevels].*0;
total_momentum    = [1:Nlevels].*0;
m                 = [1:Nlevels].*0;
gm                = [1:Nlevels].*0;

for i=1:Nlevels
	energy(i)         = levels(i).energy;
	ang_momentum(i)   = levels(i).ang_momentum;
	total_momentum(i) = levels(i).total_momentum;
	m(i)              = levels(i).m;
	gm(i)             = levels(i).gm;
endfor
H0=diag(energy)*hbar;


dipole_elements.left   = zeros(Nlevels);
dipole_elements.right  = zeros(Nlevels);
dipole_elements.linear = zeros(Nlevels);
% define dipole elements for transition from level j->k
for j=1:Nlevels
	for k=1:Nlevels
		if ( abs(ang_momentum(j) - ang_momentum(k)) == 1)
			%transition allowed for L =L' +/- 1
			% but they correct within a factor, but not sign
			if ( ( H0(j,j) < H0(k,k)) )
				de=dipole_elementRb87D1line(...
					total_momentum(j),m(j), ...
					total_momentum(k),m(k)  ...
				);
				dipole_elements.left(j,k)   = de.left;
				dipole_elements.right(j,k)  = de.right;
				dipole_elements.linear(j,k) = de.linear;
			endif	
		endif
	endfor
endfor
dipole_elements.linear+=dipole_elements.linear';
dipole_elements.left+=dipole_elements.left';
dipole_elements.right+=dipole_elements.right';

maximum_dipole_elements=abs(dipole_elements.linear) + abs(dipole_elements.left) + abs(dipole_elements.right);

%defasing matrix
g_deph=0;
g_dephasing=zeros(Nlevels);
% dephasing only for non zero dipole elements (am I right?)
g_dephasing=g_deph*(abs(maximum_dipole_elements) != 0);

% decay matrix g(i,j) correspnds to decay from i-->j
gamma=6;
g_decay=zeros(Nlevels);
for i=1:Nlevels
	for j=1:Nlevels
		if (  H0(i,i) > H0(j,j) )
			% only upper levels decaying and decay is always positive
			g_decay(i,j)=gamma * abs( maximum_dipole_elements(i,j) );
		endif
	endfor
endfor

% ground level mixing need to be artificial
gamma_hpf=.0001;
for i=1:Nlevels
	for j=1:Nlevels
		% it would be better to introduce a level corresponding to the bath
		% but this slows calculations
		% So 
		% ground hyperfine are equally mixed together
		if (  (abs( H0(i,i) - H0(j,j)) == w_hpf_ground*hbar ) || ...
			    (abs( H0(i,i) - H0(j,j)) == 0 ) )
			% i and j correspond to 2 ground levels
			% we cannot decay to itself
			if ( i != j )	
				% total eight ground levels F=2 and F=1 so we decay in 7 channels
				g_decay(i,j)+=gamma_hpf/7; 
			endif
		endif
	endfor
endfor

% apply B field Zeeman splitting
energy+=B_field * ( gm .* m);
% convert frequency to energy units
H0=diag(energy)*hbar;

% vim: ts=2:sw=2:fdm=indent