function Meinel(adj);
% Calculate [O3], [OH] for v'=1 to 9, Intensity for v'=1 to 9

% *************** MSIS background data and rates ********************
% Read in data from MSIS background atmosphere and the rate constants using
% program previously calculated from AllA.m

[alt,Temp,N2,O2,O,H,M,NO,CO2]=textread('background.txt','%f %f %f %f %f %f %f %f %f','delimiter',',');
n=1:length(Temp);

O(n)=O(n)/adj;

R9 = dlmread('rates9.txt'); % Temp,k1,k2a,k2b,k3,k4,k5a,k5b,k6
R8 = dlmread('rates8.txt');
R7 = dlmread('rates7.txt');
R6 = dlmread('rates6.txt');
R5 = dlmread('rates5.txt');
R4 = dlmread('rates4.txt');
R3 = dlmread('rates3.txt');
R2 = dlmread('rates2.txt');
R1 = dlmread('rates1.txt');

alpha9=0.44; % Yield of H + O3 reaction going into v' = 9 (Klenerman and Smith 1987)
alpha8=0.28;
alpha7=0.19;
alpha6=0.09;
alpha5=0; % Need to change OH5 
alpha4=0; % Need to change OH4 

% *************** Einstein Coefficents *********************
% Read in data produced by file AllA.m

A9 = dlmread('A9.txt'); % matrix containing Temp,Av'-->v'-1, Av'-->v'-2, Av'-->v'-3, Av'-->v'-4, Av'-->v'-5, Av'
A8 = dlmread('A8.txt');
A7 = dlmread('A7.txt');
A6 = dlmread('A6.txt');
A5 = dlmread('A5.txt');
A4 = dlmread('A4.txt');
A3 = dlmread('A3.txt');
A2 = dlmread('A2.txt');
A1 = dlmread('A1.txt');

% **************** Calculate Concentrations *********************

% Calculate [O3] for use in [OH] calculation
% [O3] = ([O][O2](k2a*[N2]+k2b*[O2]))/(k1[H]+k3[NO]+k6[O])

O3(n)=((O(n).*O2(n).*(R9(n,3).*N2(n)+R9(n,4).*O2(n)))./(R9(n,2).*H(n)+R9(n,5).*NO(n)+R9(n,9).*O(n)));
O3=transpose(O3);

dlmwrite('O3.txt',O3);

% Calculate [OH] to find Iv'v" 

% [OH9]=(alpha*k1*[O3]*[H])/(Av'+ k4[O]+k5a[O2]+k5b[CO2])

OH9(n)=(alpha9.*R9(n,2).*O3(n).*H(n))./(A9(n,7)+R9(n,7).*O2(n)+R9(n,8).*CO2(n)+R9(n,6).*O(n));
OH9=transpose(OH9);

% [OH8]=(alpha*k1*[O3]*[H]+OH(9)*(k5a(9)*O2+K5b(9)*CO2)+A98*OH(9))/(Av'+ k4[O]+k5a[O2]+k5b[CO2])

OH8(n)=(alpha8.*R8(n,2).*O3(n).*H(n)+OH9(n).*(R9(n,7).*O2(n)+R9(n,8).*CO2(n))+OH9(n).*A9(n,2))./(A8(n,7)+R8(n,7).*O2(n)+R8(n,8).*CO2(n)+R8(n,6).*O(n));
OH8=transpose(OH8);

% [OH7]=(alpha*k1*[O3]*[H]+OH(8)*(k5a(8)*O2+k5b(8)*CO2)+A97*OH(9)+A87*OH8)/(Av'+ k4[O]+k5a[O2]+k5b[CO2])

OH7(n)=(alpha7.*R7(n,2).*O3(n).*H(n)+OH8(n).*(R8(n,7).*O2(n)+R8(n,8).*CO2(n))+OH9(n).*A9(n,3)+OH8(n).*A8(n,2))./(A7(n,7)+R7(n,7).*O2(n)+R7(n,8).*CO2(n)+R7(n,6).*O(n));
OH7=transpose(OH7);

% [OH6]=(alpha*k1*[O3]*[H]+OH(7)*(k5a(7)*O2+k5b(7)*CO2)+A96*OH(9)+A86*OH8+A76*OH7)/(Av'+ k4[O]+k5a[O2]+k5b[CO2])

OH6(n)=(alpha6.*R6(n,2).*O3(n).*H(n)+OH7(n).*(R7(n,7).*O2(n)+R7(n,8).*CO2(n))+OH9(n).*A9(n,4)+OH8(n).*A8(n,3)+OH7(n).*A7(n,2))./(A6(n,7)+R6(n,7).*O2(n)+R6(n,8).*CO2(n)+R6(n,6).*O(n));
OH6=transpose(OH6);

% [OH5]=(OH(6)*(k5a(6)*O2+k5b*CO2)+A95*OH9+A85*OH8+A75*OH7+A65*OH6)/(Av'+ k4[O]+k5a[O2]+k5b[CO2])

OH5(n)=(alpha5.*R5(n,2).*O3(n).*H(n)+OH6(n).*(R6(n,7).*O2(n)+R6(n,8).*CO2(n))+OH9(n).*A9(n,5)+OH8(n).*A8(n,4)+OH7(n).*A7(n,3)+OH6.*A6(n,2))./(A5(n,7)+R5(n,7).*O2(n)+R5(n,8).*CO2(n)+R5(n,6).*O(n));
OH5=transpose(OH5);

% OH5(n)=(OH6(n).*(R6(n,7).*O2(n)+R6(n,8).*CO2(n))+OH9(n).*A9(n,5)+OH8(n).*A8(n,4)+OH7(n).*A7(n,3)+OH6.*A6(n,2))./(A5(n,7)+R5(n,7).*O2(n)+R5(n,8).*CO2(n)+R5(n,6).*O(n));
% OH5=transpose(OH5);

% [OH4]=(OH(5)*(k5a(5)*O2+k5b(5)*CO2)+A94*OH9+A84*OH8+A74*OH7+A64*OH6+A54*OH5)/(Av'+ k4[O]+k5a[O2]+k5b[CO2])

OH4(n)=(alpha4.*R4(n,2).*O3(n).*H(n)+OH5(n).*(R5(n,7).*O2(n)+R5(n,8).*CO2(n))+OH9(n).*A9(n,6)+OH8(n).*A8(n,5)+OH7(n).*A7(n,4)+OH6.*A6(n,3)+OH5.*A5(n,2))./(A4(n,6)+R4(n,7).*O2(n)+ R4(n,8).*CO2(n)+R4(n,6).*O(n));
OH4=transpose(OH4);

% OH4(n)=(OH5(n).*(R5(n,7).*O2(n)+R5(n,8).*CO2(n))+OH9(n).*A9(n,6)+OH8(n).*A8(n,5)+OH7(n).*A7(n,4)+OH6.*A6(n,3)+OH5.*A5(n,2))./(A4(n,6)+R4(n,7).*O2(n)+ R4(n,8).*CO2(n)+R4(n,6).*O(n));
% OH4=transpose(OH4);

% [OH3]=(OH(4)*(k5a(4)*O2+k5b(4)*CO2)+A83*OH8+A73*OH7+A63*OH6+A53*OH5+A43*OH4)/(Av'+ k4[O]+k5a[O2]+k5b[CO2])

OH3(n)=(OH4(n).*(R4(n,7).*O2(n)+R4(n,8).*CO2(n))+OH8(n).*A8(n,6)+OH7(n).*A7(n,5)+OH6.*A6(n,4)+OH5.*A5(n,3)+OH4.*A4(n,2))./(A3(n,5)+R3(n,7).*O2(n)+ R3(n,8).*CO2(n)+R3(n,6).*O(n));
OH3=transpose(OH3);

% [OH2]=(OH(3)*(k5a(3)*O2+k5b(3)*CO2)+A72*OH7+A62*OH6+A52*OH5+A42*OH4+A32*OH3)/(Av'+ k4[O]+k5a[O2]+k5b[CO2])

OH2(n)=(OH3(n).*(R3(n,7).*O2(n)+R3(n,8).*CO2(n))+OH7(n).*A7(n,6)+OH6.*A6(n,5)+OH5.*A5(n,4)+OH4.*A4(n,3)+OH3.*A3(n,2))./(A2(n,4)+R2(n,7).*O2(n)+R2(n,8).*CO2(n)+R2(n,6).*O(n));
OH2=transpose(OH2);

% [OH1]=(OH(2)*(k5a(2)*O2+k5n(2)*CO2)+A61*OH6+A51*OH5+A41*OH4+A31*OH3)/(Av'+ k4[O]+k5a[O2]+k5b[CO2])

OH1(n)=(OH2(n).*(R2(n,7).*O2(n)+R2(n,8).*CO2(n))+OH6.*A6(n,6)+OH5.*A5(n,5)+OH4.*A4(n,4)+OH3.*A3(n,3)+OH2(n).*A2(n,2))./(A1(n,3)+R1(n,7).*O2(n)+R1(n,8).*CO2(n)+R1(n,6).*O(n));
OH1=transpose(OH1);

OH=[OH9,OH8,OH7,OH6,OH5,OH4,OH3,OH2,OH1];
dlmwrite('OH.txt',OH);

% ************ Calculate intensities ****************

% Calculate Iv'v" = OHv'*Av'v", Note Delta is the max Dv that can come from
% the upper state

table=[];
Delta=5;
for z=[2:Delta+2];
    I=(OH9(n).*A9(n,z));
    table=[table,I];
end
dlmwrite('I9.txt',table);

table=[];
Delta=5;
for z=[2:Delta+2];
    I=(OH8(n).*A8(n,z));
    table=[table,I];
end
dlmwrite('I8.txt',table);

table=[];
Delta=5;
for z=[2:Delta+2];
    I=(OH7(n).*A7(n,z));
    table=[table,I];
end
dlmwrite('I7.txt',table);

table=[];
Delta=5;
for z=[2:Delta+2];
    I=(OH6(n).*A6(n,z));
    table=[table,I];
end
dlmwrite('I6.txt',table);

table=[];
Delta=5;
for z=[2:Delta+2];
    I=(OH5(n).*A5(n,z));
    table=[table,I];
end
dlmwrite('I5.txt',table);

table=[];
Delta=4;
for z=[2:Delta+2];
    I=(OH4(n).*A4(n,z));
    table=[table,I];
end
dlmwrite('I4.txt',table);

table=[];
Delta=3;
for z=[2:Delta+2];
    I=(OH3(n).*A3(n,z));
    table=[table,I];
end
dlmwrite('I3.txt',table);

table=[];
Delta=2;
for z=[2:Delta+2];
    I=(OH2(n).*A2(n,z));
    table=[table,I];
end
dlmwrite('I2.txt',table);

table=[];
Delta=1;
for z=[2:Delta+2];
    I=(OH1(n).*A1(n,z));
    table=[table,I];
end
dlmwrite('I1.txt',table);
 
%************** Create table of intensities **********

rayleightable;

% produces txt file 'IntTable.txt' which has Intensities for v'-1, v'-2, v'-3, v'-4, v'-5,
% total intensity for each band, from v'=1 at the top to v'=9 at the bottom.
% units are kRayleighs

%*********** Create table of Temperatures *********************

TempCalc;

% *************************************************************************

function[Int] = rayleigh(filename,delta)
% Converts an intensity band into Rayleighs. Integrates over the
% altitudes of the whole band to give a single intensity. Returns total intensity 
% for all v' bands and the intensity of each v' to v" band, in the order v'-5, v'-4 etc. 

% Finding the intensity flux of each transistion

I = dlmread(filename); % Read in volume emission rate for each band

a=1*10^5; % 1km = 10000cm
b=1*10^6; % to change to Rayleighs (10^6*ph/cm2*sec*sphere)

%m=(1:(delta+1));
m=(delta+1):-1:1;
last=length(I); 
n=1:last; % Set up indexing system  
interval=I(n,m).*a;
Int=sum(interval); % integrate over altitude, with intervals of 1km 
Int=Int./b; %convert to Rayleighs
Int=Int./1000; % convert to K Rayleighs

% *************************************************************************

function rayleightable

tab9 = rayleigh('I9.txt',5); % converting intensities to rayleighs
tab8 = rayleigh('I8.txt',5);
tab7 = rayleigh('I7.txt',5);
tab6 = rayleigh('I6.txt',5);
tab5 = rayleigh('I5.txt',5);
tab4 = rayleigh('I4.txt',4);
tab4 = [tab4,0];
tab3 = rayleigh('I3.txt',3);
tab3 = [tab3,0,0];
tab2 = rayleigh('I2.txt',2);
tab2 = [tab2,0,0,0];
tab1 = rayleigh('I1.txt',1);
tab1 = [tab1,0,0,0,0];

table=[tab1;tab2;tab3;tab4;tab5;tab6;tab7;tab8;tab9];
dlmwrite('IntTable.txt',table);


%***********************************************************************

function TempCalc

% *************** MSIS background data and rates ********************
% Read in data from MSIS background atmosphere and the rate constants using program constants.m

[alt,T,N2,O2,O,H,M,NO,CO2]=textread('background.txt','%f %f %f %f %f %f %f %f %f','delimiter',',');
n=1:length(T);

% Read in OH concentrations for the various v' levels

OH = dlmread('OH.txt'); % OH9, OH8, OH7, OH6, OH5, OH4, OH3, OH2, OH1

% Calculate the [OH] over the whole altitude range ie. integral of [OH]
% from alt max to min

a=1*10^5; % 1km = 10000cm
a = a * abs( alt(2)-alt(1) );  % cm in the altitude interval
n=1:(length(OH)); % Set up indexing system  
IntOH= [];

for m=1:9
    interval=OH(n,m).*a;
    Int=sum(interval); % integrate over altitude, with intervals of 1km 
    IntOH = [IntOH,Int];
end

% Calculate the [OH]*T over the whole altitude range ie. integral of [OH]
% from alt max to min

IntOHT= [];

for m=1:9
    interval=OH(n,m).*a.*T;
    Int=sum(interval); % integrate over altitude, with intervals of 1km 
    IntOHT = [IntOHT,Int];
end

% Calculate the average Temperature for each v' level

Temperature= [];

for m=1:9
    avg=IntOHT(m)./IntOH(m);
    Temperature = [Temperature;m,avg];
end

dlmwrite('Temp.txt',Temperature);