maxima related code into separate dir

Ignore-this: 460324d098e59fdff6dd25e2ed0aac2a

darcs-hash:20110414025404-067c0-3f766887d6911a2822b4a0b97f9a8423766438b5

3 files changed, 936 insertions, 0 deletions

diff --git a/maxima_code/session.out b/maxima_code/session.out
new file mode 100644
index 0000000..723e4e6
--- /dev/null
+++ b/maxima_code/session.out
@@ -0,0 +1,347 @@
+
+Starts dribbling to session.out (2006/7/5, 21:13:38).
+(%i2) 
+dist(z):=(matrix([1,z],[0,1]))$
+
+(%i3) lens(f):=(matrix([1,0],[-1/f,1]))$
+
+(%i4) 
+abcd(optics_system):=block([tmp:matrix([1,0],[0,1])], (for i:1 thru length(optics_system) do tmp:optics_system[i].tmp ), return(tmp))$ 
+
+prop(q0,abcd):=(abcd[1][1]*q0+abcd[1][2])/(abcd[2][1]*q0+abcd[2][2])$
+
+(%i6) 
+wr2q(w,r):=1/(1/r-%i*lamb/%pi/(w^2))$
+
+(%i7) q2w(q):=sqrt(-lamb/%pi/imagpart(1/q))$
+
+(%i8) q2r(q):=block([rp:(realpart(1/q))],
+	(if (rp = 0)
+		then
+		(print("#dbg: R=1/0, will substitute R with something big"),
+		return(1e100))
+		else
+			return(1/rp)
+	)  
+) $
+
+(%i9) 
+
+/* find q0 for such q */
+q2w0(q):=block([w,r,w0],
+	w:q2w(q),r:q2r(q), 
+	w0:(w/sqrt(1+(%pi*w^2/lamb/r)^2)), 
+	return (w0) 
+	)$
+
+(%i10) q2z0(q):=block([w,r],
+	w:q2w(q),r:q2r(q), return(r/(1+(lamb*r/%pi/w^2)^2))
+	)$
+
+(%i11) q2q0(q):=%i*%pi*q2w0(q)^2/lamb$
+
+(%i12) 
+q0f2z_for_colim_lens(q0,f):=f+%pi*q2w(q0)^2/lamb$
+
+(%i13) 
+qf2colim_lens_posit(qstart,f):= block([q0,z0,zd],
+	q0:q2q0(qstart),
+	z0:q2z0(qstart),
+	zd:q0f2z_for_colim_lens(q0,f),
+	return(zd-z0)
+)$
+
+(%i14) 
+
+check_real(x):=if featurep(x[1],real) then if featurep(x[2],real) then if featurep(x[3],real) then true else false $
+
+(%i15) 
+check_posit(x):=if (part(x[1],2)>0) then if (part(x[2],2)>0) then if (part(x[3],2)>0) then true else false $
+
+(%i16) 
+check_sm1(x):=if (part(x[1],2)<ltotal) then if (part(x[2],2)<ltotal) then if (part(x[3],2)<ltotal) then true else false $
+
+(%i17) 
+select_realsolutions(x):=block([temp:[]], for i:1 thru length(x) do (if check_real(x[i]) then temp:float(append(temp,[x[i]])) ), return(temp) ) $ 
+
+select_posit_solutions(x):=block([temp:[]], for i:1 thru length(x) do (if check_posit(x[i]) then temp:float(append(temp,[x[i]])) ), return(temp) ) $ 
+
+select_sm1_solutions(x):=block([temp:[]], for i:1 thru length(x) do (if check_sm1(x[i]) then temp:float(append(temp,[x[i]])) ), return(temp) ) $ 
+
+/*=========================================================*/
+/* 
+   match one Gauss beam parameter q0(w0,r0)
+   to another one qf(wf,qf)
+   with the help of two known lens.
+   distance between q0 and qf is 'ltotal'
+   q0>|----z1-----(f1)------z2-------(f2)------z3--------|>qf
+      ^                                                  ^
+      |-------------------ltotal-------------------------|
+   the system should looks like above.
+ */  
+/*=========================================================*/
+qffql2z1z2z3(q0,f1,f2,qf,ltotal):=block(
+	[
+	 m1,m2,m3,m4,m5,z1,z2,z3,
+	 optics, mfull, 
+	 q0r,qoi, qfr,qfi,
+	 eqt,eq1,solut,
+	 eq1,eq2,eq3,
+	 eqn1, eqn2,eqn3,
+	 eqnev1, eqnev2, eqnev3,
+	 slt, rsol
+	 ],
+	 
+	m1:dist(z1),
+	m2:lens(f1),
+	m3:dist(z2),
+	m4:lens(f2),
+	m5:dist(z3),
+
+	optics:[m1,m2, m3, m4, m5],
+	mfull:abcd(optics),
+
+	eqt:qf=prop(q0,mfull),
+	
+	q0r:float(realpart(q0)),
+	q0i:float(imagpart(q0)),
+	qfr:float(realpart(qf)),
+	qfi:float(imagpart(qf)),
+	
+	eq1:ev(eqt, q0=q0r+q0i*%i, qf=qfr+qfi*%i),
+
+	solut:solve([eq1], [z3]),
+
+	eqn1:imagpart(solut[1]),
+	eqn2:realpart(solut[1]),
+	eqn3:z1+z2+z3=ltotal,
+
+	eqnev1:ev(eqn1,numer),
+	eqnev2:ev(eqn2,numer),
+	eqnev3:ev(eqn3,numer),
+	
+
+	slt:solve([eqnev1,eqnev2,eqnev3],[z1,z2,z3]),
+	
+	print("slt=",slt),
+
+	rsol:select_realsolutions(slt),
+	rsol:select_posit_solutions(rsol),
+	rsol:select_sm1_solutions(rsol),
+	
+	print("rsol=",rsol),
+
+
+	return(rsol)
+);
+
+(%o20) qffql2z1z2z3(q0, f1, f2, qf, ltotal) := 
+block([m1, m2, m3, m4, m5, z1, z2, z3, optics, mfull, q0r, qoi, qfr, qfi, eqt, 
+eq1, solut, eq1, eq2, eq3, eqn1, eqn2, eqn3, eqnev1, eqnev2, eqnev3, slt, 
+rsol], m1 : dist(z1), m2 : lens(f1), m3 : dist(z2), m4 : lens(f2), 
+m5 : dist(z3), optics : [m1, m2, m3, m4, m5], mfull : abcd(optics), 
+eqt : qf = prop(q0, mfull), q0r : float(realpart(q0)), 
+q0i : float(imagpart(q0)), qfr : float(realpart(qf)), 
+qfi : float(imagpart(qf)), eq1 : ev(eqt, q0 = q0r + q0i %i, 
+qf = qfr + qfi %i), solut : solve([eq1], [z3]), eqn1 : imagpart(solut ), 
+								     1
+eqn2 : realpart(solut ), eqn3 : z1 + z2 + z3 = ltotal, 
+		     1
+eqnev1 : ev(eqn1, numer), eqnev2 : ev(eqn2, numer), eqnev3 : ev(eqn3, numer), 
+slt : solve([eqnev1, eqnev2, eqnev3], [z1, z2, z3]), print("slt=", slt), 
+rsol : select_realsolutions(slt), rsol : select_posit_solutions(rsol), 
+rsol : select_sm1_solutions(rsol), print("rsol=", rsol), return(rsol))
+(%i21) /*=========================================================*/
+/*
+	Zero ewaste to zero waste transformatin with lense
+	q0|------------(f)-----------|qf0
+	  ^             ^            ^
+	  |----z--------|            |
+	  |                          |
+	  |-------------l------------|
+	  we are lookking for z1,l
+*/
+	q0fq02z1(q0,f,qf):=block([f0,z1,w0tmp:q2w(q0),wftmp:q2w(qf)],
+		f0:%pi*w0tmp*wftmp/lamb,
+		z1:f+w0tmp/wftmp*sqrt(f^2-f0^2),
+		return(z1)
+		)$
+
+(%i22) 	q0fq02z2(q0,f,qf):=block([f0,z2,w0tmp:q2w(q0),wftmp:q2w(qf)],
+		f0:%pi*w0tmp*wftmp/lamb,
+		z2:f+wftmp/w0tmp*sqrt(f^2-f0^2),
+		return(z2)
+		)$
+
+(%i23) 	q0fq02l(q0,f,qf):=q0fq02z1(q0,f,qf)+q0fq02z2(q0,f,qf)$
+
+(%i24) /*=========================================================*/
+/*
+   for back propogation we need to change R -> -R
+*/
+q2bq(q):=wr2q(q2w(q),-q2r(q));
+
+(%o24) 		       q2bq(q) := wr2q(q2w(q), - q2r(q))
+(%i25) /*=========================================================*/
+/*=========================================================*/
+/*=========================================================*/
+/*=========================================================*/
+
+
+
+/*=========================================================*/
+/* lets solve final waste to parrallel beam, 
+   (actually back propogation)
+*/
+/*=========================================================*/
+/* 
+   this code tryes to find configuration of 3 given lenses with focal
+   length f1,f2, f_col 
+   in such a way that starting from with initial waste (w0) 
+   and radius (r0)
+   after propogation through lens 1 (f1 focal length) 
+   and 2 (f2 focal length) 
+   we have a collimated region with length (dist_betw_f2_f_col)
+   and then after passing through the lense 3 we at some distance
+   match to the final waste and radius (wend, rend).
+   Total distance between starting and final lense position will be
+   calculated in to the 'ltotal'
+   Beam have a wave length 'lamb'.
+   Everything calculated in meters.
+*/   
+/*======== Initial condition ==============================*/
+ltotal:.85$
+
+(%i26) lamb:1.064e-6$
+
+(%i27) /* starting waste */
+w0:0.675e-3$
+
+(%i28) r0:1e100$
+
+(%i29) /* final waste */	
+wend:0.023e-3$
+
+(%i30) rend:1e100$
+
+(%i31) 
+/* lenses focal length */
+f1:.025$
+
+(%i32) f2:.055$
+
+(%i33) f_col:0.07118$
+
+(%i34) 
+/* collimated region size */
+dist_betw_f2_f_col:.30$
+
+(%i35) /*=========================================================*/
+
+optics_after_parral_beam:[lens(-.111),dist(.05)]$
+
+(%i36) 
+	q0:ev(wr2q(w0,r0),numer)$
+
+(%i37) 	qend:wr2q(wend,rend)$
+
+(%i38) /*
+qf:prop(qend,abcd(reverse(optics_after_parral_beam)))$
+float(realpart(qf))$
+float(imagpart(qf))$
+*/
+qf:qend$
+
+(%i39) 
+/*finally position of collimated lens */
+z_col:float(qf2colim_lens_posit(qf,f_col))$
+
+(%i40) if (z_col <0) 
+	then (print("it is not possible to get collimated beam  z_col=",z_col),
+	quit() 
+) $
+
+(%i41) qf:prop(qf,abcd([dist(z_col),lens(f_col)]))$
+
+(%i42) 
+
+
+qf:prop(qf,abcd([dist(dist_betw_f2_f_col),lens(f2)]) )$
+
+(%i43) 
+/* below something wrong */
+dist2waste:ev(-q2z0(qf),numer)$
+
+(%i44) print( "===============");d
+=============== 
+(%o44) 			        ===============
+ist2waste;
+(%o45) 			       0.05407340741495
+(%i46) rad_now:ev(q2r(qf),numer)$
+
+(%i47) qf:ev(q2q0(qf),numer)$
+
+(%i48) /*===========================================*/
+
+
+space_between_wastes:ev(q0fq02l(q0,f1,qf), numer)$
+
+(%i49) dist2first_lens:ev(q0fq02z1(q0,f1,qf),numer)$
+
+(%i50) dist2first_lens_tmp:ev(q0fq02z2(q0,f1,qf),numer)$
+
+(%i51) 
+z1:ev(dist2first_lens,numer)$
+
+(%i52) z2:ev(space_between_wastes+dist2waste,numer)$
+
+(%i53) z3:ev(z2+dist_betw_f2_f_col,numer)$
+
+(%i54) 
+ltotal:z3+z_col$
+
+(%i55) 
+/*
+qf:wr2q( -q2r(qf), q2w(qf) )$
+*/ 
+
+
+/*
+rsol:qffql2z1z2z3(q0,f1,f2,qf,ltotal-z_col)$
+	print("rsol=",rsol),
+*/
+
+
+	with_stdout("answ.txt", 
+		print("##########################################"),
+		print( "lambda=",lamb),
+		print( "Ltot=",  ltotal),
+		print( "r0=",  r0),
+		print( "w0=",  w0),
+		print(),
+		print( "lns1.abcd=abcd_lens(", f1 , ")" ),
+		print( "lns1.x=",  z1   ),
+		print( "lns2.abcd=abcd_lens(", f2 , ")" ),
+		print( "lns2.x=", (z2) ),
+		print( "lns3.abcd=abcd_lens(", f_col , ")" ),
+		print( "lns3.x=", (z3) ),
+		print("wf=", wend),
+		print("rf=", rend),
+		/*
+		qfinal:prop(q0,ev(mfull)),
+		wfinal:float(q2w(qfinal)),
+		rfinal:float(q2r(qfinal)),
+		print("#final waste=", wfinal),
+		print("#I was aiming to final waste"),
+		print("#final radius=", rfinal),
+		print("#I was aiming to final radius"),
+		*/
+		print("##########################################")
+	)$
+
+(%i56) 	system("cat answ.txt")$
+
+(%i57) 	system("xterm -e octave abcd.m") $
+
+(%i58) 
+quit()$
diff --git a/maxima_code/w2wl.max b/maxima_code/w2wl.max
new file mode 100644
index 0000000..158c900
--- /dev/null
+++ b/maxima_code/w2wl.max
@@ -0,0 +1,288 @@
+writefile("session.out")$
+
+dist(z):=(MATRIX([1,z],[0,1]))$
+lens(f):=(MATRIX([1,0],[-1/f,1]))$
+
+abcd(optics_system):=block([tmp:MATRIX([1,0],[0,1])], (for i:1 thru length(optics_system) do tmp:optics_system[i].tmp ), return(tmp))$ 
+
+prop(q0,abcd):=(abcd[1][1]*q0+abcd[1][2])/(abcd[2][1]*q0+abcd[2][2])$
+
+wr2q(w,r):=1/(1/r-%I*lamb/%PI/(w^2))$
+q2w(q):=sqrt(-lamb/%PI/imagpart(1/q))$
+q2r(q):=block([rp:(realpart(1/q))],
+	(if (rp = 0)
+		then
+		(print("#dbg: R=1/0, will substitute R with something big"),
+		return(1e100))
+		else
+			return(1/rp)
+	)  
+) $
+
+
+/* find q0 for such q */
+q2w0(q):=block([w,r,w0],
+	w:q2w(q),r:q2r(q), 
+	w0:(w/sqrt(1+(%Pi*w^2/lamb/r)^2)), 
+	return (w0) 
+	)$
+q2z0(q):=block([w,r],
+	w:q2w(q),r:q2r(q), return(r/(1+(lamb*r/%Pi/w^2)^2))
+	)$
+q2q0(q):=%I*%Pi*q2w0(q)^2/lamb$
+
+q0f2z_for_colim_lens(q0,f):=f+%Pi*q2w(q0)^2/lamb$
+
+qf2colim_lens_posit(qstart,f):= block([q0,z0,zd],
+	q0:q2q0(qstart),
+	z0:q2z0(qstart),
+	zd:q0f2z_for_colim_lens(q0,f),
+	return(zd-z0)
+)$
+
+
+check_real(x):=if featurep(x[1],real) then if featurep(x[2],real) then if featurep(x[3],real) then true else false $
+
+check_posit(x):=if (part(x[1],2)>0) then if (part(x[2],2)>0) then if (part(x[3],2)>0) then true else false $
+
+check_sm1(x):=if (part(x[1],2)<ltotal) then if (part(x[2],2)<ltotal) then if (part(x[3],2)<ltotal) then true else false $
+
+select_realsolutions(x):=block([temp:[]], for i:1 thru length(x) do (if check_real(x[i]) then temp:float(append(temp,[x[i]])) ), return(temp) ) $ 
+
+select_posit_solutions(x):=block([temp:[]], for i:1 thru length(x) do (if check_posit(x[i]) then temp:float(append(temp,[x[i]])) ), return(temp) ) $ 
+
+select_sm1_solutions(x):=block([temp:[]], for i:1 thru length(x) do (if check_sm1(x[i]) then temp:float(append(temp,[x[i]])) ), return(temp) ) $ 
+
+/*=========================================================*/
+/* 
+   match one Gauss beam parameter q0(w0,r0)
+   to another one qf(wf,qf)
+   with the help of two known lens.
+   distance between q0 and qf is 'ltotal'
+   q0>|----z1-----(f1)------z2-------(f2)------z3--------|>qf
+      ^                                                  ^
+      |-------------------ltotal-------------------------|
+   the system should looks like above.
+ */  
+/*=========================================================*/
+qffql2z1z2z3(q0,f1,f2,qf,ltotal):=block(
+	[
+	 m1,m2,m3,m4,m5,z1,z2,z3,
+	 optics, mfull, 
+	 q0r,qoi, qfr,qfi,
+	 eqt,eq1,solut,
+	 eq1,eq2,eq3,
+	 eqn1, eqn2,eqn3,
+	 eqnev1, eqnev2, eqnev3,
+	 slt, rsol
+	 ],
+	 
+	m1:dist(z1),
+	m2:lens(f1),
+	m3:dist(z2),
+	m4:lens(f2),
+	m5:dist(z3),
+
+	optics:[m1,m2, m3, m4, m5],
+	mfull:abcd(optics),
+
+	eqt:qf=prop(q0,mfull),
+	
+	q0r:float(realpart(q0)),
+	q0i:float(imagpart(q0)),
+	qfr:float(realpart(qf)),
+	qfi:float(imagpart(qf)),
+	
+	eq1:ev(eqt, q0=q0r+q0i*%I, qf=qfr+qfi*%I),
+
+	solut:solve([eq1], [z3]),
+
+	eqn1:imagpart(solut[1]),
+	eqn2:realpart(solut[1]),
+	eqn3:z1+z2+z3=ltotal,
+
+	eqnev1:ev(eqn1,numer),
+	eqnev2:ev(eqn2,numer),
+	eqnev3:ev(eqn3,numer),
+	
+
+	slt:solve([eqnev1,eqnev2,eqnev3],[z1,z2,z3]),
+	
+	print("slt=",slt),
+
+	rsol:select_realsolutions(slt),
+	rsol:select_posit_solutions(rsol),
+	rsol:select_sm1_solutions(rsol),
+	
+	print("rsol=",rsol),
+
+
+	return(rsol)
+);
+/*=========================================================*/
+/*
+	Zero ewaste to zero waste transformatin with lense
+	q0|------------(f)-----------|qf0
+	  ^             ^            ^
+	  |----z--------|            |
+	  |                          |
+	  |-------------l------------|
+	  we are lookking for z1,l
+*/
+	q0fq02z1(q0,f,qf):=block([f0,z1,w0tmp:q2w(q0),wftmp:q2w(qf)],
+		f0:%Pi*w0tmp*wftmp/lamb,
+		z1:f+w0tmp/wftmp*sqrt(f^2-f0^2),
+		return(z1)
+		)$
+	q0fq02z2(q0,f,qf):=block([f0,z2,w0tmp:q2w(q0),wftmp:q2w(qf)],
+		f0:%Pi*w0tmp*wftmp/lamb,
+		z2:f+wftmp/w0tmp*sqrt(f^2-f0^2),
+		return(z2)
+		)$
+	q0fq02l(q0,f,qf):=q0fq02z1(q0,f,qf)+q0fq02z2(q0,f,qf)$
+/*=========================================================*/
+/*
+   for back propogation we need to change R -> -R
+*/
+q2bq(q):=wr2q(q2w(q),-q2r(q));
+/*=========================================================*/
+/*=========================================================*/
+/*=========================================================*/
+/*=========================================================*/
+
+
+
+/*=========================================================*/
+/* lets solve final waste to parrallel beam, 
+   (actually back propogation)
+*/
+/*=========================================================*/
+
+optics_after_parral_beam:[lens(-.111),dist(.05)]$
+
+ltotal:4.85$
+lamb:1.064e-6$
+/* modecleaner waste */
+w0:25.63e-6$
+r0:1e100$
+	q0:ev(wr2q(w0,r0),numer)$
+
+q:ev(prop(q0,abcd([dist(.05),lens(-.111)])),numer);
+dist2waste:ev(-q2z0(q),numer);
+
+quit();
+
+
+/* OPA waste */	
+wend:37.09e-6$
+rend:1e100$
+	qend:wr2q(wend,rend)$
+/*
+qf:prop(qend,abcd(reverse(optics_after_parral_beam)))$
+float(realpart(qf))$
+float(imagpart(qf))$
+*/
+qf:qend$
+f_col:.203$
+
+/*finally position of collimated lens */
+z_col:float(qf2colim_lens_posit(qf,f_col))$
+if (z_col <0) 
+	then (print("it is not possible to get collimated beam  z_col=",z_col),
+	quit() 
+) $
+qf:prop(qf,abcd([dist(z_col),lens(f_col)]))$
+
+
+
+f1:.075$
+f2:.075$
+
+dist_betw_f2_f_col:.40$
+
+qf:prop(qf,abcd([dist(dist_betw_f2_f_col),lens(f2)]) )$
+
+/* below something wrong */
+dist2waste:ev(-q2z0(qf),numer)$
+rad_now:ev(q2r(qf),numer)$
+qf:ev(q2q0(qf),numer)$
+/*===========================================*/
+
+
+space_between_wastes:ev(q0fq02l(q0,f1,qf), numer)$
+dist2first_lens:ev(q0fq02z1(q0,f1,qf),numer)$
+dist2first_lens_tmp:ev(q0fq02z2(q0,f1,qf),numer)$
+
+z1:ev(dist2first_lens,numer)$
+z2:ev(space_between_wastes+dist2waste,numer)$
+z3:ev(z2+dist_betw_f2_f_col,numer)$
+
+ltotal:z3+z_col$
+
+/*
+qf:wr2q( -q2r(qf), q2w(qf) )$
+*/ 
+
+
+/*
+rsol:qffql2z1z2z3(q0,f1,f2,qf,ltotal-z_col)$
+	print("rsol=",rsol),
+*/
+
+
+	with_stdout("answ.txt", 
+		print("##########################################"),
+		print( "lambda=",lamb),
+		print( "Ltot=",  ltotal),
+		print( "r0=",  r0),
+		print( "w0=",  w0),
+		print(),
+		print( "lns1.abcd=abcd_lens(", f1 , ")" ),
+		print( "lns1.x=",  z1   ),
+		print( "lns2.abcd=abcd_lens(", f2 , ")" ),
+		print( "lns2.x=", (z2) ),
+		print( "lns3.abcd=abcd_lens(", f_col , ")" ),
+		print( "lns3.x=", (z3) ),
+		print("wf=", wend),
+		print("rf=", rend),
+		/*
+		qfinal:prop(q0,ev(mfull)),
+		wfinal:float(q2w(qfinal)),
+		rfinal:float(q2r(qfinal)),
+		print("#final waste=", wfinal),
+		print("#I was aiming to final waste"),
+		print("#final radius=", rfinal),
+		print("#I was aiming to final radius"),
+		*/
+		print("##########################################")
+	)$
+	system("cat answ.txt")$
+	system("xterm -e octave abcd.m") $
+
+quit()$
+
+/*
+solve([eqn1,eqn2,eqn3],[z1,z2,z3]);
+*/
+
+
+
+/*
+
+
+z3:ltotal-z1-z2$
+fabcd:ev(mfull)$
+f(z1,z2):=ev(float(abs(qf)-abs(prop(q0,fabcd))))$
+fr(z1,z2):=abs(f(z1,z2));
+fi(z1,z2):=imagpart(f(z1,z2));
+fr(.2,.4);
+fi(.2,.4);
+
+plot3d(f(z1,z2), [z1,0,1],[z2,0,1], ['colour_z,true], [plot_format,gnuplot],
+[run_viewer,false]);
+
+printreal(x):=if featurep(x[1],real) then if featurep(x[2],real) then if featurep(x[3],real) then x;
+
+printreal(slt[4]);float(ev(realpart(solut),number));
+float(ev(imagpart(solut),number));
+*/
diff --git a/maxima_code/waste2waste.max b/maxima_code/waste2waste.max
new file mode 100644
index 0000000..c431352
--- /dev/null
+++ b/maxima_code/waste2waste.max
@@ -0,0 +1,301 @@
+writefile("session.out")$
+
+dist(z):=(matrix([1,z],[0,1]))$
+lens(f):=(matrix([1,0],[-1/f,1]))$
+
+abcd(optics_system):=block([tmp:matrix([1,0],[0,1])], (for i:1 thru length(optics_system) do tmp:optics_system[i].tmp ), return(tmp))$ 
+
+prop(q0,abcd):=(abcd[1][1]*q0+abcd[1][2])/(abcd[2][1]*q0+abcd[2][2])$
+
+wr2q(w,r):=1/(1/r-%i*lamb/%pi/(w^2))$
+q2w(q):=sqrt(-lamb/%pi/imagpart(1/q))$
+q2r(q):=block([rp:(realpart(1/q))],
+	(if (rp = 0)
+		then
+		(print("#dbg: R=1/0, will substitute R with something big"),
+		return(1e100))
+		else
+			return(1/rp)
+	)  
+) $
+
+
+/* find q0 for such q */
+q2w0(q):=block([w,r,w0],
+	w:q2w(q),r:q2r(q), 
+	w0:(w/sqrt(1+(%pi*w^2/lamb/r)^2)), 
+	return (w0) 
+	)$
+q2z0(q):=block([w,r],
+	w:q2w(q),r:q2r(q), return(r/(1+(lamb*r/%pi/w^2)^2))
+	)$
+q2q0(q):=%i*%pi*q2w0(q)^2/lamb$
+
+q0f2z_for_colim_lens(q0,f):=f+%pi*q2w(q0)^2/lamb$
+
+qf2colim_lens_posit(qstart,f):= block([q0,z0,zd],
+	q0:q2q0(qstart),
+	z0:q2z0(qstart),
+	zd:q0f2z_for_colim_lens(q0,f),
+	return(zd-z0)
+)$
+
+
+check_real(x):=if featurep(x[1],real) then if featurep(x[2],real) then if featurep(x[3],real) then true else false $
+
+check_posit(x):=if (part(x[1],2)>0) then if (part(x[2],2)>0) then if (part(x[3],2)>0) then true else false $
+
+check_sm1(x):=if (part(x[1],2)<ltotal) then if (part(x[2],2)<ltotal) then if (part(x[3],2)<ltotal) then true else false $
+
+select_realsolutions(x):=block([temp:[]], for i:1 thru length(x) do (if check_real(x[i]) then temp:float(append(temp,[x[i]])) ), return(temp) ) $ 
+
+select_posit_solutions(x):=block([temp:[]], for i:1 thru length(x) do (if check_posit(x[i]) then temp:float(append(temp,[x[i]])) ), return(temp) ) $ 
+
+select_sm1_solutions(x):=block([temp:[]], for i:1 thru length(x) do (if check_sm1(x[i]) then temp:float(append(temp,[x[i]])) ), return(temp) ) $ 
+
+/*=========================================================*/
+/* 
+   match one Gauss beam parameter q0(w0,r0)
+   to another one qf(wf,qf)
+   with the help of two known lens.
+   distance between q0 and qf is 'ltotal'
+   q0>|----z1-----(f1)------z2-------(f2)------z3--------|>qf
+      ^                                                  ^
+      |-------------------ltotal-------------------------|
+   the system should looks like above.
+ */  
+/*=========================================================*/
+qffql2z1z2z3(q0,f1,f2,qf,ltotal):=block(
+	[
+	 m1,m2,m3,m4,m5,z1,z2,z3,
+	 optics, mfull, 
+	 q0r,qoi, qfr,qfi,
+	 eqt,eq1,solut,
+	 eq1,eq2,eq3,
+	 eqn1, eqn2,eqn3,
+	 eqnev1, eqnev2, eqnev3,
+	 slt, rsol
+	 ],
+	 
+	m1:dist(z1),
+	m2:lens(f1),
+	m3:dist(z2),
+	m4:lens(f2),
+	m5:dist(z3),
+
+	optics:[m1,m2, m3, m4, m5],
+	mfull:abcd(optics),
+
+	eqt:qf=prop(q0,mfull),
+	
+	q0r:float(realpart(q0)),
+	q0i:float(imagpart(q0)),
+	qfr:float(realpart(qf)),
+	qfi:float(imagpart(qf)),
+	
+	eq1:ev(eqt, q0=q0r+q0i*%i, qf=qfr+qfi*%i),
+
+	solut:solve([eq1], [z3]),
+
+	eqn1:imagpart(solut[1]),
+	eqn2:realpart(solut[1]),
+	eqn3:z1+z2+z3=ltotal,
+
+	eqnev1:ev(eqn1,numer),
+	eqnev2:ev(eqn2,numer),
+	eqnev3:ev(eqn3,numer),
+	
+
+	slt:solve([eqnev1,eqnev2,eqnev3],[z1,z2,z3]),
+	
+	print("slt=",slt),
+
+	rsol:select_realsolutions(slt),
+	rsol:select_posit_solutions(rsol),
+	rsol:select_sm1_solutions(rsol),
+	
+	print("rsol=",rsol),
+
+
+	return(rsol)
+);
+/*=========================================================*/
+/*
+	Zero ewaste to zero waste transformatin with lense
+	q0|------------(f)-----------|qf0
+	  ^             ^            ^
+	  |----z--------|            |
+	  |                          |
+	  |-------------l------------|
+	  we are lookking for z1,l
+*/
+	q0fq02z1(q0,f,qf):=block([f0,z1,w0tmp:q2w(q0),wftmp:q2w(qf)],
+		f0:%pi*w0tmp*wftmp/lamb,
+		z1:f+w0tmp/wftmp*sqrt(f^2-f0^2),
+		return(z1)
+		)$
+	q0fq02z2(q0,f,qf):=block([f0,z2,w0tmp:q2w(q0),wftmp:q2w(qf)],
+		f0:%pi*w0tmp*wftmp/lamb,
+		z2:f+wftmp/w0tmp*sqrt(f^2-f0^2),
+		return(z2)
+		)$
+	q0fq02l(q0,f,qf):=q0fq02z1(q0,f,qf)+q0fq02z2(q0,f,qf)$
+/*=========================================================*/
+/*
+   for back propogation we need to change R -> -R
+*/
+q2bq(q):=wr2q(q2w(q),-q2r(q));
+/*=========================================================*/
+/*=========================================================*/
+/*=========================================================*/
+/*=========================================================*/
+
+
+
+/*=========================================================*/
+/* lets solve final waste to parrallel beam, 
+   (actually back propogation)
+*/
+/*=========================================================*/
+/* 
+   this code tryes to find configuration of 3 given lenses with focal
+   length f1,f2, f_col 
+   in such a way that starting from with initial waste (w0) 
+   and radius (r0)
+   after propogation through lens 1 (f1 focal length) 
+   and 2 (f2 focal length) 
+   we have a collimated region with length (dist_betw_f2_f_col)
+   and then after passing through the lense 3 we at some distance
+   match to the final waste and radius (wend, rend).
+   Total distance between starting and final lense position will be
+   calculated in to the 'ltotal'
+   Beam have a wave length 'lamb'.
+   Everything calculated in meters.
+*/   
+/*======== Initial condition ==============================*/
+ltotal:.85$
+lamb:1.064e-6$
+/* starting waste */
+w0:0.675e-3$
+r0:1e100$
+/* final waste */	
+wend:0.023e-3$
+rend:1e100$
+
+/* lenses focal length */
+f1:.025$
+f2:.055$
+f_col:0.07118$
+
+/* collimated region size */
+dist_betw_f2_f_col:.30$
+/*=========================================================*/
+
+optics_after_parral_beam:[lens(-.111),dist(.05)]$
+
+	q0:ev(wr2q(w0,r0),numer)$
+	qend:wr2q(wend,rend)$
+/*
+qf:prop(qend,abcd(reverse(optics_after_parral_beam)))$
+float(realpart(qf))$
+float(imagpart(qf))$
+*/
+qf:qend$
+
+/*finally position of collimated lens */
+z_col:float(qf2colim_lens_posit(qf,f_col))$
+if (z_col <0) 
+	then (print("it is not possible to get collimated beam  z_col=",z_col),
+	quit() 
+) $
+qf:prop(qf,abcd([dist(z_col),lens(f_col)]))$
+
+
+
+qf:prop(qf,abcd([dist(dist_betw_f2_f_col),lens(f2)]) )$
+
+/* below something wrong */
+dist2waste:ev(-q2z0(qf),numer)$
+print( "===============");dist2waste;
+rad_now:ev(q2r(qf),numer)$
+qf:ev(q2q0(qf),numer)$
+/*===========================================*/
+
+
+space_between_wastes:ev(q0fq02l(q0,f1,qf), numer)$
+dist2first_lens:ev(q0fq02z1(q0,f1,qf),numer)$
+dist2first_lens_tmp:ev(q0fq02z2(q0,f1,qf),numer)$
+
+z1:ev(dist2first_lens,numer)$
+z2:ev(space_between_wastes+dist2waste,numer)$
+z3:ev(z2+dist_betw_f2_f_col,numer)$
+
+ltotal:z3+z_col$
+
+/*
+qf:wr2q( -q2r(qf), q2w(qf) )$
+*/ 
+
+
+/*
+rsol:qffql2z1z2z3(q0,f1,f2,qf,ltotal-z_col)$
+	print("rsol=",rsol),
+*/
+
+
+	with_stdout("answ.txt", 
+		print("##########################################"),
+		print( "lambda=",lamb),
+		print( "Ltot=",  ltotal),
+		print( "r0=",  r0),
+		print( "w0=",  w0),
+		print(),
+		print( "lns1.abcd=abcd_lens(", f1 , ")" ),
+		print( "lns1.x=",  z1   ),
+		print( "lns2.abcd=abcd_lens(", f2 , ")" ),
+		print( "lns2.x=", (z2) ),
+		print( "lns3.abcd=abcd_lens(", f_col , ")" ),
+		print( "lns3.x=", (z3) ),
+		print("wf=", wend),
+		print("rf=", rend),
+		/*
+		qfinal:prop(q0,ev(mfull)),
+		wfinal:float(q2w(qfinal)),
+		rfinal:float(q2r(qfinal)),
+		print("#final waste=", wfinal),
+		print("#I was aiming to final waste"),
+		print("#final radius=", rfinal),
+		print("#I was aiming to final radius"),
+		*/
+		print("##########################################")
+	)$
+	system("cat answ.txt")$
+	system("xterm -e octave propagation.m") $
+
+quit()$
+
+/*
+solve([eqn1,eqn2,eqn3],[z1,z2,z3]);
+*/
+
+
+
+/*
+
+
+z3:ltotal-z1-z2$
+fabcd:ev(mfull)$
+f(z1,z2):=ev(float(abs(qf)-abs(prop(q0,fabcd))))$
+fr(z1,z2):=abs(f(z1,z2));
+fi(z1,z2):=imagpart(f(z1,z2));
+fr(.2,.4);
+fi(.2,.4);
+
+plot3d(f(z1,z2), [z1,0,1],[z2,0,1], ['colour_z,true], [plot_format,gnuplot],
+[run_viewer,false]);
+
+printreal(x):=if featurep(x[1],real) then if featurep(x[2],real) then if featurep(x[3],real) then x;
+
+printreal(slt[4]);float(ev(realpart(solut),number));
+float(ev(imagpart(solut),number));
+*/