disp('gasspec.m')
%  Reads GASEX98 data set.....generates spectra of H2O and CO2 from Licor and OakRidge

if hhh>9,
   if mmm>9,
      zz=[num2str(ddd),num2str(hhh),num2str(mmm),'.dat'];
   else
      zz=[num2str(ddd),num2str(hhh),'0',num2str(mmm),'.dat'];
   end;
else
   if mmm>9,
      zz=[num2str(ddd),'0',num2str(hhh),num2str(mmm),'.dat'];
   else
      zz=[num2str(ddd),'0',num2str(hhh),'0',num2str(mmm),'.dat'];
   end;
end;

if odddrive==1,
   filename=[drive,'\gasex98\',zz];fid=fopen(filename,'r');
else
   filename=[drive,'\gasex\',zz];fid=fopen(filename,'r');
end
filename
clear zz

[c_to_v,dt,fc,rotor,deviation,gains,offsets,R,G] = const(20);
fc = 1./60;

%************************************************************
%   Read in the data from the CD or HD
%************************************************************

[yd,Sonic,Met,GPS,Campbell,YSI,Licor,Imet,Waves,Licorsea,Acoeff]=readgas(fid,ddd,hhh,mmm);
Waves=Waves(1:1200);
clear Acoeff

LicorseaCO2=Licorsea(2,:);

[Corrected,platform,deg_rate,Ubulk,Vbulk,Ustream] = voltophys(gains,offsets,c_to_v,Sonic,deviation,rotor);
clear rotor

Tsonic=Sonic(5,:);Tsonic=Tsonic./50;
u=Sonic(2,:);u=u./100;
v=Sonic(3,:);v=v./100;
w=Sonic(4,:);w=w./100;
u=-u;

Uu=u.*cos(deviation)-v.*sin(deviation);   %+ from stern to bow (from South to North)
Vv=v.*cos(deviation)+u.*sin(deviation);   %+ from starboard to port (from East to West)
Ustream=mean(sqrt(Uu.^2 + Vv.^2));
Ubulk=mean(Uu);   %mean Ubow
Vbulk=mean(Vv);   %mean Vbow

u=Uu;
v=Vv;
clear Uu Vv

Tsonic=(Tsonic.*Tsonic+(u.*u+w.*w)./2+v.*v)./403-273.15;

Tsonic=resample(Tsonic,12000,12500);

u=resample(u,12000,12500);
v=resample(v,12000,12500);
w=resample(w,12000,12500);

platform1=resample(platform(1,:),12000,12500);
platform2=resample(platform(2,:),12000,12500);
platform3=resample(platform(3,:),12000,12500);
clear platform

deg_rate1=resample(deg_rate(1,:),12000,12500);
deg_rate2=resample(deg_rate(2,:),12000,12500);
deg_rate3=resample(deg_rate(3,:),12000,12500);
clear deg_rate

if length(u)>12000,
   u=u(1:12000);
end
if length(v)>12000,
   v=v(1:12000);
end
if length(w)>12000,
   w=w(1:12000);
end
if length(Tsonic)>12000,
   Tsonic=Tsonic(1:12000);
end
if length(deg_rate1)>12000,
   deg_rate1=deg_rate1(1:12000);
end
if length(deg_rate2)>12000,
   deg_rate2=deg_rate2(1:12000);
end
if length(deg_rate3)>12000,
   deg_rate3=deg_rate3(1:12000);
end
if length(platform1)>12000,
   platform1=platform1(1:12000);
end
if length(platform2)>12000,
   platform2=platform2(1:12000);
end
if length(platform3)>12000,
   platform3=platform3(1:12000);
end

platform=[platform1 
   platform2
   platform3];
clear platform1 platform2 platform3

deg_rate=[deg_rate1 
   deg_rate2
   deg_rate3];
clear deg_rate1 deg_rate2 deg_rate3

%************************************************************
%  Get the mean T, Tsea, and RH from the WHOI system.
%************************************************************

[tair,rh]=meantrh(gains,offsets,c_to_v,Sonic);

rh=resample(rh,12000,12500);
tair=resample(tair,12000,12500);

clear Sonic scomp ccomp

%************************************************************
%  Get the ophir and the  Licor
%************************************************************
[LicorCO2,LicorH2O,pair]=co2aircalc(Licor,1000);

Ophir=Met(9:10,:);
[qophir,rhoair,qflag,sig_clear] = getophir(Ophir,tair,rh,mean(pair));

refclear=Met(9,:);
refabsorb=Met(10,:);

%************************************************************
%  Get the Oak Ridge
%************************************************************
   
oakCO2 = Met(3,:)*5/2047;
oakH2O = Met(4,:)*5/2047;

oakCO2 = oakCO2.*96.9+441;            %mg/m3
oakH2O = oakH2O.*oakH2O.*1.145+oakH2O.*1.888+0.43;  %g/m3

oakCO2ref = Met(5,:)*5/2047;
oakH2Oref = Met(6,:)*5/2047;

if length(oakCO2)>12000,
   oakCO2=oakCO2(1:12000);
end
if length(oakH2O)>12000,
   oakH2O=oakH2O(1:12000);
end
if length(LicorCO2)>12000,
   LicorCO2=LicorCO2(1:12000);
end
if length(LicorH2O)>12000,
   LicorH2O=LicorH2O(1:12000);
end

oakH2O=oakH2O./1000;      %kg/m3
oakCO2=oakCO2./1000000;   %kg/m3

LicorCO2=LicorCO2./1000000.*1.964.*273.15./293.15;   %kg/m3
LicorH2O=LicorH2O./1000;    %kg/m3

rc=oakCO2-mean(oakCO2)+mean(oakCO2).*(1.61.*(oakH2O-mean(oakH2O))./1.25+(1+1.61.*mean(oakH2O)./1.25).*(Tsonic-mean(Tsonic))./293.15);

%************************************************************
%  Get the ship speed from the GPS  
%************************************************************
   
yday=fix(yd(1));
clear yd
[yd2,gpsshipn,gpsshipe,cogm,sogm,latm,longm]=gps_vels(GPS,yday);
     
%clear GPS cogm sogm
   
%************************************************************
%  Get the Campbell info   
%************************************************************
   
[Tcamp,Tair,RH,Tsea,IR1,IR2,Solar1,Solar2] = GetCbell(Campbell);
if isempty(Solar1),
   Solar1=0;
end;
if isempty(Solar2),
   Solar2=0;
end;
if isempty(IR1),
   IR1=0;
end;
if isempty(IR2),
   IR2=0;
end;
if mean(IR2)>1000,
   IR2=IR1;
end
if isempty(RH),
   RH=0;
end;
if isempty(Tair),
   Tair=0;
end;
if isempty(Tcamp),
   Tcamp=0;
end;
if isempty(Tsea),
   Tsea=0;
end;

clear Campbell

badgps=0;
%************************************************************
%  Get the IMET info   
%************************************************************
[ringyro,pcodespeed,rwind,rdir,pair,tair3,rh3,Solar3,salt3,rainrate,precip,tsea3, ...
      pcodelong,pcodelat,pcodeshipn,pcodeshipe,badgps] = GetIMET(Imet,yd2,ddd,hhh,mmm);

for ii=2:length(tair3);
   if tair3(ii)>25,
      tair3(ii)=tair3(ii-1);
   end;
   if tair3(ii)<10,
      tair3(ii)=tair3(ii-1);
   end;
end

%clear Imet 
   
%************************************************************
%  Get the YSI info   
%************************************************************
   
[tsea2,salt2,DO,PH] = GetYSI(YSI);
for ii=1:length(tsea2)
   if abs(tsea2(ii))>1000,
      tsea2(ii)=15;
   end
end

clear YSI   
   
%************************************************************
%  Speed up the heading, speed, and scalar data
%************************************************************
   
[yd2,gpsshipn,gpsshipe,cogm,sogm,latm,longm]=gps_vels(GPS,yday);

  
if badgps==0,
   ringyro = unwrap(ringyro);
   [gyro] = faster2(ringyro,20,1,12000);
   [lat] = faster2(pcodelat,20,1,12000);
   [long] = faster2(pcodelong,20,1,12000);
   [shipn] = faster2(pcodeshipn,20,1,12000);
   [shipe] = faster2(pcodeshipe,20,1,12000);
   %[shipspd] = faster2(pcodespeed,20,1,12000);
   clear pcodespeed
   clear ringyro
   cog = atan3(shipe,shipn)*180/pi;
   sog = sqrt(shipe.^2 + shipn.^2);
else
   gyro=(1:12000).*0;
   [lat] = faster2(latm,20,1,12000);
   [long] = faster2(longm,20,1,12000);
   [shipn] = faster2(gpsshipn,20,1,12000);
   [shipe] = faster2(gpsshipe,20,1,12000);
   for ii=1:length(shipn),
      if abs(shipn(ii))>15,
         shipn(ii)=0;
      end
   end
   for ii=1:length(shipe),
      if abs(shipe(ii))>15,
         shipe(ii)=0;
      end
   end
   cog = atan3(shipe,shipn)*180/pi;
   sog = sqrt(shipe.^2 + shipn.^2);
end;


bf = ones(1,10*20)/(10*20);          %10 second filter
af = 1;

north = filtfilt(bf,af,shipn);
east = filtfilt(bf,af,shipe);

%************************************************************
%  Compute the Euler Angles using UCI code
%************************************************************
   
clear shipspd speed

[bhi,ahi,bf2,af2,tau2,spinup2]=filtcoef(20,fc);
clear spinup2 af2 bf2

[euler,dr]=angles(ahi,bhi,20,platform,deg_rate,gyro);    %!!!!!!!!!!

leuler=length(euler);
bike = std(euler(2,100:leuler-100))/std(euler(1,100:leuler-100));

for ii=1:length(gyro);
   if gyro(ii)<0,
      gyro(ii)=gyro(ii)+2.*pi;
   end
end

nrel = sog.*cos(gyro);
erel = sog.*sin(gyro);
  
clear deg_rate speed 
     
%************************************************************
%  Compute the platform velocity and position
%************************************************************
   
   [uvwplat,xyzplat]=accels(bhi,ahi,20,platform,euler);      %!!!!!!!!!!
   
clear ahi bhi  
%  I AM USING THE GPS VELOCITY FOR THE TIME BEING

uvwplat(1,:) = uvwplat(1,:); %+ north;
uvwplat(2,:) = uvwplat(2,:); %+ east;
heave = xyzplat(3,:);
clear platform xyzplat
%************************************************************
%  Compute the platform velocity and position
%************************************************************

clear Corrected
sonics=[u
      v
      w];
%clear u v w
[uvw] = sonic(sonics,dr,euler,uvwplat,R);       %!!!!!!!!!!
clear R

uvw(1,:)=uvw(1,:)+uvwplat(1,:);
uvw(2,:)=uvw(2,:)+uvwplat(2,:);
uvw(3,:)=uvw(3,:)+uvwplat(3,:);
%clear uvwplat
  
clear euler dr sonics 
  
reldir1 = atan3(Vbulk,Ubulk)*180/pi;
relwind = sqrt(Vbulk*Vbulk+Ubulk*Ubulk);
unorth = mean(uvw(1,:));						%ship speed not removed
uwest = mean(uvw(2,:));
wdir = atan3(uwest, -unorth)*180/pi;
wspd = sqrt(uwest*uwest + unorth*unorth);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

mpair=mean(pair);
if abs(mpair)>2000,
   mpair=1005;
end

Tetl=mean(Tair);				%mean air temperature from Campbell
Twhoi=mean(tair);				%mean air temperature from WHOI
Timet=mean(tair3);			%mean air temperature from IMET

RHetl=mean(RH);					%mean RH from Campbell
RHwhoi=mean(rh)*100;				%mean RH from WHOI
RHimet=mean(rh3);					%mean RH from IMET

Qbar=0.622*RHwhoi/100*6.112*exp((17.67*Twhoi)/(Twhoi+243.5))/1013.25*1000;

rho=mpair*100/(287*(Twhoi+273.15)*(1+0.61*Qbar));
if rho>2 | rho<1,
   rho=1.25;
end;

Tseaship=mean(tsea3);
Tseaysi=mean(tsea2);
Tseasnake=mean(Tsea);

Esat=6.112.*exp(17.67.*Tseaysi./(Tseaysi+243.5))*0.98;
Qsea=Esat*0.622/(mpair-0.378*Esat)*1000;

mean_oakCO2=mean(oakCO2./rho);
med_oakCO2=median(oakCO2./rho);
std_oakCO2=std(oakCO2./rho);
min_oakCO2=min(oakCO2./rho);
max_oakCO2=max(oakCO2./rho);

mean_rc=mean(rc./rho);
med_rc=median(rc./rho);
std_rc=std(rc./rho);
min_rc=min(rc./rho);
max_rc=max(rc./rho);

mean_oakH2O=mean(oakH2O./rho);
med_oakH2O=median(oakH2O./rho);
std_oakH2O=std(oakH2O./rho);
min_oakH2O=min(oakH2O./rho);
max_oakH2O=max(oakH2O./rho);

mean_LicorCO2=mean(LicorCO2./rho);
med_LicorCO2=median(LicorCO2./rho);
std_LicorCO2=std(LicorCO2./rho);
min_LicorCO2=min(LicorCO2./rho);
max_LicorCO2=max(LicorCO2./rho);

mean_LicorH2O=mean(LicorH2O./rho);       %  kg/kg
med_LicorH2O=median(LicorH2O./rho);      %  kg/kg
std_LicorH2O=std(LicorH2O./rho);       %  kg/kg
min_LicorH2O=min(LicorH2O./rho);       %  kg/kg
max_LicorH2O=max(LicorH2O./rho);       %  kg/kg

mean_LicorseaCO2=mean(LicorseaCO2./rho);
med_LicorseaCO2=median(LicorseaCO2./rho);
std_LicorseaCO2=std(LicorseaCO2./rho);
min_LicorseaCO2=min(LicorseaCO2./rho);
max_LicorseaCO2=max(LicorseaCO2./rho);

mean_OphirH2O=mean(qophir./1000./rho);
med_OphirH2O=median(qophir./1000./rho);
std_OphirH2O=std(qophir./1000./rho);
min_OphirH2O=min(qophir./1000./rho);
max_OphirH2O=max(qophir./1000./rho);

mean_absorb=mean(refabsorb);
med_absorb=median(refabsorb);
std_absorb=std(refabsorb);
min_absorb=min(refabsorb);
max_absorb=max(refabsorb);

mean_clear=mean(refclear);
med_clear=median(refclear);
std_clear=std(refclear);
min_clear=min(refclear);
max_clear=max(refclear);

mean_Tsonic=mean(Tsonic);
med_Tsonic=median(Tsonic);
std_Tsonic=std(Tsonic);
min_Tsonic=min(Tsonic);
max_Tsonic=max(Tsonic);

north=mean(north);         %ship speed to north
east=mean(east);           %ship speed to north

mean_sog=mean(sog);        %mean ship speed
if mean_sog<1,
   mean_sog=0;
end
std_sog=std(sog);				%std ship speed
max_sog=max(sog);				%max ship speed
min_sog=min(sog);				%min ship speed

mean_cog=atan2(mean(sin(cog./180.*pi)),mean(cos(cog./180.*pi))).*180./pi;  	%mean course in degrees
std_cog=std(cog);				%std course
min_cog=min(cog);				%min course
max_cog=max(cog);				%max course
if mean_cog<0,
   mean_cog=mean_cog+360;
end
if max_cog<0,
   max_cog=max_cog+360;
end
if min_cog<0,
   min_cog=min_cog+360;
end

Rs=(median(Solar1)+median(Solar2))/2;
std_Rs=std((Solar1 + Solar2)/2);
min_Rs=(min(Solar1)+min(Solar2))/2;
max_Rs=(max(Solar1)+max(Solar2))/2;

Rs1=median(Solar1);
std_Rs1=std(Solar1);
min_Rs1=min(Solar1);
max_Rs1=max(Solar1);

Rs2=median(Solar2);
std_Rs2=std(Solar2);
min_Rs2=min(Solar2);
max_Rs2=max(Solar2);

med_Rss=median(Solar3);
std_Rss=std(Solar3);
min_Rss=min(Solar3);
max_Rss=max(Solar3);

Rl=(median(IR1)+median(IR2))/2;
std_Rl=std((IR1 + IR2)/2);
min_Rl=min((IR1 + IR2)/2);
max_Rl=max((IR1 + IR2)/2);

Rl1=median(IR1);
std_Rl1=std(IR1);
min_Rl1=min(IR1);
max_Rl1=max(IR1);

Rl2=median(IR2);
std_Rl2=std(IR2);
min_Rl2=min(IR2);
max_Rl2=max(IR2);

mean_wave=mean(Waves);        %
med_wave=median(Waves);       %
std_wave=std(Waves);				%
max_wave=max(Waves);				%
min_wave=min(Waves);				%

uvw(1,:)=uvw(1,:)+mean_sog.*cos(mean_cog/180*pi)+mean(uvwplat(1,:));    %total correction
uvw(2,:)=uvw(2,:)-mean_sog.*sin(mean_cog/180*pi)+mean(uvwplat(2,:));
uvw(3,:)=uvw(3,:)+mean(uvwplat(3,:));

[U,V,W,alpha,beta]=rotcorr(u,v,w);			%uncorrected streamwise rotation
[Uu,Vv,Ww,a,b]=rotcorr(uvw(1,:),uvw(2,:),uvw(3,:));			%corrected streamwise rotation

Ubarstr=mean(U);
Ustr=sqrt(mean(u).^2 + mean(v).^2 + mean(w).^2);
Ubarcorr=sqrt(mean(uvw(1,:)).^2 + mean(uvw(2,:)).^2 + mean(uvw(3,:)).^2);

reldir = atan3(mean(v),mean(u))*180/pi;
relwind = sqrt(mean(u).^2 + mean(v).^2 +mean(w).^2);
if reldir<0,
   reldir=reldir+360;
end;

unorth = mean(uvw(1,:));                     %corrected
uwest = mean(uvw(2,:));                      %corrected

wdir = atan3(uwest, -unorth)*180/pi;		   %true wind direction
if wdir<0,
   wdir=wdir+360;
end

mean_head=atan2(mean(sin(gyro./180.*pi)),mean(cos(gyro./180.*pi)))*180/pi;	%mean heading
std_head=std(gyro./180.*pi)*180/pi;		%std heading
min_head=min(gyro./180.*pi)*180/pi;		%min heading	
max_head=max(gyro./180.*pi)*180/pi;		%max heading
if mean_head<0,
   mean_head=mean_head+360;
end
if max_head<0,
   max_head=max_head+360;
end
if min_head<0,
   min_head=min_head+360;
end

wu=cov(Uu,Ww);					      %u'w'
wv=cov(Vv,Ww);					      %v'w'
wt=cov(Ww,Tsonic);		         %w't'
wq=cov(Ww,qophir./1000./rho);		%w'q'
wlh=cov(Ww,LicorH2O./rho);       % m/s * g/kg
wlc=cov(Ww,LicorCO2./rho);       % m/s * g/kg
woh=cov(Ww,oakH2O./rho);
woc=cov(Ww,rc./rho);

Sp1=std(uvwplat(1,:));
Sp2=std(uvwplat(2,:));
Sp3=std(uvwplat(3,:));

rain=mean(rainrate);
if abs(rain)<0.1,
   rain=0;
end;

Lat=mean(lat);
Lon=mean(long);

Us=0;zi=600;zu=19;zq=18;zt=18;
if isnan(Ubarcorr);
   Ubarcorr=0.5;
end

xship=[Ubarcorr Us Tseaysi Twhoi Qsea Qbar Rs Rl rain zi mpair zu zt zq];
yship=coul_lkb(xship);