%% this calcualte the two different medium
%%thermal fluctuation from sine heat
%% Tara 2010_10_23

%coating compound(SiO2 and Ta2O5) thermal properties
k1 = 2.22; %W/mK
C1 = 1.83e9; %J/m^3K
thermal_co1 = 5.8e-6; %(this is calculated outside the code) 
                       % use the result from Evans' paper
%substrate thermal properties, SiO2
k2 = 1.38;
C2=1.64e9;
thermal_co2 = 0.51e-6; %(substrate, bulk, thermal exp )

%some constants for calculation, not wavelength
lambda = sqrt(k2*C2/(k1*C1));

Q= (lambda +1)/(lambda -1);

%ratio between thermal conductivity and heat capacity
%not thermal expansion coeff
alpha1 = k1/C1;
alpha2 = k2/C2;

%coating thickness, based on 300 ppm transmission
d=4.4e-6;

%stepsize inside the coating, 50 step
stepsize = d/50;
x=[0:stepsize:d];

%stepsize in the substrate, make it 10 times the coaing thickness
y =[0:stepsize:10*d];

%choose frequency of interest
f = 10;
w = 2*pi*f;

QQ= Q*exp(2*d*sqrt(i*w/alpha1));
QI= 1/QQ;


%constant C1, C2, C3 as described in the report
const1 = (1/(1+QI))*sqrt(alpha1/(i*w))...
         *(1/k1);
 
const2 = -(1/(1+QQ))*sqrt(alpha1/(i*w))...
         *(1/k1);    

const3 = const1*(exp(-d*sqrt(i*w/(alpha1))))...
         +const2*(exp(d*sqrt(i*w/(alpha1))))   ;
     
     %TF inside the coating
   TF = const1.*exp(-x*sqrt(i*w/alpha1))...
          +const2.*exp(+x*sqrt(i*w/alpha1));
     
     %TF inside the substrate
   TF2 = const3.*exp(-y*sqrt(i*w/alpha2));

% % %TF inside the coating
% % TF = (1/k1).*sqrt(alpha1/(w*i)).*(...
% %      exp(-x.*sqrt(w*i/alpha1)) / (1 + QI  )...
% %     -exp( x.*sqrt(w*i/alpha1)) / (1 + QQ  ));
% % 
% % 
% % 
% % %TF inside the substrate
% % TF2= (1/k1).*sqrt(alpha1/(w*i)).*(...
% %      exp(-d*sqrt(w*i/alpha1)) / (1 + QI  )...
% %     -exp( d*sqrt(w*i/alpha1)) / (1 + QQ  )).*...
% %       exp(-y.*sqrt(i*w/alpha2));

result(:,1)=real(TF);
result(:,2)=imag(TF);

result2(:,1)=real(TF2);
result2(:,2)=imag(TF2);


%calculate magnitude and phase of the TF
mag= sqrt(result(:,1).^2 + result(:,2).^2);
phi= (180/pi) *atan2(result(:,2),result(:,1));


mag2= sqrt(result2(:,1).^2 + result2(:,2).^2);
phi2= (180/pi) *atan2(result2(:,2),result2(:,1));

%concatenate the results from 1 and 2nd medium
magcat=cat(1,mag,mag2(2:end,:));
phicat=cat(1,phi,phi2(2:end,:));

%%plot magcat() and phicat -> TF.png


%create x axis for plotting
xx= [0:stepsize:stepsize*(length(magcat)-1)];

%create time step, 100 steps in one cycle for ploting
t= [0:1/(100*f):1/f];

%create time matrix for calculating
lengtht=length(t);
tt= ones(length(magcat),1) * t;

%convert vectors of magnitude and phase into matrixes
%for calculating
magcat = magcat*ones(1,lengtht);
phicat = phicat*ones(1,lengtht);

%TF at each timestep will be written in 
%timeresp 
%to convert to T, multiplying by the excitation
timeresp= magcat.*cos(w.*tt-phicat.*pi./180);

%plot time resp -> 

%%Excitation parameter
Pin = 10e-3; %power in 10mW
rG = 292e-6; %spotsize
intensity = Pin/(pi*rG^2);
RIN = 1e-4; %
Finesse = 1e4;
abs = 5e-6;  %absorption coeff 

%excitation amplitude is
Exc = (Pin/(pi*rG^2)) * (Finesse/pi)*abs*RIN;

%real temperature response will be
Temp_resp= Exc*timeresp;


%create the themal expansion coeff at each depth
thermal_vec1 = thermal_co1*ones(length(x),1); 
thermal_vec2 = thermal_co2*ones(length(y)-1,1);
thermal_vec = cat(1,thermal_vec1,thermal_vec2);
%create the matrix of thermal expansion coeff for calculation
thermal_matrix= thermal_vec*ones(1,lengtht);





% dz = SUM[ dT * thermal exp coef * stepsize ]
%    = SUM[ Tem_resp *thermal exp coef * stepsize]  
dz = stepsize*sum(Temp_resp.*thermal_matrix,1);

%pick the maximum amplitude 
max_z=max(dz);



L=0.2035;
lamb=1064*10^-9;
c= 3e8;



%envelope
env = (Pin/(pi*rG^2)) * (Finesse/pi)*abs*RIN...
         *magcat;
       
     % df due to TE, x2 for two mirrors, x2 for pk-pk
     %Thermal expansion-contraction.
df = 4*max_z*c/(lamb*L) 

%%%%%%%%%%
%%%plot%%%
%%%%%%%%%%


% Transfer function of the heat at 10 Hz.
% it's not a bode plot
hold off
hold on
plot(xx,env)
plot(xx,-env)
plot(xx,Temp_resp(:,1),...
     xx,Temp_resp(:,25),...
     xx,Temp_resp(:,51),...
     xx,Temp_resp(:,76))

 %axis([0 20e-6 -1.4e-7 1.4e-7])
 axis tight
 axis([0 20e-6 -1.4e-7 1.4e-7])

 grid on
 %axes1 = axes('Parent',figure1,'FontSize',30);
 
 title('Temp profile in the mirror at different times, for 10Hz',...
       'FontSize',14)
xlabel('Depth[m]',...
       'FontWeight','normal',...
       'Color','black','FontSize',18)
   
%% I look up for this command for 2 hrs to change tick label fontsize   
set(gca,'FontSize',20)
ylabel('Temperature [Kelvin]',...
       'FontWeight','normal',...
       'Color','black','FontSize',18)

legh = legend('Temp at time = 0',...
              'Temp at time = 0.025 s',...
              'Temp at time = 0.05 s',...
              'Temp at time = 0.075 s'...
              );
set(legh,'FontSize',13,'FontWeight','normal');

hold off
print -dpng TF.png


%%%%%%%%%%%%%%%%%%%
%%plot 2, dz vs time
%%%%%%%%%%%%%%%%%%
plot(t,dz)

axis tight
 axis([0 0.1 -3e-18 3e-18])

 grid on
 %axes1 = axes('Parent',figure1,'FontSize',20);
 
 title('surface position change(TE) vs time','FontSize',14)
xlabel('Time [s]',...
       'FontWeight','normal',...
       'Color','black','FontSize',18)
   

set(gca,'FontSize',20)
ylabel('Position [m]',...
       'FontWeight','normal',...
       'Color','black','FontSize',18)



print -dpng dz.png