CalcTheoOp

 calculated theoretical zenith opacity based on weather station data

0001 function [TheoOp] = CalcTheoOp(d)
0002 % calculated theoretical zenith opacity based on weather station data
0003 
0004 To = mean(d.array.weather.airTemperature) + 273.15; % atmospheric temp
0005 T = To;
0006 Po = mean(d.array.weather.pressure); % atmospheric pressure
0007 Ppascals = Po*100;
0008 
0009 RelH = mean(d.array.weather.relativeHumidity);
0010 RelH = RelH/100; %get it out of % form
0011 
0012 L = zeros(1,94);
0013 
0014 for  j = 1222:200:20022
0015    
0016      H = j;
0017      V = 5;
0018 
0019 
0020 if    H < 11000
0021       % Temeprature decreasing via lapse rate
0022       T = T - (0.0065*200);
0023       % barometric formula
0024       P1 = Ppascals*(T/(T-(0.0065*(H-1222))))^((9.81*0.0289644)/(-8.3143*0.0065));
0025       P1 = P1*0.01;
0026           
0027       q = 1 - (T/647.096);
0028       % water vapour pressure formula (at saturation) gives wvp in
0029       % millibars
0030       wvp = 10000*22.064*exp((647.096/T)*((-7.85951783*q) + ...
0031             (1.84408259*(q^1.5)) + (-11.7866497*(q^3)) + ...
0032             (22.6807411*(q^3.5)) + (-15.9618719*(q^4)) + (1.80122502*(q^7.5))));
0033       wvpHg = wvp * 0.750063755; % convert from mbar to mmHg
0034       wvpHg = wvpHg*exp(-H*1000/6); % take from sea level to whatever height we're at
0035       
0036       % sea level water vapour density (saturation density * rel humidity
0037       % ), p120 Allen astrophysical quantities
0038      R = RelH*(2.886*10^(-4)*wvpHg/To) * 10^6; % factor 10^6 is to convert from g/cm3 to g/m3
0039                
0040       DeltaVh2o = 2.96*(P1/1013)*((300/T)^0.626)*(1 +(0.018*R*T/P1));     
0041      
0042       Kh2o = R*V*V*DeltaVh2o*(T^(-3/2))*((7.18*exp(-644/T)/T)*(1/((494.40190-V*V)^2 ...
0043              + 4*V*V*DeltaVh2o*DeltaVh2o)) + (2.77*10^(-8))); 
0044                   
0045      if P1 >= 333
0046          Gama = 0.59;
0047       elseif P1 < 333 
0048              Gama = 0.59*(1+0.0031*(333-P1));
0049       end
0050       
0051       Gama2 = Gama*(P1/1013)*((300/T)^0.85);
0052       
0053       con = (10^(6))*(log10(exp(1)));
0054    
0055     factor = (1/((V-60)^2 + Gama2^2)) + (1/(V^2 + Gama2^2)); 
0056     Ko2 = 0.011*V^2*(P1/1013)*((300/T)^2)*Gama2*factor/con;
0057             
0058       Both = Kh2o + Ko2;
0059 
0060       BothDiff = Both*20000;
0061 
0062 elseif  H > 11000
0063     
0064         T = T;
0065        
0066         P2 = P1*100*exp((-0.0289644*9.81*200)/(8.3144*T));
0067         P2 = P2*0.01;
0068     
0069         q = 1 - (T/647.096);
0070 
0071         wvp = 0.01*22064000*exp((647.096/T)*((-7.85951783*q) + ...
0072             (1.84408259*(q^1.5)) + (-11.7866497*(q^3)) + ...
0073             (22.6807411*(q^3.5)) + (-15.9618719*(q^4)) + (1.80122502*(q^7.5))));
0074         wvpHg = wvp* 0.750063755;
0075         wvpHg = wvpHg*exp(-H*1000/6); 
0076       
0077         R = RelH*(2.886*10^(-4)*wvpHg/T) * 10^6;
0078               
0079         DeltaVh2o = 2.96*(P2/1013)*((300/T)^0.626)*(1 +(0.018*R*T)/P2);
0080         
0081         Kh2o = R*V*V*DeltaVh2o*(T^(-3/2))*((7.18*exp(-644/T))/T* ... 
0082         (1/((494.40190-V*V)^2 + 4*V*V*DeltaVh2o*DeltaVh2o)) + (2.77*10^(-8)));
0083     
0084       if P2 >= 333
0085          Gama = 0.59;
0086       elseif P2 < 333 
0087              Gama = 0.59*(1+0.0031*(333-P2));
0088       end
0089       
0090       Gama2 = Gama*(P2/1013)*((300/T)^0.85);
0091       
0092       con = (10^(6))*(log10(exp(1)));
0093       
0094     factor = (1/((V-60)^2 + Gama2^2)) + (1/(V^2 + Gama2^2)); 
0095     Ko2 = 0.011*V^2*(P2/1013)*((300/T)^2)*Gama2*factor/con;          
0096    
0097     Both = Kh2o + Ko2;
0098         
0099         BothDiff = Both*20000; 
0100 end 
0101  
0102 L(1,((j-1222)/200)+1) = BothDiff;
0103 
0104 end
0105 % Sum everything but the very first element
0106 
0107 TheoOp = sum(L) - L(1);
0108 
0109 
0110 end