1;
useful_constants;

%  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        

B_field=.3;

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=.0; % gyro magnetic ration for exited level % CHECKME


%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);
for j=1:Nlevels
	for k=1:Nlevels
		if ( abs(ang_momentum(j) - ang_momentum(k)) == 1)
			%transition allowed for L =L' +/- 1
			% incorrect Clebsch-Gordan coefficients but should show some ideas
			% but they correct within a factor, but not sign
			if ( ((m(j)- m(k))==1) && ( H0(j,j) > H0(k,k)) )
				dipole_elements.left(j,k)=1;
			endif	
			if ( ((m(j)- m(k))==-1) && ( H0(j,j) > H0(k,k)) )
				dipole_elements.right(j,k)=1;
			endif	
			if ( ((m(j)- m(k))==0) && ( H0(j,j) > H0(k,k)) )
				dipole_elements.linear(j,k)=1;
				if ( (m(j)==0) && ( total_momentum(j) == 1 ) && (total_momentum(k) == 1 ) )
					% remember m=0 ->m'=0 is forbidden when F=1 and F'=1
					% so we drop it back to zero
					dipole_elements.linear(j,k)=0;
				else
				endif
			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=dipole_elements.linear + dipole_elements.left + 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=.001;
for i=1:Nlevels
	for j=1:Nlevels
		if (  abs( H0(i,i) - H0(j,j)) == w_hpf_ground*hbar )
			% ground hyperfine are mixed together
			g_decay(i,j)=gamma_hpf;
		endif
	endfor
endfor

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