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()$