this code is for PSR calculation, added pseudo code to not recalculate atom properties matrices if they already saved

1 files changed, 150 insertions, 0 deletions

diff --git a/psr/rb87_D1_line.m b/psr/rb87_D1_line.m
new file mode 100644
index 0000000..ecdd6a3
--- /dev/null
+++ b/psr/rb87_D1_line.m
@@ -0,0 +1,150 @@
+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=0.0;
+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