
close all

for kk = 1:2
    figure(1)

    % static Lens
    if kk == 1
        ITMstat = -1/80E3;
        ITMnameSTR = 'LHO ITMX';
    else
        ITMstat = 1.7E-6;
        ITMnameSTR = 'LHO ITMY';

    end

    % work out the TCS settings for O2
    RH_SUBdef = -9E-6;
    CO2_SUBdef = 2.5E-5; %was 6.23E-5
    SelfSUBdef = 4.87E-4;

    RH_SURFdef = 1.07E-6;
    SelfSURFdef = -3.60E-5;

    % arm Power gain from input laser power
    armPowerG = 1.20E5/24;

    % absorption
    absNom = 5E-7;
    nomPower = 20;
    %nomLens = absNom*armPowerG*nomPower*SelfSUBdef;
    nomLens = 1/50E3;
    %nomSurf = absNom*armPowerG*nomPower*SelfSURFdef;

    % use a range of absorptions to account for uncertainty in this level
    for absAct = 2.0E-7:0.5E-8:3E-7

        if absAct == 2.5E-7
            LW = 3;
        else
            LW = 1;
        end

        % select a range of input power levels
        pIn = 1:1:60;
        selfLens = SelfSUBdef*(absAct*armPowerG*pIn);
        staticLens = ITMstat*(pIn*0 +1);
        designLens = nomLens*(pIn*0 +1);

        totalLensErr = designLens - (selfLens+staticLens);
        % get the CO2 and RH powers
        CO2power = totalLensErr/CO2_SUBdef;
        CO2power = CO2power.*(CO2power>0);
        RHpower = totalLensErr/RH_SUBdef;
        RHpower = RHpower.*(RHpower>0);

        % get the CO2 and RH lenses
        CO2Lens = CO2power*CO2_SUBdef;
        RHLens = RHpower*RH_SUBdef;


        % plot all the results
        subplot(3,2,[0,2]+kk)
        plot(pIn, (selfLens+ staticLens)*1E6, 'k', 'LineWidth', LW);
        legSTR{1} = 'Self absorption + static lens';
        hold all
        plot(pIn, CO2Lens*1E6,  'r', 'LineWidth', LW);
        legSTR{2} = 'CO2 Lens';
        plot(pIn, RHLens*1E6,  'b', 'LineWidth', LW);
        legSTR{3} = 'RH Lens';
        plot(pIn, designLens*1E6, 'g--', 'LineWidth', 2);
        legSTR{4} = '50km Optimum Les';
        grid on
        xlabel('Input Power (W)')
        ylabel('Substrate Lens (uD)')
        title([ITMnameSTR, ' single-pass substrate lens'])
        axis([0 60 -30 50])
        plot([50, 50], [-30, 50], 'k--', 'LineWidth', 2);
        legSTR{5} = 'O2 power';
        legend(legSTR)


        subplot(3,2,4+kk)
        plot(pIn, CO2power, 'r','LineWidth', LW);
        hold all
        plot(pIn, RHpower, 'b','LineWidth', LW);
        grid on
        xlabel('Input Power (W)')
        ylabel('Actuator Power (W)')
        legend('CO2 central heating', 'RH power')
        title('TCS actuator power levels')
        axis([0 60 0 2])
        plot([50, 50], [0, 2], 'k--', 'LineWidth', 2);


    end

end

orient tall 
print('-dpdf', [mfilename(), '_LHO.pdf'])

