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:4.85$
lamb:1.064e-6$
/* starting waste */
w0:0.30e-3$
r0:1e100$
/* final waste */	
wend:37.09e-6$
rend:1e100$

/* lenses focal length */
f1:.025$
f2:.075$
f_col:0.7118$

/* 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)$
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));
*/