summaryrefslogtreecommitdiff
path: root/psr/psr.m
blob: 563179107f8747db930f3fab5cc4745b96764a40 (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
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=-200.0;
detuning_p_max=-detuning_p_min;
detuning_p_max=1000;
detuning_freq=linspace(detuning_p_min,detuning_p_max,N_detun_steps);
%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=-5/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);
	% we add required ellipticity
	E_field_lab_pos_freq.x*=exp(I*psi_el);
	E_field_lab_pos_freq.y*=exp(-I*psi_el);
	E_field_lab_pos_freq.z=E_field_zero;

	E_field_pos_freq=xyz_lin2atomic_axis_polarization(theta, E_field_lab_pos_freq);


fprintf (stderr, "tuning laser in forloop to set conditions vs detuning\n");
fflush (stderr);
for detuning_p_cntr=1:length(detuning_freq);
	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);
	detuning_p=detuning_freq(detuning_p_cntr);
	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
[xi_linear, xi_left, xi_right]=parcellfun(2, @susceptibility_steady_state_at_freq, problems_cell_array);
%[xi_linear, xi_left, xi_right]=cellfun( @susceptibility_steady_state_at_freq, problems_cell_array);

%save '/tmp/relative_transmission_vs_detuning.mat' detuning_freq relative_transmission_vs_detuning;
save '/tmp/xi_vs_detuning.mat' detuning_freq xi_linear xi_left xi_right  E_field_pos_freq E_field_probe ;

output_psr_results_vs_detuning;

elapsed_time = etime (clock (), t0)

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