function AllA(filename)
%Program to calculate all Einstien A Co-efficents and output as files

% **************** Add directories where files are stored ****************


% Directory where the Einstein coefficients to be used are stored.


% Einstein Co-efficent data from Langhoff(1986) data
    % REF: Stephen R. Langhoff, H-J Werner and P. Rosmus, Theoretical
    % Transistion probabilities for the OH Meinel System, [J. Mol. Spec.,
    % 118, 507-529 (1986)]

addpath C:\Users\espy\Documents\MATLAB\OH_Model\LanghoffData\Delta1; 
addpath C:\Users\espy\Documents\MATLAB\OH_Model\LanghoffData\Delta2;
addpath C:\Users\espy\Documents\MATLAB\OH_Model\LanghoffData\Delta3;
addpath C:\Users\espy\Documents\MATLAB\OH_Model\LanghoffData\Delta4;
addpath C:\Users\espy\Documents\MATLAB\OH_Model\LanghoffData\Delta5;
% Directory where the rate coefficient files are stored
addpath C:\Users\espy\Documents\MATLAB\OH_Model\rates;

% Use function Av to read in the data
Av(filename,9,'NO.txt','rates9.txt','A9.txt','v9to8.txt','v9to7.txt','v9to6.txt','v9to5.txt','v9to4.txt');
Av(filename,8,'NO.txt','rates8.txt','A8.txt','v8to7.txt','v8to6.txt','v8to5.txt','v8to4.txt','v8to3.txt');
Av(filename,7,'NO.txt','rates7.txt','A7.txt','v7to6.txt','v7to5.txt','v7to4.txt','v7to3.txt','v7to2.txt');
Av(filename,6,'NO.txt','rates6.txt','A6.txt','v6to5.txt','v6to4.txt','v6to3.txt','v6to2.txt','v6to1.txt');
Av(filename,5,'NO.txt','rates5.txt','A5.txt','v5to4.txt','v5to3.txt','v5to2.txt','v5to1.txt','v5to0.txt');
Av(filename,4,'NO.txt','rates4.txt','A4.txt','v4to3.txt','v4to2.txt','v4to1.txt','v4to0.txt');
Av(filename,3,'NO.txt','rates3.txt','A3.txt','v3to2.txt','v3to1.txt','v3to0.txt');
Av(filename,2,'NO.txt','rates2.txt','A2.txt','v2to1.txt','v2to0.txt');
A1v(filename,1,'NO.txt','v1to0.txt');


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

function Av(filename1,v,NO,filename2, filename3,varargin);
% Use T and the function EinsteinATemp to calculate data for Av' --> v" and output 
% as filename3. Also calculates rate constant data for v' and outputs as
% filename2. 
%Inputs are: 
%   name of file where backround data stored (MSIS),
%   v' level,
%   NO concentration, 
%   file to store rate coefficients in,
%   file to store Einstein Co-efficents in, 
%   filename for the Transistion probabilities for each transistion. 

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

constants(filename1,v,NO);
[z,Temp,N2,O2,O,H,M,NO,CO2]=textread('background.txt','%f %f %f %f %f %f %f %f %f','delimiter',',');

rate = dlmread('rates.txt');
dlmwrite(filename2,rate);
delete('rates.txt');

% ************* Calc the Einstein A values using EinsteinATemp *********

table=[];
 
 for indices=[1:length(varargin)];
     Av=EinsteinATemp(Temp,char(varargin(indices)));
     Av=transpose(Av);
     table=[table,Av];
 end
  
A = table; % matrix containing Temp,Av'-->v'-1, Av'-->v'-2, Av'-->v'-3, Av'-->v'-4, Av'-->v'-5

%  Calculate Av'=Av'-->v'-1 + Av'-->v'-2 + Av'-->v'-3 + Av'-->v'-4 + Av'-->v'-5

A=transpose(A);
Av=sum(A);
Av=transpose(Av);
A=transpose(A);
Aall=[Temp,A,Av];

dlmwrite(filename3,Aall);

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

function A1v(filename1,v,NOfilename,varargin);
% Use T and the function EinsteinATemp to calculate data for A10
% A1 and output as file A1.txt

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

constants(filename1,v,NOfilename);
[z,Temp,N2,O2,O,H,M,NO,CO2]=textread('background.txt','%f %f %f %f %f %f %f %f %f','delimiter',',');

rate = dlmread('rates.txt');
dlmwrite('rates1.txt',rate);
delete('rates.txt');

table=[];

 for indices=[1:length(varargin)];
     Av=EinsteinATemp(Temp,char(varargin(indices)));
     Av=transpose(Av);
     table=[table,Av];
 end
 
%  dlmwrite('data.txt',table);
%  delete('data.txt');
 
A = table; % matrix containing Temp, A10

% Calculate Av'=A1=A10 

A1=[Temp,table,table];

dlmwrite('A1.txt',A1);

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

function constants(filename,v,NOfilename)
% Calculates/Stores data which is constant between the various vibrational levels, eg. rate
% constants, background atmosphere. Input is:
% MSIS background atmosphere
% v' level
% [NO]
% Output is in the form of .txt files:
% 'background.txt' outputs z,Temp,N2,O2,O,H,M,NO,CO2 concentrations in molecules cm-3 
% 'rates.txt' ouputs Temp,k1,k2a,k2b,k3,k4,k5a, k5b all units cm^3 s-1

% Reactions used in this model:

% H + O3 --> HO* + O2        k1 Ref: Nicolet 1971
%                               update 30/08/06 now using JPL (2006) values
%                               instead of Nicolet
% O + O2 + N2 --> O3 + N2    k2a Ref: IUPAC
% O + O2 + O2 --> O3 + O2    k2b Ref: IUPAC
% NO + O3 --> NO2 + O2       k3 Ref: IUPAC
% OH* + O --> H + O2         k4 Ref: Lopez-Moreno 1987
% OH* + O2 --> *OH + O2      k5a Ref: Dodd (1991) and Copeland et all (1993-1997)
% OH* + CO2 --> *OH + CO2    k5b Ref: as above, see Copeland.xls for
% details
% O + O3 --> 2O2              k6 Ref: JPL


% addpath C:\Matlab6p5\Work\OH_Model_run\Rates
% ****************** BACKGROUND ATMOSPHERE FROM MSIS ******************

% Read in data from MSIS background atmosphere

[z, N2, O2, O, H, Temp, M,massdens]=textread(filename,'%f %f %f %f %f %f %f %f');

% ****************** Set up indexing vector ****************** 

last = length(Temp);
n=(1:last);

% set [NO] 

[alt,NO]=textread(NOfilename,'%f %f');
NO(n)=NO(n).*M(n);

% Set [CO2] - values based on linear plot valid for 70km - 110km (see file
% Copeland.xls for details) y=mx+c y=alt, m = -0.125, c = 115, x=CO2

c=115;
m=-0.125;

CO2(n)=(z(n)-c)/m;

for y=1:last;
    if CO2(y)<0
        CO2(y)=0;
    end;
end;

CO2=transpose(CO2); % CO2 mixing ratio in ppm

ppm=1*10^-6;

CO2(n)=ppm.*M(n).*CO2(n);

% ****************** Rate Constants data  ******************

% H + O3 --> HO + O2            k1
% k1 = 1.4e-10 * exp(-470/T) (JPL 2006)

a=1.4*10^-10;
aa=470;

k1(n)=a.*exp(-aa./Temp(n));
k1=transpose(k1);

% O + O2 + N2 --> O3 + N2     k2a
% k2a = 5.6*10^-34(T/300)^-2.6 cm^6 s-1

b=5.6*10^-34;
c=-2.6;

k2a(n)=b.*((Temp(n)/300).^c);
k2a=transpose(k2a);

% O + O2 + O2 --> O3 + O2     k2b
% k2b = 6.0*10^-34(T/300)^-2.6 cm^6 s-1

d=6.0*10^-34;

k2b(n)=d.*((Temp(n)/300).^c);
k2b=transpose(k2b);

% NO + O3 --> NO2 + O2      k3
% k3 = 1.4*10^-12*exp(-1310/T) cm^3 s-1 

e=1.4*10^-12;
f=1310;

k3(n)=e.*exp(-f./Temp(n));
k3=transpose(k3);

% OH* + O --> H + O2        k4 
% k4 taken from Lopez-Moreno J.J., Rodrigo R., Moreno F., Lopez-Puertas M.,
% Molina A. [Planet. Space Sci. 35 (8), 1029-1038, 1987]
% This equation ignores Temperature change and just perturbs due to v'
% Read rate constant for each level v'=1 to v'=9 from data file 
 
[level,k4]=textread('rateOOH.txt','%f %f','headerlines',1,'delimiter',',');

k4=k4(v);
k4(n)=k4;
k4=transpose(k4);

% OH* + O2 --> OH + O2        k5a  cm^3 s-1
% OH* + CO2 --> OH + CO2        k5b  cm^3 s-1
% Read rate constant for each level v'=1 to v'=9 from data file 
 
[level,k5a,k5b]=textread('rateO2CO2.txt','%f %f %f','headerlines',1,'delimiter',',');

k5a=k5a(v);
k5a(n)=k5a;
k5a=transpose(k5a);

k5b=k5b(v);
k5b(n)=k5b;
k5b=transpose(k5b);

% O + O3 --> 2O2              k6 Ref: JPL
% k6= 8.0*10^-12 exp(-2060/T)

g=8.0*10^-12;
h=2060;

k6(n)=g.*exp(-h./Temp(n));
k6=transpose(k6);

% ****************** OUTPUT TO FILES ****************** 

% Background atmosphere

table=[z,Temp,N2,O2,O,H,M,NO,CO2];
dlmwrite('background.txt',table);

% Rate Constants at level v'

rate=[Temp,k1,k2a,k2b,k3,k4,k5a,k5b,k6];
dlmwrite('rates.txt',rate);

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

function[Tot]= EinsteinA(T,filename);
% A function which takes the inputs (Temperature, filename) and produces a
% Av'v" value corresponding to the Temperature given and the
% information given in the file. The file should show the following: v',
% v", J', P1, Q1, R1, P2, Q2, R2, E1, E2. There are two header lines. 

% Read in data from v'x'to'y'.txt plus the energy levels for v=x (ie v' v'' J' P1 Q1 R1 P2 Q2 R2 E1 E2) 

[vx, vy, J, P1, Q1, R1, P2, Q2, R2, E1, E2]=textread(filename,'%f %f %f %f %f %f %f %f %f %f %f','headerlines',2);

% Set up contstants
% k=1.3806505*10^-23; %Boltzmann's const
k=0.695; % in wavenumbers

% Find the A J' --> J'' Einstein Co-efficients
% AJj(k)=P1(n)+Q1(n)+R1(n)+P2(n)+Q2(n)+R2(n);

% For E1
last=length(P1);
part=0;
n=1;
for y=1:last;
    x(y)=2.*J(n)+1;  % (2J'+1)
    dist(y)=x(y)*exp(-(E1(n)./(k*T))); % (2J'+1)*exp(-Ej'/kT)
    part=part+dist(n); % sum of j' (J'+1)*exp(-Ej'/kT)
    n=n+1;
end
part;

last=length(P1);
n=1;
for y=1:last;
    AJj(y)=P1(n)+Q1(n)+R1(n); % sum j'-->j'' 
    pop(y)=(dist(n)/part).*AJj(n); % (2J'+1)*exp(-Ej'/kT)/ (sum of j' (J'+1)*exp(-Ej'/kT))* (sum J' -- J'')
n=n+1;
end

% For E2
last=length(P2);
part2=0;
n=1;
for y=1:last;
    x2(y)=2.*J(n)+1; 
    dist2(y)=x2(y)*exp(-(E2(n)./(k*T)));
    part2=part2+dist(n);
    n=n+1;
end
part2;

last=length(P1);
n=1;
for y=1:last;
    AJj2(y)=P2(n)+Q2(n)+R2(n);
    pop2(y)=(dist2(n)/part2).*AJj2(n);
n=n+1;
end

% calc Av'v"

AVv=sum(pop);
AVv2=sum(pop2);

Tot=AVv+AVv2;

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

function[AvvT]= EinsteinATemp(Temp,filename);
% A function which produces a vector AvvT which corresponds to Av'v" with 
% a temperature, Temp. Input needed is (Temp, filename) ie the filename where 
% the Einstein Coefficents are stored (see function EinsteinA for format).

% Calculating Av'v" for various temperatures 

last = length(Temp);
n=1;
for y=1:last;
    AvvT(y)=EinsteinA(Temp(n),filename);
    n=n+1;
end