% refl=load('TEK00003.CSV');
% trans=load('TEK00004.CSV');
% ramp=load('TEK00005.CSV');
cavity_scan=dlmread('escape_eff_after_M1_better_MM\TEK00005.CSV', ',', 1, 0);
%cavity_blocked_data=dlmread('1064_data_before_m1_swap\TEK00002.CSV', ',', 1, 0);
dark_level_data=dlmread('escape_eff_after_M1_better_MM\TEK00002.CSV', ',', 1, 0);

%for some reason there is a single point which is at index 8982 which has a
%trans value of 0.812, this is not in the CSV file when I look at it so I
%am going to delete it???  Also, in the CSV files I keep having to delete
%rows which have a random charachter in them, doesn't seem like a problem

%cavity_scan(6310,:)=[];

dark_level_refl=median(dark_level_data(:,3));
dark_level_trans=median(dark_level_data(:,2));
%cavity_blocked_refl=median(cavity_blocked_data(:,2))-dark_level_refl;

time=cavity_scan(:,1);
refl=cavity_scan(:,3)-dark_level_refl;
trans=cavity_scan(:,2)-dark_level_trans;
ramp=cavity_scan(:,4);

figure(10)
subplot(3,1,1)
plot(time,trans, '.')%, time_00,pks_00,'r*')
ylim([0,Inf])
subplot(3,1,2)
plot(time, refl)
subplot(3,1,3)
plot(time, ramp)%, time_00, pzt_volts_00, 'r*')

% %get rid of data before ramp turns around
% cut_index=2336;
% time(1:cut_index,:)=[];
% refl(1:cut_index,:)=[];
% trans(1:cut_index,:)=[];
% ramp(1:cut_index,:)=[];
% 
% % %get rid of last FSR since it doesn't include the hole HOM peak
% cut_above_index=find(time>0.0457,1);
% time(cut_above_index:end,:)=[];
% refl(cut_above_index:end,:)=[];
% trans(cut_above_index:end,:)=[];
% ramp(cut_above_index:end,:)=[];

%find locations of 00 peaks in trans data, locs_00 are in units of time on
%scope
[pks_00, locs_00,p_00,w_00]=findpeaks(trans,...
    'MinPeakProminence',0.5,...
    'MinPeakDistance', 100);

time_00=time(locs_00);
pzt_volts_00=ramp(locs_00);

figure(1)
subplot(3,1,1)
plot(time,trans, '.', time_00,pks_00,'r*')
ylim([-Inf,Inf])
subplot(3,1,2)
plot(time, refl)
subplot(3,1,3)
plot(time, ramp, time_00, pzt_volts_00, 'r*')
%% do quadratic fit of time (PZT voltage) to frequency

FSR=853.452; %MHz from T1700104 page 4
frequency_00=FSR*[0:1:length(pks_00)-1]';

[fit_freq, gof_freq]=fit(time_00,frequency_00, 'poly2');
freq_vector=fit_freq(time);

[fit_freq_volts, gof_freq]=fit(pzt_volts_00,frequency_00, 'poly2');

figure(2)
plot(fit_freq,time_00,frequency_00, 'o')
xlabel('Voltage to thorlabs PZT driver [V]')
ylabel('OPO resonant frequency [Hz]')

%% find HOM peaks
[pks_hom,locs_hom]=findpeaks(smooth(trans),...
    'MinPeakHeight',0.01,...
        'MinPeakProminence',0.03,...
    'MinPeakDistance', 50);
delete_index=find(pks_hom>0.2);
pks_hom(delete_index)=[];
locs_hom(delete_index)=[];
frequency_hom=freq_vector(locs_hom);

%find a way to add up all homs in each FSR (must be a better way)
% hom_1=pks_hom(1:4:16);
% hom_2=pks_hom(2:4:16);
% hom_3=pks_hom(3:4:16);
% hom_4=pks_hom(4:4:16);
% hom_content=hom_1;%+hom_2+hom_3+hom_4;

%find refl dips
[pks_refl,locs_refl,p,w]=findpeaks(-1*refl,...
    'MinPeakProminence',0.45,...
    'MinPeakDistance', 50);
frequency_refl=freq_vector(locs_refl);

figure(3)
subplot(2,1,1)
plot(freq_vector,trans, freq_vector, smooth(trans),frequency_00,pks_00,'r*',...
    frequency_hom, pks_hom,'b*')%,...
    %freq_vector(times_off_res), trans(times_off_res))
%legend('IR transmission data', '00 peaks', 'mode mismatch peaks','off resonance times')
ylim([0,1.3])
xlim([-Inf, Inf])
ylabel('M2 transmitted PD [V]')
subplot(2,1,2)
plot(freq_vector, refl, frequency_refl, -1*pks_refl, 'r*')
xlim([-Inf, Inf])
xlabel('Resonant frqeuency [MHz]')
ylabel('M1 Reflected PD [V]')
print('OPO_scan_inject_m1_after_m1_swap.png', '-dpng')

%% get awnsers about mode matching, escape efficiency
%what is higher order mode spacing?
%hom_spacing=((frequency_00-frequency_hom)/FSR)/2;
%this is the mode mismatch peak, so it is the second order hom, divide by 2

%mode matching estimate
mode_matching=pks_00./(pks_00+pks_hom);
mode_mismatch=hom_content./(pks_00+hom_content);
% 
% for ii=1:3
% hom_content=
%find standard error
mean_mode_match=mean(mode_matching)
mm_std_err=std(mode_matching)/sqrt(length(mode_matching))

%find times we were off resonance
times_off_res=find(trans<0.02);
refl_off_resonance=median(refl(times_off_res));
refl_dips=-1.*pks_refl/refl_off_resonance;

% %estimate mode matching by integrating over peak heights
% %for each peak find where you will stop integrating
% 
% for ii=1:length(locs_00)
%     %times off resonsnace is a list of indeces of when the cavity is not on
%     %resoaonce, find the next one after each peak
%     next_off_res_00(ii)=times_off_res(find(times_off_res>locs_00(ii), 1));
%     last_off_res_00(ii)=times_off_res(find(times_off_res<locs_00(ii), 1, 'last'));
%     integrate_00(ii)=trapz(trans(last_off_res_00(ii):next_off_res_00(ii)));
% end
% 
% 
% for ii=1:length(locs_hom)
%     %times off resonsnace is a list of indeces of when the cavity is not on
%     %resoaonce, find the next one after each peak
%     next_off_res_hom(ii)=times_off_res(find(times_off_res>locs_hom(ii), 1));
%     last_off_res_hom(ii)=times_off_res(find(times_off_res<locs_hom(ii), 1, 'last'));
%     integrate_hom(ii)=trapz(trans(last_off_res_hom(ii):next_off_res_hom(ii)));
% end
% 
% smooth_trans=smooth(trans, 10);
% figure(11)
% plot(time,trans, '.', time, smooth_trans, time_00,pks_00,'r*', ...
%     time(next_off_res_00),trans(next_off_res_00),'b*',...
%         time(last_off_res_00),trans(last_off_res_00),'g*',...
%         time(locs_hom),pks_hom,'r*', ...
%     time(next_off_res_hom),trans(next_off_res_hom),'b*',...
%         time(last_off_res_hom),trans(last_off_res_hom),'g*')
% 
% mode_matching_integrated=integrate_00./(integrate_00+integrate_hom);
% mean_mode_match_int=mean(mode_matching_integrated)
% mm_int_std_err=std(mode_matching_integrated)/sqrt(length(mode_matching_integrated))
% 

figure(1)
subplot(3,1,1)
plot(time,trans, '.', time_00,pks_00,'r*')
ylim([0,Inf])
subplot(3,1,2)
plot(time, refl, time(times_off_res), refl(times_off_res))
subplot(3,1,3)
plot(time, ramp, time_00, pzt_volts_00, 'r*')

mean_refl_dips=mean(refl_dips)
dips_std_err=std(refl_dips)/sqrt(length(refl_dips))

%mirror reflectivity from E1800011
% From Fabrice:I believe the LHO Mirrors are #37, 11, 6 , 8 for M1 to M4 respectively. 

r_m2=sqrt(0.99818636);  %range of values is +/-0.000004
r_m1=sqrt(0.87298077); %range is +/-0.0007

%how much of reflected power is due to mode mismatch
P_mismatch=refl_off_resonance.*mode_mismatch;

refl_dips_00=(-1.*pks_refl-P_mismatch)./(refl_off_resonance-P_mismatch);
mean_refl_dips_00=mean(refl_dips_00)
dips_std_err_00=std(refl_dips_00)/sqrt(length(refl_dips_00))

syms rloss
for ii=1:length(refl_dips_00)
    s=solve((((r_m1-rloss)/(1-r_m1*rloss))/((r_m1+rloss)/(1+r_m1*rloss)))^2-refl_dips_00(ii), rloss);
    %find the solution which has Rloss^2 less than 1 (the only physical one)
    indx_high=find(double(s).^2<0.9999);
    indx=indx_high(find(double(s(indx_high)).^2>0.1));
    Rloss(ii)=double(s(indx).^2);
end

mean_Rloss=mean(Rloss)
Rloss_std_err=std(Rloss)/sqrt(length(Rloss))

%Rloss_m2_dip inlcudes front coupler and other cavity losses other than m2
loss_not_M2=1-Rloss./r_m2.^2;
mean_loss_not_M2=mean(loss_not_M2)
loss_not_M2_std_err=std(loss_not_M2)/sqrt(length(loss_not_M2))

%Rloss used for escape efficency should include m2
%make this loss measurement into a escape efficiency
%front coupler transmission 
Tf=1-r_m1^2;
Tloss=1-Rloss'.*r_m2^2;
%eta_esc=Tf./(Tf+Tloss)

eta_esc=log(r_m1^2)./log(r_m1^2.*Rloss);

mean_eta_esc=mean(eta_esc)
eta_esc_std_err=std(eta_esc)/sqrt(length(eta_esc))

best_possible_eta_esc= log(r_m1^2)./log(r_m1^2.*r_m2.^2)


%try using blocked cavity data




