summaryrefslogtreecommitdiff
path: root/abcd.m
blob: ac9b1db5d9319b963a1c5d44b3a628999c6c19c2 (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
function fs_abcd = abcd_free_space( distance)
	fs_abcd=[1, distance; 0,1];
endfunction

function lens_abcd =abcd_lens(focal_distance)
	lens_abcd = [1, 0; -1/focal_distance, 1];
endfunction

function qnew=q_afteer_element(q_old,abcd)
	qnew=(q_old*abcd(1,1)+abcd(1,2))/(q_old*abcd(2,1)+abcd(2,2));
endfunction


function optics = arrange_optics_along_x(optics_unsorted)
% arrange optics in proper order so it x position increases with number
	N=length(optics_unsorted);

	% assign x positions
	x=zeros(1,N);
	for i=1:N
		x(i)=optics_unsorted{i}.x;
	end

	[xs,indx]=sort(x);
	cntr=1;
	for i=indx
		optics{cntr}=optics_unsorted{i};
		cntr=cntr+1;
	end
end

function q = prop_forward(x_pos, q_in, x_in,  optics_elements)
% calculate the 'q' parameter of the Gaussian beam propagating through optical
% 'optics_elements' only in the positive direction  along 'x' axis at points 'x_pos'
%  takes the gaussian beam with initial q_in parameter at x_in
%
% all x_pos must be to the right of x_in
% x_pos must be monotonic!
	if (any(x_pos < x_in))
		error('all beam positions must be to the right of the x_in');
	end

	optics_elements=arrange_optics_along_x(optics_elements);

	% Forward propagation to the right of x_in
	Np=length(x_pos); % number of 'x' points
	Nel=length(optics_elements) ;
	q=0*x_pos; % q vector initialization
	q_last_calc=q_in;
	x_last_calc=x_in; % the furthest calculated point
	for i=1:Np
		x_pos_i=x_pos(i);
		for k=1:length(optics_elements) 
			% iterates through optics_elements to make sure 
			% we take them in account for the beam propagation
			el=optics_elements{k};
			if ( (x_last_calc < el.x) && (el.x <= x_pos_i) )
				abcd=abcd_free_space(el.x-x_last_calc);
				q_last_calc=q_afteer_element(q_last_calc,abcd);
				q_last_calc=q_afteer_element(q_last_calc,el.abcd);
				x_last_calc=el.x;
			endif
		endfor
		if (x_pos_i > x_last_calc);
			abcd=abcd_free_space(x_pos_i-x_last_calc);
			q_last_calc=q_afteer_element(q_last_calc,abcd);
			x_last_calc=x_pos_i;
		endif
		q(i)=q_last_calc;
	endfor
end

function q = prop(x_pos, q_in, x_in,  optics_elements)
% calculate the 'q' parameter of the Gaussian beam propagating through optical
% 'optics_elements' array along 'x' axis at points 'x_pos'
%  takes the gaussian beam with initial q_in parameter at x_in
% x_pos must be monotonic!

	q=0*x_pos; % q vector initialization

	if any(x_pos >= x_in)
		% Forward propagation to the right of x_in
		q(x_pos >= x_in) = prop_forward(x_pos(x_pos>=x_in), q_in, x_in,  optics_elements);
	end

	if any(x_pos < x_in)
		% Backward propagation part the left of x_in
		% do it as forward propagation of the reverse beam
		x_backw=x_pos(x_pos<x_in);
		% now let's reflect the beam with respect to x_in
		% and solve the problem as forward propagating.
		x_backw=x_in-x_backw;
		% now we need to flip x positions
		x_backw=fliplr(x_backw);
		% reflected beam means inverted radius of curvature or real part of q parameter
		q_in_backw = -real(q_in) + 1i*imag(q_in);
		optics_elements_backw=optics_elements;
		% we need to flip all optics elements around x_in as well
		for i=1:length(optics_elements_backw)
			optics_elements_backw{i}.x=x_in-optics_elements_backw{i}.x;
		end

		q_backw = prop_forward(x_backw, q_in_backw, 0,  optics_elements_backw);
		% now we need to flip the radius of curvature again
		q_backw =  -real(q_backw) + 1i*imag(q_backw);

		% final assignment of the backwards propagating beam 
		% which we need to flip back
		q(x_pos<x_in) = fliplr(q_backw);
	end

endfunction

function waste =q2waste(q, lambda)
	for i=1:size(q,2)
		waste(i)=sqrt (-lambda/pi/imag(1/q(i)));
	endfor
endfunction

function radius =q2radius(q, lambda)
	for i=1:size(q,2)
		radius(i)=(1/real(1/q(i)));
	endfor
endfunction

function q=waste_r2q(waste,R,lambda)
	q=1/(1/R-I*lambda/pi/waste/waste);
endfunction