% CavityRingdown.m
% Eqn 41 of 
% "Doppler-induced dynamics of fields in Fabry–Perot
% cavities with suspended mirrors", Malik Rakhmanov (2000).
% http://www.ligo.caltech.edu/docs/P/P000017-A.pdf

clear all

% read in ringdown timeseries:
at = importdata('tek00000.csv');
% time vector
tt = at(:,1)+22.e-6;    % s
% voltage at transmission:
Vt = at(:,2);
ttd = tt*1.e6;          % us

% constants
c = 2.9979e8;   % speed of light in meters
lambda = 1.064e-6;   % laser wavelength, meters
k = 2*pi/lambda;     % laser wavenumber, m^-1
Amp = 1.e0;  % amplitude of ringdown, arb units
% fudge for the Doppler field amplitude
AmpD = 2.;

% 40m mode cleaner:
L = 13.54;   % Cavity length in meters
% mirror power loss, assuming 2 mirrors
Loss = 15e-5;
% mirror velocity, assumed constant
vel = 320.e-6;   % m/s
% delay of ringdown wrt measured time series
Dt = 2;  % us

% power transmissivity of input and output coupling mirrors
T1 = 0.002;
T2 = T1;
R1 = 1-T1-Loss;
R2 = 1-T2-Loss;
r1 = sqrt(R1);
r2 = sqrt(R2);
t1 = sqrt(T1);
t2 = sqrt(T2);

% transit time
T = L/c;

% amplitude folding time:
tau = 2*T/abs(log(r1*r2));

% finesse:
Finesse = pi*sqrt(r1*r2)/(1.-r1*r2);

% time vector for prediction:
dt = 100.e-9;
tmax = 1000.e-6;
t = dt:dt:tmax;
td = t*1.e6-Dt;

% time-dependent velocity:
dmx = 20.e-6;
pwr = 0.2;
%vel = vel*(dmx./min(t,dmx)).^pwr;

% Adiabatic field
C = t1*t2*Amp ./ (1.-r1*r2*exp(-2*i*k*vel.*t));

% Doppler modulation:
D0 = t1*t2*Amp*sqrt(i*pi./(2*k*vel*T)).*exp(i*T./(2*k*vel*tau^2));
D = D0.*exp(-t/tau-i*k*vel.*t.^2/(2*T));

% total field
E = C + AmpD*D;

%E = t1*t2*Amp ./ (1.-r1*r2*exp(-2*i*k*vel*t)) + ...
%    D0*exp(-i*k*vel*t.^2/(2*T));

% voltage on DC photodiode at transmission
V = abs(E).^2;

% Divide out the exponential decay, normalize average:
R = exp(2*(t-Dt*1.e-6)/tau);
tc = find(td > 20 & td < 50);
VRavg = mean(V(tc).*R(tc));
V = V/VRavg;
VR = V.*R;

% do the same, for the raw data:
Rt = exp(2*tt/tau);
ttc = find(ttd > 20 & ttd < 50);
VRtavg = mean(Vt(ttc).*Rt(ttc));
Vt = Vt/VRtavg;
VRt = Vt.*Rt;

% compute the chisq:
DVR = VRt - interp1(td,VR,ttd);
chisq = sum(DVR(ttc).^2)/length(ttc);
dstr = sprintf('%s %5.1f %5.1f %7.2f','vel, Dt, chisq = ',...
    vel(end)*1e6,Dt,chisq);
disp(dstr)

% plots:

figure(1)
semilogy(td,V)
xlabel('time, \mu s')
grid on

figure(2)
plot(td,V)
xlabel('time, \mu s')
grid on

figure(3)
plot(ttd,Vt,'r',td,V,'b')
axis([0 50 0 2])
grid on
xlabel('time, \mu s')
title(dstr)

figure(4)
plot(ttd,VRt,'r',td,VR,'b')
axis([0 50 0 2])
grid on
xlabel('time, \mu s')
title(dstr)

figure(5)
plot(t,vel)
xlabel('time, \mu s')
ylabel('velocity(t)')