summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--psr/basis_transformation.m57
-rw-r--r--psr/dipole_elementRb87D1line.m177
-rw-r--r--psr/field_description.m14
-rw-r--r--psr/output_results.m31
-rw-r--r--psr/output_xi_results.m46
-rw-r--r--psr/psr.m149
-rw-r--r--psr/rb87_D1_line.m150
-rw-r--r--psr/useful_constants.m6
-rw-r--r--psr/useful_functions.m337
9 files changed, 967 insertions, 0 deletions
diff --git a/psr/basis_transformation.m b/psr/basis_transformation.m
new file mode 100644
index 0000000..82945fe
--- /dev/null
+++ b/psr/basis_transformation.m
@@ -0,0 +1,57 @@
+1;
+
+% matrix of circular to linear transformation
+% [x, y, z]' = lin2circ * [r, l, z]'
+function transformation_matrix = circ2lin()
+ transformation_matrix = ...
+ [ ...
+ [ 1/sqrt(2), 1/sqrt(2), 0]; ...
+ [-1i/sqrt(2), 1i/sqrt(2), 0]; ...
+ [ 0, 0, 1] ...
+ ];
+endfunction
+
+
+% matrix of linear to circular transformation
+% [r, l, z]' = lin2circ * [x, y, z]'
+function transformation_matrix = lin2circ()
+ transformation_matrix = ...
+ [ ...
+ [ 1/sqrt(2), 1i/sqrt(2), 0]; ...
+ [ 1/sqrt(2), -1i/sqrt(2), 0]; ...
+ [ 0, 0, 1] ...
+ ];
+endfunction
+% linear basis rotation
+% x axis untouched
+% z and y rotated by angle theta around 'x' axis
+% [x_new, y_new, z_new]' = oldlin2newlin * [x_old, y_old, z_old]'
+function oldlin2newlin_m = oldlin2newlin(theta)
+ oldlin2newlin_m = [ ...
+ [ 1, 0, 0]; ...
+ [ 0, cos(theta), -sin(theta)]; ...
+ [ 0, sin(theta), cos(theta)]...
+ ];
+endfunction
+
+% rotate x polarized light by angle phi around
+% light propagation axis (Z)
+function [E_field_x, E_field_y] = rotXpolarization(phi, E_field_linear)
+ % important negative frequency behave as they rotate in opposite direction
+ E_field_x=cos(phi)*E_field_linear;
+ E_field_y=sin(phi)*E_field_linear;
+endfunction
+
+% transform x,y,z linearly polarized light in the lab/light system coordinate
+% to left, right, linear along z atom system of coordinate
+% atom magnetic field is along new axis Z wich is at angle theta with respect to
+% light propagation direction
+function E_field_pos_freq=xyz_lin2atomic_axis_polarization(theta, E_field_lab_pos_freq)
+ coord_transf_m = lin2circ() * oldlin2newlin(theta);
+ E_field_pos_freq.right = coord_transf_m(1,1)*E_field_lab_pos_freq.x + coord_transf_m(1,2)*E_field_lab_pos_freq.y + coord_transf_m(1,3)*E_field_lab_pos_freq.z;
+ E_field_pos_freq.left = coord_transf_m(2,1)*E_field_lab_pos_freq.x + coord_transf_m(2,2)*E_field_lab_pos_freq.y + coord_transf_m(2,3)*E_field_lab_pos_freq.z;
+ E_field_pos_freq.linear = coord_transf_m(3,1)*E_field_lab_pos_freq.x + coord_transf_m(3,2)*E_field_lab_pos_freq.y + coord_transf_m(3,3)*E_field_lab_pos_freq.z;
+endfunction
+
+
+% vim: ts=2:sw=2:fdm=indent
diff --git a/psr/dipole_elementRb87D1line.m b/psr/dipole_elementRb87D1line.m
new file mode 100644
index 0000000..9235274
--- /dev/null
+++ b/psr/dipole_elementRb87D1line.m
@@ -0,0 +1,177 @@
+function d=dipole_elementRb87D1line(Fl,ml,Fu,mu)
+% Fl, ml are F and m quantum numbers of lower state
+% Fu, mu are F and m quantum numbers of upper state
+% F is total momentum and m is projection
+ d.left = 0; %default return value
+ d.linear = 0; %default return value
+ d.right = 0; %default return value
+ if ( mu==(ml+1) )
+ % sigma plus polarization
+ % ------ Fl=2 -> Fu=2 --------
+ if ( (ml==-2) & (Fl==2) & (Fu==2) )
+ d.right=sqrt(1/6);
+ endif
+ if ( (ml==-1) & (Fl==2) & (Fu==2) )
+ d.right=sqrt(1/4);
+ endif
+ if ( (ml== 0) & (Fl==2) & (Fu==2) )
+ d.right=sqrt(1/4);
+ endif
+ if ( (ml== 1) & (Fl==2) & (Fu==2) )
+ d.right=sqrt(1/6);
+ endif
+ if ( (ml== 2) & (Fl==2) & (Fu==2) )
+ d.right=0;
+ endif
+ % ------ Fl=2 -> Fu=1 --------
+ if ( (ml==-2) & (Fl==2) & (Fu==1) )
+ d.right=sqrt(1/2);
+ endif
+ if ( (ml==-1) & (Fl==2) & (Fu==1) )
+ d.right=sqrt(1/4);
+ endif
+ if ( (ml== 0) & (Fl==2) & (Fu==1) )
+ d.right=sqrt(1/12);
+ endif
+ if ( (ml== 1) & (Fl==2) & (Fu==1) )
+ d.right=0;
+ endif
+ if ( (ml== 2) & (Fl==2) & (Fu==1) )
+ d.right=0;
+ endif
+ % ------ Fl=1 -> Fu=2 --------
+ if ( (ml==-1) & (Fl==1) & (Fu==2) )
+ d.right = -sqrt(1/12);
+ endif
+ if ( (ml== 0) & (Fl==1) & (Fu==2) )
+ d.right = -sqrt(1/4);
+ endif
+ if ( (ml== 1) & (Fl==1) & (Fu==2) )
+ d.right = -sqrt(1/2);
+ endif
+ % ------ Fl=1 -> Fu=1 --------
+ if ( (ml==-1) & (Fl==1) & (Fu==1) )
+ d.right = -sqrt(1/12);
+ endif
+ if ( (ml== 0) & (Fl==1) & (Fu==1) )
+ d.right = -sqrt(1/12);
+ endif
+ if ( (ml== 1) & (Fl==1) & (Fu==1) )
+ d.right = 0;
+ endif
+ endif
+ if ( mu==(ml+0) )
+ % pi polarization
+ % ------ Fl=2 -> Fu=2 --------
+ if ( (ml==-2) & (Fl==2) & (Fu==2) )
+ d.linear=-sqrt(1/3);
+ endif
+ if ( (ml==-1) & (Fl==2) & (Fu==2) )
+ d.linear=-sqrt(1/12);
+ endif
+ if ( (ml== 0) & (Fl==2) & (Fu==2) )
+ d.linear=0;
+ endif
+ if ( (ml== 1) & (Fl==2) & (Fu==2) )
+ d.linear=sqrt(1/12);
+ endif
+ if ( (ml== 2) & (Fl==2) & (Fu==2) )
+ d.linear=sqrt(1/3);
+ endif
+ % ------ Fl=2 -> Fu=1 --------
+ if ( (ml==-2) & (Fl==2) & (Fu==1) )
+ d.linear = 0;
+ endif
+ if ( (ml==-1) & (Fl==2) & (Fu==1) )
+ d.linear=sqrt(1/4);
+ endif
+ if ( (ml== 0) & (Fl==2) & (Fu==1) )
+ d.linear=sqrt(1/3);
+ endif
+ if ( (ml== 1) & (Fl==2) & (Fu==1) )
+ d.linear=sqrt(1/4);
+ endif
+ if ( (ml== 2) & (Fl==2) & (Fu==1) )
+ d.linear = 0;
+ endif
+ % ------ Fl=1 -> Fu=2 --------
+ if ( (ml==-1) & (Fl==1) & (Fu==2) )
+ d.linear = sqrt(1/4);
+ endif
+ if ( (ml== 0) & (Fl==1) & (Fu==2) )
+ d.linear = sqrt(1/2);
+ endif
+ if ( (ml== 1) & (Fl==1) & (Fu==2) )
+ d.linear = sqrt(1/4);
+ endif
+ % ------ Fl=1 -> Fu=1 --------
+ if ( (ml==-1) & (Fl==1) & (Fu==1) )
+ d.linear = sqrt(1/12);
+ endif
+ if ( (ml== 0) & (Fl==1) & (Fu==1) )
+ d.linear = 0;
+ endif
+ if ( (ml== 1) & (Fl==1) & (Fu==1) )
+ d.linear = -sqrt(1/12);
+ endif
+ endif
+ if ( mu==(ml-1) )
+ % sigma minus polarization
+ % ------ Fl=2 -> Fu=2 --------
+ if ( (ml==-2) & (Fl==2) & (Fu==2) )
+ d.left = 0;
+ endif
+ if ( (ml==-1) & (Fl==2) & (Fu==2) )
+ d.left = -sqrt(1/6);
+ endif
+ if ( (ml== 0) & (Fl==2) & (Fu==2) )
+ d.left = -sqrt(1/4);
+ endif
+ if ( (ml== 1) & (Fl==2) & (Fu==2) )
+ d.left = -sqrt(1/4);
+ endif
+ if ( (ml== 2) & (Fl==2) & (Fu==2) )
+ d.left = -sqrt(1/6);
+ endif
+ % ------ Fl=2 -> Fu=1 --------
+ if ( (ml==-2) & (Fl==2) & (Fu==1) )
+ d.left = 0;
+ endif
+ if ( (ml==-1) & (Fl==2) & (Fu==1) )
+ d.left = 0;
+ endif
+ if ( (ml== 0) & (Fl==2) & (Fu==1) )
+ d.left = sqrt(1/12);
+ endif
+ if ( (ml== 1) & (Fl==2) & (Fu==1) )
+ d.left = sqrt(1/4);
+ endif
+ if ( (ml== 2) & (Fl==2) & (Fu==1) )
+ d.left = sqrt(1/2);
+ endif
+ % ------ Fl=1 -> Fu=2 --------
+ if ( (ml==-1) & (Fl==1) & (Fu==2) )
+ d.left = -sqrt(1/2);
+ endif
+ if ( (ml== 0) & (Fl==1) & (Fu==2) )
+ d.left = -sqrt(1/4);
+ endif
+ if ( (ml== 1) & (Fl==1) & (Fu==2) )
+ d.left = -sqrt(1/12);
+ endif
+ % ------ Fl=1 -> Fu=1 --------
+ if ( (ml==-1) & (Fl==1) & (Fu==1) )
+ d.left = 0;
+ endif
+ if ( (ml== 0) & (Fl==1) & (Fu==1) )
+ d.left = sqrt(1/12);
+ endif
+ if ( (ml== 1) & (Fl==1) & (Fu==1) )
+ d.left = sqrt(1/12);
+ endif
+ endif
+endfunction
+
+
+
+% vim: ts=2:sw=2:fdm=indent
diff --git a/psr/field_description.m b/psr/field_description.m
new file mode 100644
index 0000000..738ebe3
--- /dev/null
+++ b/psr/field_description.m
@@ -0,0 +1,14 @@
+1;
+
+%EM field definition
+Ed=10.12; %drive
+Edc=conj(Ed);
+Ep=0.8*Ed; %probe
+Epc=conj(Ep);
+Em=-Ep; % opposite sideband (resulting from EOM modulation of drive)
+Emc=conj(Em);
+%wd=w13;
+%wp=w12;
+%wm=wd-(wp-wd);
+%light_positive_freq = [wp, wd, wp-wd];
+E_field = [Ep, Ed, 0 ];
diff --git a/psr/output_results.m b/psr/output_results.m
new file mode 100644
index 0000000..eca6085
--- /dev/null
+++ b/psr/output_results.m
@@ -0,0 +1,31 @@
+1;
+
+
+load '/tmp/relative_transmission_vs_detuning.mat' ;
+
+figure(1);
+hold off;
+zoom_factor=1;
+plot(detuning_freq, zoom_factor*(relative_transmission_vs_detuning), '-');
+title("relative transmission");
+xlabel("two photon detuning");
+
+%load 'xi_vs_detuning.mat' ;
+
+%figure(1);
+ %hold off;
+ %plot(detuning_freq, imag(xi_linear), '-1;linear;');
+ %hold on;
+ %plot(detuning_freq, imag(xi_left), '-2;left;');
+ %plot(detuning_freq, imag(xi_right), '-3;right;');
+ %title("probe absorption");
+ %hold off;
+%figure(2);
+ %hold off;
+ %plot(detuning_freq, real(xi_linear), '-1;linear;');
+ %hold on;
+ %plot(detuning_freq, real(xi_left), '-2;left;');
+ %plot(detuning_freq, real(xi_right), '-3;right;');
+ %title("probe dispersion");
+ %hold off;
+
diff --git a/psr/output_xi_results.m b/psr/output_xi_results.m
new file mode 100644
index 0000000..7dc3855
--- /dev/null
+++ b/psr/output_xi_results.m
@@ -0,0 +1,46 @@
+1;
+
+
+load '/tmp/xi_vs_detuning.mat' ;
+
+figure(1);
+hold off;
+plot(detuning_freq, real(xi_left-xi_right), '-');
+title("differential real xi");
+xlabel("two photon detuning");
+
+figure(2);
+hold off;
+plot(detuning_freq, imag(xi_left-xi_right), '-');
+title("differential imag xi");
+xlabel("two photon detuning");
+
+figure(3);
+hold off;
+plot(detuning_freq, imag(xi_left), '-', detuning_freq, imag(xi_right), '-');
+title("imag xi");
+xlabel("two photon detuning");
+
+figure(4);
+hold off;
+plot(detuning_freq, real(xi_left), '-', detuning_freq, real(xi_right), '-');
+title("real xi");
+xlabel("two photon detuning");
+
+%figure(1);
+ %hold off;
+ %plot(detuning_freq, imag(xi_linear), '-1;linear;');
+ %hold on;
+ %plot(detuning_freq, imag(xi_left), '-2;left;');
+ %plot(detuning_freq, imag(xi_right), '-3;right;');
+ %title("probe absorption");
+ %hold off;
+%figure(2);
+ %hold off;
+ %plot(detuning_freq, real(xi_linear), '-1;linear;');
+ %hold on;
+ %plot(detuning_freq, real(xi_left), '-2;left;');
+ %plot(detuning_freq, real(xi_right), '-3;right;');
+ %title("probe dispersion");
+ %hold off;
+
diff --git a/psr/psr.m b/psr/psr.m
new file mode 100644
index 0000000..7c8f604
--- /dev/null
+++ b/psr/psr.m
@@ -0,0 +1,149 @@
+1;
+clear all;
+t0 = clock (); % we will use this latter to calculate elapsed time
+
+
+% load useful functions;
+useful_functions;
+
+% some physical constants
+useful_constants;
+
+basis_transformation; % load subroutines
+
+% load atom energy levels and decay description
+rb87_D1_line;
+%four_levels_with_polarization;
+%four_levels;
+%three_levels;
+%two_levels;
+
+% load EM field description
+field_description;
+
+%Nfreq=length(modulation_freq);
+
+
+
+%tune probe frequency
+detuning_p=0;
+N_detun_steps=100;
+%detuning_p_min=-B_field*gmg*4; % span +/-4 Zeeman splitting
+detuning_p_min=-10040.0;
+detuning_p_max=-detuning_p_min;
+detuning_freq=zeros(1,N_detun_steps+1);
+kappa_p =zeros(1,N_detun_steps+1);
+kappa_m =zeros(1,N_detun_steps+1);
+detun_step=(detuning_p_max-detuning_p_min)/N_detun_steps;
+
+fprintf (stderr, "calculating atom properties\n");
+fflush (stderr);
+pfile='rb87_D1_line.m'; % the parent file from which L0_and_polarization_submatrices calculated
+cfile='L0m_and_polarizability_calculated.mat'; % the child file to which calculated matrices writen
+[s, err, msg] = stat (pfile);
+if(err)
+ %file does not exist
+ disp('Big troubles are coming, no file to define Hamiltonian)');
+ msg=cstrcat('File: ', pfile, ' is missing...exiting');
+ disp(msg);
+ return;
+else
+ pfile_mtime=s.mtime;
+endif
+[s, err, msg] = stat (cfile);
+if(err)
+ %file does not exist
+ cfile_mtime=0;
+else
+ cfile_mtime=s.mtime;
+endif;
+if ( cfile_mtime >= pfile_mtime)
+ % matrices already calculated and up to date, all we need to load them
+ load(cfile);
+ else
+ % calculate E_field independent properties of the atom
+ % to be used as sub matrix templates for Liouville operator matrix
+ [L0m, polarizability_m]=L0_and_polarization_submatrices( ...
+ Nlevels, ...
+ H0, g_decay, g_dephasing, dipole_elements ...
+ );
+ save(cfile, 'L0m', 'polarizability_m');
+ endif
+elapsed_time = etime (clock (), t0);
+fprintf (stderr, "elapsed time so far is %.3f sec\n",elapsed_time);
+fflush (stderr);
+
+global atom_properties;
+atom_properties.L0m=L0m;
+atom_properties.polarizability_m=polarizability_m;
+atom_properties.dipole_elements=dipole_elements;
+
+%light_positive_freq = [wp];
+E_field_drive = [0 ];
+E_field_probe = [Ep ];
+E_field_zero = [0 ];
+E_field_lab_pos_freq.linear = E_field_zero + (1.00000+0.00000i)*E_field_probe + (1.00000+0.00000i)*E_field_drive;
+%E_field_lab_pos_freq.right = E_field_zero + (0.00000+0.00000i)*E_field_probe + (0.00000+0.00000i)*E_field_drive;
+%E_field_lab_pos_freq.left = E_field_zero + (0.00000+0.00000i)*E_field_probe + (0.00000+0.00000i)*E_field_drive;
+
+% phi is angle between linear polarization and axis x
+phi=pi*2/8;
+% theta is angle between lab z axis (light propagation direction) and magnetic field axis (z')
+theta=0;
+% psi_el is the ellipticity parameter (phase difference between left and right polarization)
+psi_el=-3/180*pi;
+
+% we define light as linearly polarized
+% where phi is angle between light polarization and axis x
+ % only sign of modulation frequency is important now
+ % we define actual frequency later on
+ [E_field_lab_pos_freq.x, E_field_lab_pos_freq.y] = rotXpolarization(phi, E_field_lab_pos_freq.linear);
+ E_field_lab_pos_freq.z=E_field_zero;
+
+ E_field_pos_freq=xyz_lin2atomic_axis_polarization(theta, E_field_lab_pos_freq);
+ E_field_pos_freq.left*=exp(I*psi_el);
+
+
+fprintf (stderr, "tuning laser in forloop to set conditions vs detuning\n");
+fflush (stderr);
+for detuning_p_cntr=1:N_detun_steps+1;
+ wp0=w_pf1-w_sf2; %Fg=2 -> Fe=1
+ %wd=w_pf1-w_hpf_ground;
+ detuning_p=detuning_p_min+detun_step*(detuning_p_cntr-1);
+ wp=wp0+detuning_p;
+ light_positive_freq=[ wp];
+ % we calculate dc and negative frequiencies as well as amplitudes
+ [modulation_freq, E_field] = ...
+ light_positive_frequencies_and_amplitudes2full_set_of_modulation_frequencies_and_amlitudes(...
+ light_positive_freq, E_field_pos_freq);
+ freq_index=freq2index(wp,modulation_freq);
+
+ atom_field_problem.E_field = E_field;
+ atom_field_problem.modulation_freq = modulation_freq;
+ atom_field_problem.freq_index = freq_index;
+
+ problems_cell_array{detuning_p_cntr}=atom_field_problem;
+
+
+ %kappa_p(detuning_p_cntr)=susceptibility_steady_state_at_freq( atom_field_problem);
+ detuning_freq(detuning_p_cntr)=detuning_p;
+endfor
+
+save '/tmp/problem_definition.mat' problems_cell_array atom_properties detuning_freq ;
+fprintf (stderr, "now really hard calculations begin\n");
+fflush (stderr);
+% once we define all problems the main job is done here
+%kappa_p=cellfun( @susceptibility_steady_state_at_freq, problems_cell_array);
+%kappa_p=parcellfun(2, @susceptibility_steady_state_at_freq, problems_cell_array);
+[xi_linear, xi_left, xi_right]=parcellfun(2, @susceptibility_steady_state_at_freq, problems_cell_array);
+%relative_transmission_vs_detuning=parcellfun(2, @total_relative_transmission, problems_cell_array);
+%relative_transmission_vs_detuning=cellfun(@total_relative_transmission, problems_cell_array);
+
+save '/tmp/xi_vs_detuning.mat' detuning_freq xi_linear xi_left xi_right ;
+%save '/tmp/relative_transmission_vs_detuning.mat' detuning_freq relative_transmission_vs_detuning;
+
+output_xi_results;
+
+elapsed_time = etime (clock (), t0)
+
+% vim: ts=2:sw=2:fdm=indent
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
diff --git a/psr/useful_constants.m b/psr/useful_constants.m
new file mode 100644
index 0000000..82d1d1a
--- /dev/null
+++ b/psr/useful_constants.m
@@ -0,0 +1,6 @@
+1;
+
+% some physical constants
+hbar=1;
+im_one=0+1i;
+el_charge=1;
diff --git a/psr/useful_functions.m b/psr/useful_functions.m
new file mode 100644
index 0000000..eb3e880
--- /dev/null
+++ b/psr/useful_functions.m
@@ -0,0 +1,337 @@
+1;
+
+function ret=decay_total(g_decay,i)
+% calculate total decay for particular level taking in account all branches
+ ret=sum(g_decay(i,:));
+endfunction
+
+function ret=kron_delta(i,j)
+% kroneker delta symbol
+ if ((i==j))
+ ret=1;
+ else
+ ret=0;
+ endif
+endfunction
+
+function rho=rhoOfFreq(rhoLiouville, freqIndex, Nlevels)
+% this function create from Liouville density vector
+% the density matrix with given modulation frequency
+ rho=zeros(Nlevels);
+ rho(:)=rhoLiouville((freqIndex-1)*Nlevels^2+1:(freqIndex)*Nlevels^2);
+ rho=rho.';
+endfunction
+
+function [N, rhoLiouville_w, rhoLiouville_r, rhoLiouville_c]=unfold_density_matrix(Nlevels,Nfreq)
+% unwrap density matrix to Liouville density vector and assign all possible
+% modulation frequencies as well
+% resulting vector should be Nlevels x Nlevels x length(modulation_freq)
+ N = Nfreq*Nlevels*Nlevels;
+ rho_size = Nlevels*Nlevels;
+ rhoLiouville_w=zeros(N,1);
+ rhoLiouville_r=zeros(N,1);
+ rhoLiouville_c=zeros(N,1);
+
+ w=1:Nfreq;
+ w_tmplate=(repmat(w,rho_size,1))(:);
+ rhoLiouville_w=w_tmplate;
+ r=1:Nlevels;
+ r_tmplate=(repmat(r,Nlevels,1))(:);
+ rhoLiouville_r=(repmat(r_tmplate,Nfreq,1))(:)';
+ c=(1:Nlevels)';% hold column value of rho_rc
+ rhoLiouville_c=repmat(c,Nfreq*Nlevels,1);
+endfunction
+
+
+function [L0m, polarizability_m]=L0_and_polarization_submatrices( ...
+ Nlevels, ...
+ H0, g_decay, g_dephasing, dipole_elements ...
+ )
+% create (Nlevels*Nlevels)x*(Nlevels*Nlevels)
+% sub matrices of Liouville operator
+% which repeat themselves for each modulation frequency
+% based on recipe from Eugeniy Mikhailov thesis
+ %-------------------------
+ useful_constants;
+ rho_size=Nlevels*Nlevels;
+
+ % now we create Liouville indexes list
+ [Ndummy, rhoLiouville_w_notused, rhoLiouville_r, rhoLiouville_c]=unfold_density_matrix(Nlevels,1);
+
+ kron_delta_m=eye(Nlevels);
+ % note that L0 and decay parts depend only on combination of indexes
+ % jk,mn but repeats itself for every frequency
+ L0m=zeros(rho_size); % (Nlevels^2)x(Nlevels^2) matrix
+ decay_part_m=zeros(rho_size); % (NxN)x(NxN) matrix
+ % polarization matrix will be multiplied by field amplitude letter
+ % polarization is part of perturbation part of Hamiltonian
+ polarizability_m.linear = zeros(rho_size); % (NxN)x(NxN) matrix
+ polarizability_m.left = zeros(rho_size); % (NxN)x(NxN) matrix
+ polarizability_m.right = zeros(rho_size); % (NxN)x(NxN) matrix
+ for p=1:rho_size
+ % p= j*Nlevels+k
+ % this might speed up stuff since less matrix passed back and force
+ j=rhoLiouville_r(p);
+ k=rhoLiouville_c(p);
+ for s=1:rho_size
+ % s= m*Nlevels+n
+ m=rhoLiouville_r(s);
+ n=rhoLiouville_c(s);
+
+ % calculate unperturbed part (Hamiltonian without EM field)
+ L0m(p,s)=H0(j,m)*kron_delta_m(k,n)-H0(n,k)*kron_delta_m(j,m);
+ decay_part_m(p,s)= ...
+ ( ...
+ decay_total(g_decay,k)/2 ...
+ + decay_total(g_decay,j)/2 ...
+ + g_dephasing(j,k) ...
+ )* kron_delta_m(j,m)*kron_delta_m(k,n) ...
+ - kron_delta_m(m,n)*kron_delta_m(j,k)*g_decay(m,j) ;
+ polarizability_m.linear(p,s)= ( dipole_elements.linear(j,m)*kron_delta_m(k,n)-dipole_elements.linear(n,k)*kron_delta_m(j,m) );
+ polarizability_m.left(p,s)= ( dipole_elements.left(j,m)*kron_delta_m(k,n)-dipole_elements.left(n,k)*kron_delta_m(j,m) );
+ polarizability_m.right(p,s)= ( dipole_elements.right(j,m)*kron_delta_m(k,n)-dipole_elements.right(n,k)*kron_delta_m(j,m) );
+ endfor
+ endfor
+ L0m=-im_one/hbar*L0m - decay_part_m;
+endfunction
+
+function L=Liouville_operator_matrix( ...
+ N, ...
+ L0m, polarizability_m, ...
+ E_field, ...
+ modulation_freq, rhoLiouville_w, rhoLiouville_r, rhoLiouville_c ...
+ )
+% Liouville operator matrix construction
+% based on recipe from Eugeniy Mikhailov thesis
+ %-------------------------
+ useful_constants;
+ L=zeros(N); % NxN matrix
+ Nfreq=length(modulation_freq);
+
+ % Lets be supper smart and speed up L matrix construction
+ % since it has a lot of voids.
+ % By creation of rhoLiouville we know that there are
+ % consequent chunks of rho_ij modulated with same frequency
+ % this means that rhoLiouville is split in to Nfreq chunks
+ % with length Nlevels*Nlevels=N/Nfreq
+ rho_size=N/Nfreq;
+
+ % creating building blocks of L by rho_size * rho_size
+ for w3i=1:Nfreq
+ w_iner=modulation_freq(w3i);
+ if ((w_iner == 0))
+ % calculate unperturbed part (Hamiltonian without EM field)
+ L_sub{w3i}=L0m;
+ else
+ % calculate perturbed part (Hamiltonian with EM field)
+ % in other word interactive part of Hamiltonian
+ L_sub{w3i} = ...
+ -im_one/hbar*polarizability_m.linear * E_field.linear(w3i) ...
+ -im_one/hbar*polarizability_m.left * E_field.left(w3i) ...
+ -im_one/hbar*polarizability_m.right * E_field.right(w3i) ...
+ ;
+ endif
+ endfor
+
+ % Liouville matrix operator has Nlevels*Nlevels blocks
+ % which governed by the same modulation frequency
+ for p_freq_cntr=1:Nfreq
+ p0=1+(p_freq_cntr-1)*rho_size;
+ % we guaranteed to know frequency of final and initial rhoLiouville
+ w1i=rhoLiouville_w(p0); % final
+ w_jk=modulation_freq(w1i);
+ for s_freq_cntr=1:Nfreq
+ s0=1+(s_freq_cntr-1)*rho_size;
+ w2i=rhoLiouville_w(s0); % initial
+ w_mn=modulation_freq(w2i);
+ % thus we know L matrix element frequency which we need to match
+ w_l=w_jk-w_mn;
+ % lets search this frequency in the list of available frequencies
+ % but since we not guaranteed to find it lets assign temporary 0 to Liouville matrix element
+ w3i=(w_l == modulation_freq);
+ if (any(w3i))
+ % yey, requested modulation frequency exist
+ % lets do L sub matrix filling
+ % at most we should have only one matching frequency
+ w_iner=modulation_freq(w3i);
+ L(p0:p0+rho_size-1,s0:s0+rho_size-1) = L_sub{w3i};
+
+ endif
+ endfor
+ % diagonal elements are self modulated
+ % due to rotating wave approximation
+ L(p0:p0+rho_size-1,p0:p0+rho_size-1)+= -im_one*w_jk*eye(rho_size);
+ endfor
+endfunction
+
+
+function [rhoLiouville_dot, L]=constrain_rho_and_match_L( ...
+ N, L, ...
+ modulation_freq, rhoLiouville_w, rhoLiouville_r, rhoLiouville_c)
+% now generally rhoL_dot=0=L*rhoL has infinite number of solutions
+% since we always can resclale rho vector with arbitrary constant
+% lets constrain our density matrix with some physical meaning
+% sum(rho_ii)=1 (sum of all populations (with zero modulation frequency) scales to 1
+% we will replace first row of Liouville operator with this condition
+% thus rhoLiouville_dot(1)=1
+ for i=1:N
+ w2i=rhoLiouville_w(i);
+ m=rhoLiouville_r(i);
+ n=rhoLiouville_c(i);
+ w=modulation_freq(w2i);
+ if ((w==0) & (m==n))
+ L(1,i)=1;
+ else
+ L(1,i)=0;
+ endif
+ endfor
+ rhoLiouville_dot= zeros(N,1);
+ % sum(rho_ii)=1 (sum of all populations (with zero modulation frequency) scales to 1
+ % we will replace first row of Liouville operator with this condition
+ % thus rhoLiouville_dot(1)=1
+ rhoLiouville_dot(1)=1;
+endfunction
+
+function kappa=susceptibility(wi, rhoLiouville, dipole_elements, E_field)
+% calculate susceptibility for the field at given frequency index
+ Nlevels=( size(dipole_elements.linear)(1) );
+ rho=rhoOfFreq(rhoLiouville, wi, Nlevels);
+ kappa.linear=0;
+ kappa.left=0;
+ kappa.right=0;
+
+ kappa.linear = sum( sum( transpose(dipole_elements.linear) .* rho ) );
+ kappa.left = sum( sum( transpose(dipole_elements.left) .* rho ) );
+ kappa.right = sum( sum( transpose(dipole_elements.right) .* rho ) );
+
+ kappa.linear /= E_field.linear( wi );
+ kappa.left /= E_field.left( wi );
+ kappa.right /= E_field.right( wi );
+endfunction
+
+function index=freq2index(freq, modulation_freq)
+% convert modulation freq to its index in the modulation_freq vector
+ index=[1:length(modulation_freq)](modulation_freq==freq);
+endfunction
+
+function rhoLiouville=rhoLiouville_steady_state(L0m, polarizability_m, E_field, modulation_freq)
+% calculates rhoLiouville vector assuming steady state situation and normalization of rho_ii to 1
+ Nlevels=sqrt( size(L0m)(1) );
+ Nfreq=length(modulation_freq);
+
+ % now we create Liouville indexes list
+ [N, rhoLiouville_w, rhoLiouville_r, rhoLiouville_c]=unfold_density_matrix(Nlevels,Nfreq);
+
+ % Liouville operator matrix construction
+ L=Liouville_operator_matrix(
+ N,
+ L0m, polarizability_m,
+ E_field,
+ modulation_freq, rhoLiouville_w, rhoLiouville_r, rhoLiouville_c
+ );
+
+ %use the fact that sum(rho_ii)=1 to constrain solution
+ [rhoLiouville_dot, L]=constrain_rho_and_match_L(
+ N, L,
+ modulation_freq, rhoLiouville_w, rhoLiouville_r, rhoLiouville_c);
+
+ %solving for density matrix vector
+ rhoLiouville=L\rhoLiouville_dot;
+endfunction
+
+function [xi_linear, xi_left, xi_right]=susceptibility_steady_state_at_freq( atom_field_problem)
+ % find steady state susceptibility at particular modulation frequency element
+ % at given E_field
+ global atom_properties;
+ L0m = atom_properties.L0m ;
+ polarizability_m = atom_properties.polarizability_m ;
+ dipole_elements = atom_properties.dipole_elements ;
+
+ E_field = atom_field_problem.E_field ;
+ modulation_freq = atom_field_problem.modulation_freq ;
+ freq_index = atom_field_problem.freq_index ;
+
+ rhoLiouville=rhoLiouville_steady_state(L0m, polarizability_m, E_field, modulation_freq);
+ xi=susceptibility(freq_index, rhoLiouville, dipole_elements, E_field);
+ xi_linear = xi.linear;
+ xi_right = xi.right;
+ xi_left = xi.left;
+endfunction
+
+function [dEdz_coef_linear, dEdz_coef_left, dEdz_coef_right] =dEdz_at_freq(wi, rhoLiouville, dipole_elements, E_field)
+% complex absorption coefficient
+% at given E_field frequency component
+ Nlevels=( size(dipole_elements.linear)(1) );
+ rho=rhoOfFreq(rhoLiouville, wi, Nlevels);
+
+ dEdz_coef_linear = sum( sum( transpose(dipole_elements.linear) .* rho ) );
+ dEdz_coef_left = sum( sum( transpose(dipole_elements.left) .* rho ) );
+ dEdz_coef_right = sum( sum( transpose(dipole_elements.right) .* rho ) );
+
+endfunction
+
+function [dEdz_linear, dEdz_left, dEdz_right]=dEdz( atom_field_problem)
+% complex absorption coefficient for each field
+ global atom_properties;
+ L0m = atom_properties.L0m ;
+ polarizability_m = atom_properties.polarizability_m ;
+ dipole_elements = atom_properties.dipole_elements ;
+
+ E_field = atom_field_problem.E_field ;
+ modulation_freq = atom_field_problem.modulation_freq ;
+ Nfreq=length(modulation_freq);
+
+ rhoLiouville=rhoLiouville_steady_state(L0m, polarizability_m, E_field, modulation_freq);
+ for wi=1:Nfreq
+ [dEdz_linear(wi), dEdz_left(wi), dEdz_right(wi)] = dEdz_at_freq(wi, rhoLiouville, dipole_elements, E_field);
+ endfor
+endfunction
+
+function relative_transmission=total_relative_transmission(atom_field_problem)
+% summed across all frequencies field absorption coefficient
+ global atom_properties;
+ [dEdz_linear, dEdz_left, dEdz_right]=dEdz( atom_field_problem);
+ %dEdz_total= dEdz_linear .* conj(dEdz_linear) ...
+ %+dEdz_left .* conj(dEdz_left) ...
+ %+dEdz_right .* conj(dEdz_right);
+ %total_absorption = sum(dEdz_total);
+
+ %total_absorption = |E|^2 - | E - dEdz | ^2
+ E_field.linear = atom_field_problem.E_field.linear;
+ E_field.right = atom_field_problem.E_field.right;
+ E_field.left = atom_field_problem.E_field.left;
+
+ modulation_freq = atom_field_problem.modulation_freq ;
+ pos_freq_indexes=(modulation_freq>0);
+ % WARNING INTRODUCED HERE BUT MUST BE DEFINE IN ATOM PROPERTIES FILE
+ %n_atoms - proportional to a number of interacting atoms
+ n_atoms=500;
+ transmited_intensities_vector = ...
+ abs(E_field.linear + n_atoms*(1i)*dEdz_linear).^2 ...
+ +abs(E_field.right + n_atoms*(1i)*dEdz_right).^2 ...
+ +abs(E_field.left + n_atoms*(1i)*dEdz_left).^2;
+ transmited_intensities_vector = transmited_intensities_vector(pos_freq_indexes);
+ input_intensities_vector = ...
+ abs(E_field.linear).^2 ...
+ +abs(E_field.right).^2 ...
+ +abs(E_field.left).^2;
+ input_intensities_vector = input_intensities_vector(pos_freq_indexes);
+
+ relative_transmission = sum(transmited_intensities_vector) / sum(input_intensities_vector);
+endfunction
+
+% create full list of atom modulation frequencies from positive light frequency amplitudes
+function [modulation_freq, E_field_amplitudes] = ...
+light_positive_frequencies_and_amplitudes2full_set_of_modulation_frequencies_and_amlitudes(...
+ positive_light_frequencies, positive_light_field_amplitudes )
+ % we should add 0 frequency as first element of our frequencies list
+ modulation_freq = cat(2, 0, positive_light_frequencies, -positive_light_frequencies);
+ % negative frequencies have complex conjugated light fields amplitudes
+ E_field_amplitudes.left = cat(2, 0, positive_light_field_amplitudes.left, conj(positive_light_field_amplitudes.left) );
+ E_field_amplitudes.right = cat(2, 0, positive_light_field_amplitudes.right, conj(positive_light_field_amplitudes.right) );
+ E_field_amplitudes.linear = cat(2, 0, positive_light_field_amplitudes.linear, conj(positive_light_field_amplitudes.linear) );
+endfunction
+
+
+
+% vim: ts=2:sw=2:fdm=indent
generated by cgit v1.2.3 (git 2.39.1) at 2024-11-16 10:51:09 +0000