% 考慮像散的五鏡腔(SESAM看成一個鏡子)
% 共需三個程序才能進行腔的完整分析，分別是ADhalf_l3_5mirror.m，w_l3.m和本程序
% 改變程序中的l2，得到不同的ADhalf、wsesam、wcrytal，再有另兩個程序分析可得

%有兩個凹面鏡，它們都有子午和弧矢兩個傳輸矩陣
%設R1和R2的離軸角分別為a和b

% 在本程序中，腔由一個入射平面腔鏡，晶體（不考慮熱透鏡效應），一個平凹腔鏡（r＝50mm），一個反射用的平面鏡，一個雙凸和SESAM構成
% 晶體的厚度為2mm，先不可忽略，若晶體的熱透鏡焦距為f，
% 在我的實驗的腔中存在兩個個束腰，即平面鏡處和靠近SESAM處
% 起點從晶體中間開始和SESAM算起
clear;
a=10*pi/180;      %離軸角a為15度，化為弧度
b=5*pi/180;      %離軸角b為10度，化為弧度
r1=100;           % r1為第一個凹鏡的曲率半徑
r2=50;           % r2為第二個凹鏡的曲率半徑
l1=1;            % l1為入射全反鏡與晶體表面之間的距離，在我的實驗中為1mm
l2=52;           % l2為晶體表面與平凹鏡之間的距離，在我的實驗中為24mm,激光焦點恰好在晶體上
syms l3;         % l3為平凹鏡離雙凸鏡之間的距離（中間經過反射用的平鏡）
l5=25;           % l5為雙凸鏡與SESAM之間的距離
lm=2;            % lm為晶體的長度，在我的實驗中為2mm
wl=1.033e-3;     % 激光波長為1.033um
n=1.74;          % Yb:GdCOB的折射率

T1=[1,l5
    0,1];        %T1為SESAM到R2的第一個傳輸矩陣
T2h=[1,0
    -2/r2*cos(b),1];
T2v=[1,0
    -2/r2/cos(b),1];
T3=[1,l3
    0,1];
T4h=[1,0
    -2/r1*cos(a),1];
T4v=[1,0
    -2/r1/cos(a),1];
T5=[1,l2
    0,1];
T6=[1,lm/n
    0,1];
T7=[1,l1
    0,1];
T8=[1,lm/2/n
    0,1];           % T6是晶體一半厚度的傳輸矩陣

Mh=T1*T2h*T3*T4h*T5*T6*T7*T7*T6*T5*T4h*T3*T2h*T1;
ADhalfh=(Mh(1,1)+Mh(2,2))/2
wsesamh=sqrt(wl/pi*abs(Mh(1,2))/sqrt(-1*Mh(1,2)*Mh(2,1)-((Mh(1,1)-Mh(2,2))/2)^2));

Mch=T8*T7*T7*T6*T5*T4h*T3*T2h*T1*T1*T2h*T3*T4h*T5*T8;
ADhalfch=(Mch(1,1)+Mch(2,2))/2;
wcrytalh=sqrt(wl/pi*abs(Mch(1,2))/sqrt(-1*Mch(1,2)*Mch(2,1)-((Mch(1,1)-Mch(2,2))/2)^2));

Mv=T1*T2v*T3*T4v*T5*T6*T7*T7*T6*T5*T4v*T3*T2v*T1;
ADhalfv=(Mv(1,1)+Mv(2,2))/2
wsesamv=sqrt(wl/pi*abs(Mv(1,2))/sqrt(-1*Mv(1,2)*Mv(2,1)-((Mv(1,1)-Mv(2,2))/2)^2));

Mcv=T8*T7*T7*T6*T5*T4v*T3*T2v*T1*T1*T2v*T3*T4v*T5*T8;
ADhalfcv=(Mcv(1,1)+Mcv(2,2))/2;
wcrytalv=sqrt(wl/pi*abs(Mcv(1,2))/sqrt(-1*Mcv(1,2)*Mcv(2,1)-((Mcv(1,1)-Mcv(2,2))/2)^2));
%此時可算出ADhalf=103/87-8/2175*l3
%wsesam=25/4294967296*12131099559561230^(1/2)/(-118750+2575*l3-4*l3^2)^(1/4)
%wcrytal=5/8589934592*527702830840913505^(1/2)*(abs(44436/21025+5846/1576875*l3)/(-118750+2575*l3-4*l3^2)^(1/2))^(1/2)

%再由程序ADhalf_l3_5mirror.m可畫出ADhalf和l3的關系圖，從而找出穩區。
%由程序w_l3.m可畫出SESAM和晶體上的光斑大小與l3的關系圖