

% load the OPD from SELF and MASK and determine the amount of RH power
% required for a given mask power to minimize the change in LG20 scatter
close all
clear all

load('LHO_optimum_soln.mat')


PRHnow = 0.75;
PSELFnow = 30;   % 30W input power

PListIN = 0.05:0.1:0.8;
PRHList = [];
for PCO2mask = PListIN
v = -70:5:70;

%% scale the OPD correctly
% CO2 mask
OPDCO2maskPWatt = CO2OPD/(sum(xbest(1:end-1))*0.01);

% SELF HEATING
PIFO = dPIMC; % 2W to 14W into IMC
OPDSELFPWatt = SELFOPD/PIFO;

% RH OPD
dSdP = -9E-6;  % diopters per Watt

%% get the nominal distortion at 30W
meanV = @(OPD) sum(sum(exp(-2*(X.^2+ Y.^2)/0.054^2).*OPD))/sum(sum(exp(-2*(X.^2+ Y.^2)/0.054^2)));
rmsV = @(OPD) sqrt(sum(sum(exp(-2*(X.^2+ Y.^2)/0.054^2).*(OPD - meanV(OPD)).^2))/sum(sum(exp(-2*(X.^2+ Y.^2)/0.054^2))));
selfPower = 30*40/50;          % PRG drops from 50 to 40 at 30W.
OPDSELF = OPDSELFPWatt*selfPower;  % units are [nm]

r = sqrt(X.^2 + Y.^2);
OPDRH = -dSdP*PRHnow/2*r.^2*1E9;

figure('Position', [100, 100, 1700, 1000], 'DefaultAxesFontSize', 11)
OPDnow = OPDRH+OPDSELF;
subplot(2,3,1)
contourf(X*1E3,Y*1E3,OPDSELF - meanV(OPDSELF),v)
title(['SP OPD: H1ITMY SELF @30W uncorrected. RMS = ' , num2str(round(rmsV(OPDSELF))), 'nm'])

subplot(2,3,2)
contourf(X*1E3,Y*1E3,OPDnow - meanV(OPDnow),v)
title(['SP OPD now: H1ITMY SELF @30W + RH @0.75W. RMS = ', num2str(round(rmsV(OPDnow))), 'nm'] )

%% add the nominal CO2 power
subplot(2,3,3)
OPDCO2 = OPDCO2maskPWatt*PCO2mask;
contourf(X*1E3,Y*1E3,OPDCO2 - meanV(OPDCO2), v)
title(['SP OPD CO2 mask: ', num2str(PCO2mask*1E3), 'mW, RMS = ', num2str(round(rmsV(OPDCO2))), 'nm'])

subplot(2,3,4)
contourf(X*1E3,Y*1E3,OPDCO2+OPDnow - meanV(OPDCO2+OPDnow), v)
title(['OPD Now + CO2 mask: ', num2str(PCO2mask*1E3), 'mW, RMS =' num2str(round(rmsV(OPDCO2+OPDnow))), 'nm'])


LG00 = LGpl(0,0,X,Y,0.054);
LG10 = LGpl(1,0,X,Y,0.054);

subplot(2,3,5)
contourf(X*1E3,Y*1E3,LG10.^2, 20)
title(['LG10 mode'])

% get the overlap into LG10
k = 2*pi/1064;
OLamp = @(M1,M2,OPD) sum(sum(M1.*exp(1i*k*OPD).*conj(M2)))/sqrt(sum(sum(M1.*conj(M1)))*sum(sum(M2.*conj(M2))));


OLnow = sum(sum(LG00.*exp(1i*k*OPDnow).*conj(LG10)))*...
    sum(sum(conj(LG00.*exp(1i*k*OPDnow)).*LG10))/(sum(sum(LG00.*conj(LG00)))*sum(sum(LG10.*conj(LG10))));
OLAnow = OLamp(LG00, LG10, OPDnow);

% determine the overlap with CO2 and RH
OPD2 = OPDnow + OPDCO2;
diffOL = [];
RHList = [];
for PRH = 0:0.1:10
    r = sqrt(X.^2 + Y.^2);
    OPDRH = -dSdP*PRH/2*r.^2*1E9;
    
    OPDtest = OPD2 + OPDRH;
    OLAfix = OLamp(LG00, LG10, OPDtest);
    
    diffOL = [diffOL; abs(OLAfix - OLAnow)]; 
    RHList = [RHList; PRH];
end

ind1 = find(diffOL == min(diffOL));
PRHfix = RHList(ind1);

subplot(2,3,6)
OPDRH = -dSdP*PRHfix/2*r.^2*1E9;
OPDfix = OPD2+OPDRH;
contourf(X*1E3,Y*1E3,OPDfix - meanV(OPDfix), v)
title(['OPD Now + CO2 mask: ', num2str(PCO2mask*1E3), 'mW + RH @ ', num2str(PRHnow + PRHfix), 'W total. RMS = ', ...
    num2str(round(rmsV(OPDfix))), 'nm'])

subplot(2,3,5)
contourf(X*1E3,Y*1E3,OPDRH - meanV(OPDRH), v)
title(['Optimum additional RH: ', num2str(PRHfix), 'W. RMS = ', ...
    num2str(round(rmsV(OPDRH))), 'nm'])

iList = [1,2,3,4,5,6];
for ii = iList
    subplot(2,3,ii)
    xlabel('X-coordinate (mm)')
    ylabel('Y-coordinate (mm)')
    foo = colorbar;
    ylabel(foo, 'OPD (nm)')
    grid on
    axis square
    caxis([min(v) max(v)])
end

PRHList = [PRHList; PRHfix rmsV(OPDfix)];
set(gcf, 'PaperSize', [11 17]*1.5)
orient landscape
subplot(2,3,4)
text(0.01, 1.2, ['Contour spacing = ', num2str(abs(v(2)-v(1))), 'nm'], ...
    'Units', 'normalized', 'FontSize', 12)

print('-dpdf', ['results/', mfilename(), '_CO2_power_', num2str(PCO2mask*1E3), 'mW.pdf'])

pause(0.1)
end



figure; plot(PListIN, PRHList(:, 2), '-o'); grid on


