SkyDipMethod2v2

 calculate Zenith Opacity from SkyDips using new filtered data

0001 % calculate Zenith Opacity from SkyDips using new filtered data
0002 
0003 % first apply alpha correction
0004 
0005 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0006 d = determineIndices(d);
0007 d = apparentAzEl(d);
0008 
0009 
0010 if(~isfield(d.antenna0.servo, 'equa')) 
0011   display('Calculating RA/DEC');
0012   long=-118.2822;
0013   lat=37.2339;
0014   
0015   az = d.antenna0.servo.apparent(:,1);
0016   el = d.antenna0.servo.apparent(:,2);
0017   jd=mjd2jd(d.antenna0.receiver.utc);
0018   [equa] = horiz_coo([pi/180*(az) pi/180*(el)],jd,[pi/180*(long) ...
0019         pi/180*(lat)],'e');
0020   d.antenna0.servo.equa=equa;
0021   clear equa;
0022   clear az;
0023   clear el;
0024 end
0025 
0026 s = 'FILTERED';
0027 dO = assembleAlphaStreams(d,s); 
0028 [tA,Ac,Gc,Af,Gf] = calculateAlpha(dO,s); 
0029 dO = applyAlpha(dO,tA,Ac,Gc,Af,Gf,s); 
0030 d = calculateStokes(dO,s);
0031 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0032 
0033 % now the sky dips
0034 
0035 d_dip = framecut(d, bitsearch(d.array.frame.features, 5, 'any'));
0036 
0037 Tatmos = (mean(d.array.weather.airTemperature)) + 273.15;
0038 Tamb = Tatmos - (0.0065*2000); % ave temp in sky (2 km up)
0039 
0040 Tcmb = 2.725;
0041 
0042 y = d_dip.antenna0.receiver.dataF(:,1) * 2.7; % noise diode is 2.7 K
0043 
0044 el = d_dip.antenna0.servo.apparent(:,2);
0045 
0046 q = length(el);
0047 
0048 time = d_dip.antenna0.receiver.utc - d_dip.antenna0.receiver.utc(1);
0049 time = time.*24;
0050 
0051 yy6 = spline(time, el,time);
0052 
0053 [maxtab] = peakdet(yy6,0.1,time);
0054 
0055 Divider = length(maxtab(:,1));
0056     
0057 
0058 if Divider ~= 0;
0059    NewY = (length(el))/Divider;
0060    m = round(NewY);
0061    
0062    y = y(1:m);
0063    x = 1./(sind(el(1:m)));
0064 
0065   
0066 elseif Divider == 0; 
0067        y = y;
0068        x = 1./(sind(el));
0069 end 
0070  
0071   GoodX = (x <= 2.0);
0072   
0073   n = length(y);
0074   
0075   [gradient2(n,:) intercept2(n,:)] = linfit(x(GoodX),y(GoodX));
0076   
0077     finY{n} = y;
0078     finX{n} = x;
0079   
0080     finYR1{n} = (gradient2(n,1)*x) + intercept2(n,1);
0081 
0082  % plot(x,y,'.',x,finYR1{n})
0083  % legend('Intensity Channel One')
0084  % xlabel('cosec(elevation)')
0085  % ylabel('Temperature (K)')
0086   
0087   
0088 kk = [x(GoodX),y(GoodX)];
0089 
0090 hh = kk(:,1);
0091 jj = kk(:,2);
0092 
0093 
0094  % display ('Channel One gradient and y intercept:')
0095  
0096   p = polyfit(hh,jj,1);
0097   
0098   gradientOne = p(:,1);
0099      
0100 OpacityOne = gradientOne/(Tamb - Tcmb);
0101 
0102 %getting the errors from the chi goodness of fit test
0103 
0104 len = length(y(GoodX));
0105 H1 = finYR1{n};
0106 
0107 chi_sq1 = sum(((y(GoodX,1) - H1(1:len)).^2)./H1(1:len));
0108 
0109 return
0110 
0111 
0112 
0113   
0114