function [x,y] = AVKsimulator(ZA,Esrf,Tsrf,pressure,Temperature,stateH2O,alt,nol)
%AVKsimulation simulates the {hum,delD} proxy state averaging kernels (Type 1 and Type2)
%This version has been written by Matthias Schneider on 24.1.2017

%INPUT to function
    %ZA: satellite observing angle [degree], if not avaiable please set to set to 25.0 when calling the function
    %Esrf: surface emmissivity [unitless]
    %Tsrf: surface temperature [K]
    %pressure: atmospheric pressure profile (given at <nol> altitude levels) [hPa]
    %Temperature: atmospheric temperture profile [K] (given at <nol> altitude levels)
    %stateH2O: atmospheric humidity profile ln[ppmv] (given at <nol> altitude levels)
    %alt: altitudes of the <nol> levels [m]
    %nol: number of levels [unitless]

path_reg='regularisation\'; %the regularisation files must be in this subfolder
r=1;

%Planckfuncion
wn_H2O=1250; %frequency for Planckfuncion (H2O)
wn_HDO=1250; %frequency for Planckfuncion (HDO)
BBradH2O(r)=((1.191042e-5)*wn_H2O^3)/(exp(1.4387752*wn_H2O/Tsrf(r))-1);
BBradHDO(r)=((1.191042e-5)*wn_HDO^3)/(exp(1.4387752*wn_HDO/Tsrf(r))-1);
BBradH2O_T(r)=((1.191042e-5)*wn_H2O^3)/(exp(1.4387752*wn_H2O/(Tsrf(r)+1))-1);
BBradHDO_T(r)=((1.191042e-5)*wn_HDO^3)/(exp(1.4387752*wn_HDO/(Tsrf(r)+1))-1);

for i=1:+1:nol,
   BBradH2Oatm(i)=((1.191042e-5)*wn_H2O^3)/(exp(1.4387752*wn_H2O/Temperature(i))-1);    
   BBradHDOatm(i)=((1.191042e-5)*wn_HDO^3)/(exp(1.4387752*wn_HDO/Temperature(i))-1);    
   for j=1:+1:nol,
       BBradH2Oatm_T(i,j)=((1.191042e-5)*wn_H2O^3)/(exp(1.4387752*wn_H2O/Temperature(i))-1);    
       BBradHDOatm_T(i,j)=((1.191042e-5)*wn_HDO^3)/(exp(1.4387752*wn_HDO/Temperature(i))-1);    
       if i==j, 
           BBradH2Oatm_T(i,j)=((1.191042e-5)*wn_H2O^3)/(exp(1.4387752*wn_H2O/(Temperature(i)+1))-1); 
           BBradHDOatm_T(i,j)=((1.191042e-5)*wn_HDO^3)/(exp(1.4387752*wn_HDO/(Temperature(i)+1))-1); 
       end
   end
end

%Calculation of partial columns and their variations
for i=1:+1:nol,
   for j=1:+1:nol,
       stateH2O_kmch(i,j)=0;
       if i==j, 
           stateH2O_kmch(i,j)=1.0; 
       end
   end
end
for j=1:+1:nol,
   stateH2O_pkm(1:nol,j)=stateH2O(1:nol)'+stateH2O_kmch(1:nol,j);  
end

for i=1:+1:nol-1,
   pcH2O(i)=(0.5*(((pressure(i+1)*(exp(stateH2O(i+1))*28.96/(28.96*1e6+18*exp(stateH2O(i+1)))))/(Temperature(i+1)*1.3807e-19))*1e6)*(alt(i+1)-alt(i)))+(0.5*(((pressure(i)*(exp(stateH2O(i))*28.96/(28.96*1e6+18*exp(stateH2O(i)))))/(Temperature(i)*1.3807e-19))*1e6)*(alt(i+1)-alt(i)));
   pcHDO(i)=(0.5*(((pressure(i+1)*(exp(stateH2O(i+1))*28.96/(28.96*1e6+18*exp(stateH2O(i+1)))))/(Temperature(i+1)*1.3807e-19))*1e6)*(alt(i+1)-alt(i)))+(0.5*(((pressure(i)*(exp(stateH2O(i))*28.96/(28.96*1e6+18*exp(stateH2O(i)))))/(Temperature(i)*1.3807e-19))*1e6)*(alt(i+1)-alt(i)));
   for j=1:+1:nol,
       pcH2O_pkm(i,j)=(0.5*(((pressure(i+1)*(exp(stateH2O_pkm(i+1,j))*28.96/(28.96*1e6+18*exp(stateH2O_pkm(i+1,j)))))/(Temperature(i+1)*1.3807e-19))*1e6)*(alt(i+1)-alt(i)))+(0.5*(((pressure(i)*(exp(stateH2O_pkm(i,j))*28.96/(28.96*1e6+18*exp(stateH2O_pkm(i,j)))))/(Temperature(i)*1.3807e-19))*1e6)*(alt(i+1)-alt(i)));
       pcHDO_pkm(i,j)=(0.5*(((pressure(i+1)*(exp(stateH2O_pkm(i+1,j))*28.96/(28.96*1e6+18*exp(stateH2O_pkm(i+1,j)))))/(Temperature(i+1)*1.3807e-19))*1e6)*(alt(i+1)-alt(i)))+(0.5*(((pressure(i)*(exp(stateH2O_pkm(i,j))*28.96/(28.96*1e6+18*exp(stateH2O_pkm(i,j)))))/(Temperature(i)*1.3807e-19))*1e6)*(alt(i+1)-alt(i)));
   end
end

%***************************************************************************
%****  Calculation of Jacobians for tempertures (atm+srf) and humidity  ****
%****                  Calculations according to Eq. (2)               ****
%***************************************************************************
nstr=19;
sig_comp(1:4*nstr,1:3*nol+1)=0;
for strength=1:+1:nstr,
  for i=1:+1:nol-1,
      tau_H2O(i)=0;
      tau_H2Os1(i)=0;
      tau_H2Os2(i)=0;
      tau_HDO(i)=0;
      for u=1:+1:nol,
          tau_H2O_pkm(i,u)=0;
          tau_H2Os1_pkm(i,u)=0;
          tau_H2Os2_pkm(i,u)=0;
          tau_HDO_pkm(i,u)=0;
      end

      xsection_H2O=2^(strength-1)*1e-31*(1/cos(ZA(r)*(pi()/180))); %strength of nstr H2O lines
      xsection_H2Os1=2^(strength-1)*2e-30*(1/cos(ZA(r)*(pi()/180))); %strength of nstr H2O lines
      xsection_H2Os2=2^(strength-1)*9e-29*(1/cos(ZA(r)*(pi()/180))); %strength of nstr H2O lines
      xsection_HDO=2^(strength-1)*1e-31*(1/cos(ZA(r)*(pi()/180))); %strength of nstr HDO lines
      for j=i:+1:nol-1;
         tau_H2O(i)=tau_H2O(i)+pcH2O(j)*xsection_H2O;
         tau_H2Os1(i)=tau_H2Os1(i)+pcH2O(j)*xsection_H2Os1;
         tau_H2Os2(i)=tau_H2Os2(i)+pcH2O(j)*xsection_H2Os2;
         tau_HDO(i)=tau_HDO(i)+pcHDO(j)*xsection_HDO;
         for u=1:+1:nol,
            tau_H2O_pkm(i,u)=tau_H2O_pkm(i,u)+pcH2O_pkm(j,u)*xsection_H2O;
            tau_H2Os1_pkm(i,u)=tau_H2Os1_pkm(i,u)+pcH2O_pkm(j,u)*xsection_H2Os1;
            tau_H2Os2_pkm(i,u)=tau_H2Os2_pkm(i,u)+pcH2O_pkm(j,u)*xsection_H2Os2;
            tau_HDO_pkm(i,u)=tau_HDO_pkm(i,u)+pcHDO_pkm(j,u)*xsection_HDO;
         end
      end
      dkappa_H2O(i)=pcH2O(i)*xsection_H2O;   
      dkappa_H2Os1(i)=pcH2O(i)*xsection_H2Os1;   
      dkappa_H2Os2(i)=pcH2O(i)*xsection_H2Os2;   
      dkappa_HDO(i)=pcHDO(i)*xsection_HDO;   
      for u=1:+1:nol,
         dkappa_H2O_pkm(i,u)=pcH2O_pkm(i,u)*xsection_H2O;   
         dkappa_H2Os1_pkm(i,u)=pcH2O_pkm(i,u)*xsection_H2Os1;   
         dkappa_H2Os2_pkm(i,u)=pcH2O_pkm(i,u)*xsection_H2Os2;   
         dkappa_HDO_pkm(i,u)=pcHDO_pkm(i,u)*xsection_HDO;   
      end
   end

   SueG1_H2O(r)=Esrf(r)*BBradH2O(r)*exp(-tau_H2O(1));
   SueG1_H2O_T(r)=Esrf(r)*BBradH2O_T(r)*exp(-tau_H2O(1));
   SueG2_H2O(r)=0;
   SueG1_H2Os1(r)=Esrf(r)*BBradH2O(r)*exp(-tau_H2Os1(1));
   SueG1_H2Os1_T(r)=Esrf(r)*BBradH2O_T(r)*exp(-tau_H2Os1(1));
   SueG2_H2Os1(r)=0;
   SueG1_H2Os2(r)=Esrf(r)*BBradH2O(r)*exp(-tau_H2Os2(1));
   SueG1_H2Os2_T(r)=Esrf(r)*BBradH2O_T(r)*exp(-tau_H2Os2(1));
   SueG2_H2Os2(r)=0;
   SueG1_HDO(r)=Esrf(r)*BBradHDO(r)*exp(-tau_HDO(1));
   SueG1_HDO_T(r)=Esrf(r)*BBradHDO_T(r)*exp(-tau_HDO(1));
   SueG2_HDO(r)=0;
   for u=1:+1:nol,
      SueG1_H2O_pkm(r,u)=Esrf(r)*BBradH2O(r)*exp(-tau_H2O_pkm(1,u));
      SueG2_H2O_pkm(r,u)=0;
      SueG2_H2O_Tatm(r,u)=0;
      SueG1_H2Os1_pkm(r,u)=Esrf(r)*BBradH2O(r)*exp(-tau_H2Os1_pkm(1,u));
      SueG2_H2Os1_pkm(r,u)=0;
      SueG2_H2Os1_Tatm(r,u)=0;
      SueG1_H2Os2_pkm(r,u)=Esrf(r)*BBradH2O(r)*exp(-tau_H2Os2_pkm(1,u));
      SueG2_H2Os2_pkm(r,u)=0;
      SueG2_H2Os2_Tatm(r,u)=0;
      SueG1_HDO_pkm(r,u)=Esrf(r)*BBradHDO(r)*exp(-tau_HDO_pkm(1,u));
      SueG2_HDO_pkm(r,u)=0;
      SueG2_HDO_Tatm(r,u)=0;
   end

   for i=1:+1:nol-1,
       SueG2_H2O(r)=SueG2_H2O(r)+dkappa_H2O(i)*0.5*(BBradH2Oatm(i)+BBradH2Oatm(i+1))*exp(-tau_H2O(i));
       SueG2_H2Os1(r)=SueG2_H2Os1(r)+dkappa_H2Os1(i)*0.5*(BBradH2Oatm(i)+BBradH2Oatm(i+1))*exp(-tau_H2Os1(i));
       SueG2_H2Os2(r)=SueG2_H2Os2(r)+dkappa_H2Os2(i)*0.5*(BBradH2Oatm(i)+BBradH2Oatm(i+1))*exp(-tau_H2Os2(i));
       SueG2_HDO(r)=SueG2_HDO(r)+dkappa_HDO(i)*0.5*(BBradHDOatm(i)+BBradHDOatm(i+1))*exp(-tau_HDO(i));
       for u=1:+1:nol,
          SueG2_H2O_pkm(r,u)=SueG2_H2O_pkm(r,u)+dkappa_H2O_pkm(i,u)*0.5*(BBradH2Oatm(i)+BBradH2Oatm(i+1))*exp(-tau_H2O_pkm(i,u));
          SueG2_H2O_Tatm(r,u)=SueG2_H2O_Tatm(r,u)+dkappa_H2O(i)*0.5*(BBradH2Oatm_T(i,u)+BBradH2Oatm_T(i+1,u))*exp(-tau_H2O(i));
          SueG2_H2Os1_pkm(r,u)=SueG2_H2Os1_pkm(r,u)+dkappa_H2Os1_pkm(i,u)*0.5*(BBradH2Oatm(i)+BBradH2Oatm(i+1))*exp(-tau_H2Os1_pkm(i,u));
          SueG2_H2Os1_Tatm(r,u)=SueG2_H2Os1_Tatm(r,u)+dkappa_H2Os1(i)*0.5*(BBradH2Oatm_T(i,u)+BBradH2Oatm_T(i+1,u))*exp(-tau_H2Os1(i));
          SueG2_H2Os2_pkm(r,u)=SueG2_H2Os2_pkm(r,u)+dkappa_H2Os2_pkm(i,u)*0.5*(BBradH2Oatm(i)+BBradH2Oatm(i+1))*exp(-tau_H2Os2_pkm(i,u));
          SueG2_H2Os2_Tatm(r,u)=SueG2_H2Os2_Tatm(r,u)+dkappa_H2Os2(i)*0.5*(BBradH2Oatm_T(i,u)+BBradH2Oatm_T(i+1,u))*exp(-tau_H2Os2(i));
          SueG2_HDO_pkm(r,u)=SueG2_HDO_pkm(r,u)+dkappa_HDO_pkm(i,u)*0.5*(BBradHDOatm(i)+BBradHDOatm(i+1))*exp(-tau_HDO_pkm(i,u));
          SueG2_HDO_Tatm(r,u)=SueG2_HDO_Tatm(r,u)+dkappa_HDO(i)*0.5*(BBradHDOatm_T(i,u)+BBradHDOatm_T(i+1,u))*exp(-tau_HDO(i));
       end
   end
   %Spectral signatures
   SueGtot_H2O(r)=SueG1_H2O(r)+SueG2_H2O(r);           
   SueGtot_H2O_T(r)=SueG1_H2O_T(r)+SueG2_H2O(r);
   SueGtot_H2Os1(r)=SueG1_H2Os1(r)+SueG2_H2Os1(r);           
   SueGtot_H2Os1_T(r)=SueG1_H2Os1_T(r)+SueG2_H2Os1(r);
   SueGtot_H2Os2(r)=SueG1_H2Os2(r)+SueG2_H2Os2(r);           
   SueGtot_H2Os2_T(r)=SueG1_H2Os2_T(r)+SueG2_H2Os2(r);
   SueGtot_HDO(r)=SueG1_HDO(r)+SueG2_HDO(r);           
   SueGtot_HDO_T(r)=SueG1_HDO_T(r)+SueG2_HDO(r);
   %Jacobians
   sig_H2O_T(r)=(SueGtot_H2O_T(r))-(SueGtot_H2O(r));
   sig_H2Os1_T(r)=(SueGtot_H2Os1_T(r))-(SueGtot_H2Os1(r));
   sig_H2Os2_T(r)=(SueGtot_H2Os2_T(r))-(SueGtot_H2Os2(r));
   sig_HDO_T(r)=(SueGtot_HDO_T(r))-(SueGtot_HDO(r));
   for u=1:+1:nol,
      %Spectral signatures
      SueGtot_H2O_pkm(r,u)=SueG1_H2O_pkm(r,u)+SueG2_H2O_pkm(r,u);
      SueGtot_H2O_Tatm(r,u)=SueG1_H2O(r)+SueG2_H2O_Tatm(r,u);           
      SueGtot_H2Os1_pkm(r,u)=SueG1_H2Os1_pkm(r,u)+SueG2_H2Os1_pkm(r,u);
      SueGtot_H2Os1_Tatm(r,u)=SueG1_H2Os1(r)+SueG2_H2Os1_Tatm(r,u);           
      SueGtot_H2Os2_pkm(r,u)=SueG1_H2Os2_pkm(r,u)+SueG2_H2Os2_pkm(r,u);
      SueGtot_H2Os2_Tatm(r,u)=SueG1_H2Os2(r)+SueG2_H2Os2_Tatm(r,u);           
      SueGtot_HDO_pkm(r,u)=SueG1_HDO_pkm(r,u)+SueG2_HDO_pkm(r,u);
      SueGtot_HDO_Tatm(r,u)=SueG1_HDO(r)+SueG2_HDO_Tatm(r,u);           
      %Jacobians
      sig_H2O_pkm(r,u)=(SueGtot_H2O_pkm(r,u))-(SueGtot_H2O(r));
      sig_H2O_Tatm(r,u)=(SueGtot_H2O_Tatm(r,u))-(SueGtot_H2O(r));
      sig_H2Os1_pkm(r,u)=(SueGtot_H2Os1_pkm(r,u))-(SueGtot_H2Os1(r));
      sig_H2Os1_Tatm(r,u)=(SueGtot_H2Os1_Tatm(r,u))-(SueGtot_H2Os1(r));
      sig_H2Os2_pkm(r,u)=(SueGtot_H2Os2_pkm(r,u))-(SueGtot_H2Os2(r));
      sig_H2Os2_Tatm(r,u)=(SueGtot_H2Os2_Tatm(r,u))-(SueGtot_H2Os2(r));
      sig_HDO_pkm(r,u)=(SueGtot_HDO_pkm(r,u))-(SueGtot_HDO(r));
      sig_HDO_Tatm(r,u)=(SueGtot_HDO_Tatm(r,u))-(SueGtot_HDO(r));
   end
   %creation of Jacobian Matrix
   %humidity
   sig_comp(strength,1:nol)=sig_H2O_pkm(r,1:nol);
   sig_comp(strength+nstr,1:nol)=sig_H2Os1_pkm(r,1:nol);
   sig_comp(strength+2*nstr,1:nol)=sig_H2Os2_pkm(r,1:nol);
   sig_comp(strength+3*nstr,nol+1:2*nol)=sig_HDO_pkm(r,1:nol);
   %atm temperature
   sig_comp(strength,2*nol+1:3*nol)=sig_H2O_Tatm(r,1:nol);
   sig_comp(strength+nstr,2*nol+1:3*nol)=sig_H2Os1_Tatm(r,1:nol);
   sig_comp(strength+2*nstr,2*nol+1:3*nol)=sig_H2Os2_Tatm(r,1:nol);
   sig_comp(strength+3*nstr,2*nol+1:3*nol)=sig_HDO_Tatm(r,1:nol);
   %srf temperature
   sig_comp(strength,3*nol+1)=sig_H2O_T(r);
   sig_comp(strength+nstr,3*nol+1)=sig_H2Os1_T(r);
   sig_comp(strength+2*nstr,3*nol+1)=sig_H2Os2_T(r);
   sig_comp(strength+3*nstr,3*nol+1)=sig_HDO_T(r);
end
%*******************************************************************

%***********************************
%****  Setup of Regularisation  ****
%***********************************
%reading of inv(Sa) for water vapour from the retrieval setup
Asurf(r)=alt(1);
dirconst=sprintf('inv_%i', floor(Asurf(r)/10)*10);
if (floor(Asurf(r)/10)*10)<1000,
  dirconst=sprintf('inv_0%i', floor(Asurf(r)/10)*10);
end   
if (floor(Asurf(r)/10)*10)<100,
  dirconst=sprintf('inv_00%i', floor(Asurf(r)/10)*10);
end   
if (floor(Asurf(r)/10)*10)<10,
  dirconst=sprintf('inv_000%i', floor(Asurf(r)/10)*10);
end
correl(1:2*nol,1:2*nol)=0;
fidCON=fopen([path_reg,dirconst,'\correl.inp']);
lineCON = fgetl(fidCON);
lineCON = fgetl(fidCON);
lineCON = fgetl(fidCON);
lineCON = fgetl(fidCON);
for i=1:+1:2*nol,
   lineCON = fgetl(fidCON);
   for j=1:+1:2*nol,
      ACON = sscanf(lineCON, '%f');
      correl(i,j)=ACON(j);
   end
end
fclose(fidCON);
kovar(1:2*nol,1:2*nol)=0;
fidCON=fopen([path_reg,dirconst,'\kovar.inp']);
lineCON = fgetl(fidCON);
lineCON = fgetl(fidCON);
for i=1:+1:nol,
   lineCON = fgetl(fidCON);
   for j=1:+1:nol,
      ACON = sscanf(lineCON, '%f');
      kovar(i,j)=ACON(j);
      kovar(i+nol,j+nol)=ACON(j);
   end
end
fclose(fidCON);
Sam1(1:2*nol,1:2*nol)=kovar(1:2*nol,1:2*nol)+correl(1:2*nol,1:2*nol); %inv(Sa) for atmospheric humidity        
%Calcualtion of inv(Sa) for atmospheric temperatures 
Sam1_T(1:nol,1:nol)=0;
abs_T(1:nol,1:nol)=0;
Oper_T(1:nol-1,1:nol)=0;
for i=1:+1:1,
  abs_T(i,i)=1/(1^2);
  std_T(i)=1/(0.25);
end
for i=2:+1:nol,
  abs_T(i,i)=1/(0.25^2);
  std_T(i)=1/(0.25);
end
Tcorrkm(1)=10;
Tcorrkm(2:nol)=5000;
for i=1:+1:nol-1,
   Oper_T(i,i)=1/sqrt((std_T(i)*std_T(i+1))*(1-(exp(-(((alt(i)-alt(i+1))^2)/(2*Tcorrkm(i)*Tcorrkm(i+1)))))));
   Oper_T(i,i+1)=-Oper_T(i,i);
end
smo_T(1:nol,1:nol)=Oper_T(1:nol-1,1:nol)'*Oper_T(1:nol-1,1:nol);
Sam1_T(1:nol,1:nol)=abs_T(1:nol,1:nol)+smo_T(1:nol,1:nol);  %inv(Sa) for atmospheric temperature              
%***********************************

%***********************************************
%************  Simulation of AVK ***************
%***********************************************
%setup of inv(Sa)
Sam1_comp(1:3*nol+1,1:3*nol+1)=0;
Sam1_comp(1:2*nol,1:2*nol)=Sam1(1:2*nol,1:2*nol);
Sam1_comp(2*nol+1:3*nol,2*nol+1:3*nol)=Sam1_T(1:nol,1:nol); %ref constraint for atm temperature
Sam1_comp(3*nol+1,3*nol+1)=1e-12; %ref: no constraint for surface temperature
%setup of inv(Se) 
noise_value=0.5;
noise(1:4*nstr,1:4*nstr)=0;
for strength=1:+1:4*nstr,
  noise(strength,strength)=1/(noise_value^2);
end
%Calculation according to Eq. (4)
kern_theo_all(1:3*nol+1,1:3*nol+1)=pinv(sig_comp(1:4*nstr,1:3*nol+1)'*noise(1:4*nstr,1:4*nstr)*sig_comp(1:4*nstr,1:3*nol+1)+Sam1_comp(1:3*nol+1,1:3*nol+1),1e-12)*sig_comp(1:4*nstr,1:3*nol+1)'*noise(1:4*nstr,1:4*nstr)*sig_comp(1:4*nstr,1:3*nol+1);
%Transformation to {H2O,delD} proxy state -> Type 1 avk
for i=1:+1:nol,
  for j=1:+1:2*nol,
     Traf(i,j)=0;
     if i==j, Traf(i,j)=1/2; end
     if i+nol==j, Traf(i,j)=1/2; end
  end
end
for i=1:+1:nol,
  for j=1:+1:2*nol,
     Traf(i+nol,j)=0;
     if i==j, Traf(i+nol,j)=-1; end
     if i+nol==j, Traf(i+nol,j)=1; end
  end
end
kern_theo_T1(1:2*nol,1:2*nol)=Traf(1:2*nol,1:2*nol)*kern_theo_all(1:2*nol,1:2*nol)*inv(Traf(1:2*nol,1:2*nol));
%A posteriori correction -> Type 2 avk               
CORR_retr(1:2*nol,1:2*nol)=0; %Matrix C, Eq. (14) aus Schneider et al., 2012)
CORR_retr(1:nol,1:nol)=kern_theo_T1(nol+1:2*nol,nol+1:2*nol);
CORR_retr(1+nol:2*nol,1:nol)=-kern_theo_T1(nol+1:2*nol,1:nol);
for i=1:+1:nol,
  CORR_retr(i,i+nol)=0;
  CORR_retr(i+nol,i+nol)=1;
end
kern_theo_T2(1:2*nol,1:2*nol)=CORR_retr(1:2*nol,1:2*nol)*kern_theo_T1(1:2*nol,1:2*nol);
%***********************************************

%*******************************
%****  OUTPUT from function ****
%*******************************
x=kern_theo_T1(1:2*nol,1:2*nol);
y=kern_theo_T2(1:2*nol,1:2*nol);