clear all
told = 0.048;
fold = 1/told;                              % Sample Rate 
fs = 2;
fc = 1/60;

[c_to_v,dt,fc,rotor,deviation,gains,offsets,R,G] = const(fs);

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

cddir='e:\';
dirname='gasex\';
directory = [cddir dirname];
cdx=1;
[filenames,nfiles] = readdir(directory,cdx);
[filenames] = lexsort(filenames);
disp(['Total files:  ' num2str(nfiles)]);
plotit=0;
   
%************************************************************
%   Read in the data from the CD or HD
%************************************************************

yd =[];
mdngt =[];
Sonics =[];
Ophir = [];
Oakridge = [];
GPS = [];
Campbell = [];
YSI = [];
Licor = [];
Imet = [];
Waves = [];
Licorsea = [];

ydy =[];
mdngty =[];
Sonicsy =[];
Ophiry = [];
Oakridgey = [];
GPSy = [];
Campbelly = [];
YSIy = [];
Licory = [];
Imety = [];
Wavesy = [];
Licorseay = [];

averages=3;
averages1=averages-1;
nstart=1;

%*********************************************************
%  Get the first long time series
%*********************************************************
for nfile=nstart:averages1+nstart-1,
   filename=[setstr(filenames(nfile,:)) '.dat'];
   if str2num(filename(4:7)) < 11,
      midnight=1;
   else,
      midnight=0;
   end;
   
   %************************************************************
   %   Read in the data from the CD or HD
   %************************************************************
   
   [ydx,Sonicsx,Metx,GPSx,Campbellx,YSIx,Licorx,Imetx,...
         Wavesx,Licorseax,Acoeffx]=readco2(cddir,dirname,filename);
   mdngty = [mdngty midnight]; clear midnight 
   ydy = [ydy ydx]; clear ydx
   Sonicsy =[Sonicsy Sonicsx(:,1:12500)]; clear Sonicsx
   Oakridgey = [Oakridgey Metx(3:7,:)]; 
   Ophiry = [Ophiry Metx(9:10,:)]; clear Metx
   GPSy = [GPSy GPSx]; clear GPSx
   Campbelly = [Campbelly Campbellx]; clear Campbellx
   YSIy = [YSIy YSIx(:,1:60)]; clear YSIx
   Licory = [Licory Licorx(:,1:12000)]; clear Licorx
   Imety = [Imety Imetx(:,1:300)]; clear Imetx
   Wavesy = [Wavesy Wavesx(:,1:1200)]; clear Wavesx
   Licorseay = [Licorseay Licorseax]; clear Licorseax
end

%Numer of points in time 

LL1 = 12500;    %Sonics    21 Hz
LL2 = 12000;    %Oak Ridge 20 Hz (and Ophir)
LL3 = 600;      %GPS        1 Hz
LL4 = 300;      %Campbell  .5 Hz
LL5 = 60;       %YSI       .1 Hz
LL6 = 12000;    %Licor     20 Hz
LL7 = 300;      %Imet      .5 Hz
LL8 = 1200;     %Waves      2 Hz
LL9 = 300;      %Licorsea  .5 Hz

MM1 = length(Sonicsy) + LL1;
MM2 = length(Oakridgey) + LL2;
MM3 = length(GPSy) + LL3;
MM4 = length(Campbelly) + LL4;
MM5 = length(YSIy) + LL5;
MM6 = length(Licory) + LL6;
MM7 = length(Imety) + LL7;
MM8 = length(Wavesy) + LL8;
MM9 = length(Licorseay) + LL9;

for nfile=nstart+averages1:nfiles
   filename=[setstr(filenames(nfile,:)) '.dat'];
   if str2num(filename(4:7))<11
      midnight=1;
   else
      midnight=0;
   end
   
   %************************************************************
   %   Read in the data from the CD or HD
   %************************************************************
   
   [ydx,Sonicsx,Metx,GPSx,Campbellx,YSIx,Licorx,Imetx,...
         Wavesx,Licorseax,Acoeffx]=readco2(cddir,dirname,filename);
   mdngt = [mdngty midnight]; clear midnight mdngty 
   yd = [ydy ydx]; clear ydx ydy
   Sonics =[Sonicsy Sonicsx(:,1:12500)]; clear Sonicsx Sonicsy
   Oakridge = [Oakridgey Metx(3:7,:)]; clear Oakridgey
   Ophir = [Ophiry Metx(9:10,:)]; clear Metx Ophiry
   GPS = [GPSy GPSx]; clear GPSx GPSy
   Campbell = [Campbelly Campbellx]; clear Campbellx Campbelly
   YSI = [YSIy YSIx(:,1:60)]; clear YSIx YSIy
   Licor = [Licory Licorx(:,1:12000)]; clear Licorx Licory
   Imet = [Imety Imetx(:,1:300)]; clear Imetx Imety
   Waves = [Wavesy Wavesx(:,1:1200)]; clear Wavesx Wavesy
   Licorsea = [Licorseay Licorseax]; clear Licorseax Licorseay
      
   %*********************************************************
   %  Save Data for next record
   %*********************************************************
   
   if (yd(averages)-yd(1))*86400<599*averages1
      
      ydy(1:averages1) = yd(2:averages); clear yd
      mdngty(1:averages1) = mdngt(2:averages); clear mdngt
      Sonicsy(:,1:averages1*LL1) = Sonics(:,LL1+1:MM1); clear Sonics
      Ophiry(:,1:averages1*LL2) = Ophir(:,LL2+1:MM2); clear Ophir
      Oakridgey(:,1:averages1*LL2) = Oakridge(:,LL2+1:MM2); clear Oakridge
      GPSy(:,1:averages1*LL3) = GPS(:,LL3+1:MM3); clear GPS
      Campbelly(:,1:averages1*LL4) = Campbell(:,LL4+1:MM4); clear Campbell 
      YSIy(:,1:averages1*LL5) = YSI(:,LL5+1:MM5); clear YSI
      Licory(:,1:averages1*LL6) = Licor(:,LL6+1:MM6); clear Licor
      Imety(:,1:averages1*LL7) = Imet(:,LL7+1:MM7); clear Imet
      Wavesy(:,1:averages1*LL8) = Waves(:,LL8+1:MM8); clear Waves
      Licorseay(:,1:averages1*LL9) = Licorsea(:,LL9+1:MM9); clear Licorsea
   else
      ydy(1:averages1) = yd(2:averages); 
      mdngty(1:averages1) = mdngt(2:averages);
      Wavesy(:,1:averages1*LL8) = Waves(:,LL8+1:MM8);
      Sonicsy(:,1:averages1*LL1) = Sonics(:,LL1+1:MM1);
      
   %clear Waves
   
   %************************************************************
   %  Get the mean air, sea and humidity data from WHOI system.
   %  Channels 12 and 16 are free afterwards
   %************************************************************
   
   [tair,rh] = mean1(gains,offsets,c_to_v,Sonics);
   
   %************************************************************
   %   Average the data down to fs Hz.  Compass channels need
   %   special attention.
   %************************************************************
   
   comp = (Sonics(6,:) - offsets(6))*gains(6);   
   Sonics(6,:) = cos(comp);
   Sonics(12,:) = sin(comp);
   clear comp
   
   Sonicsx = resampler(Sonics(1:15,:),fs,fold);
   clear Sonics
   Sonics = Sonicsx;
   clear Sonicsx 
   
   comp = atan3(Sonics(12,:),Sonics(6,:));
   Sonics(6,:) = comp./gains(6) + offsets(6);
   
   [row,col] = size(Sonics);
   tot = col;
   
   %************************************************************
   %   Now the Ophir 
   %************************************************************
   
   Ophirx = resampler(Ophir,fs,20);
   Ophiry(:,1:averages1*LL2) = Ophir(:,LL2+1:MM2);
   clear Ophir
   Ophir = Ophirx;
   clear Ophirx
   
   %************************************************************
   %   Now the Oak Ridge 
   %************************************************************
   
   Oakridgex = resampler(Oakridge,fs,20);
   Oakridgey(:,1:averages1*LL2) = Oakridge(:,LL2+1:MM2);
   clear Oakridge
   Oakridge = Oakridgex;
   clear Oakridgex
   
   %************************************************************
   %   Now the Licor
   %************************************************************
   
   Licorx = resampler(Licor,fs,20);
   Licory(:,1:averages1*LL6) = Licor(:,LL6+1:MM6);
   clear Licor
   Licor = Licorx;
   clear Licorx
   
   %************************************************************
   %  Get the ship speed from the GPS  
   %************************************************************
   
   yday=fix(yd(1));
   [yd2,gpsshipn,gpsshipe,cogm,sogm,latm,longm]=gps_vels(GPS,yday,mdngt);
   GPSy(:,1:averages1*LL3) = GPS(:,LL3+1:MM3);
      
   clear GPS cogm sogm
   
   %************************************************************
   %  Get the Campbell info   
   %************************************************************
   
   Campbelly(:,1:averages1*LL4) = Campbell(:,LL4+1:MM4);   
   [Tcamp,Tair,RH,Tsea,IR1,IR2,Solar1,Solar2] = GetCbell(Campbell);
   figure;plot(Tsea)
%   disp('pause');pause
   
   clear Campbell
   
   %************************************************************
   %  Get the IMET info   
   %************************************************************
   imeton=1;
   [ringyro,pcodespeed,rwind,rdir,pair,tair3,Solar3,salt3,rainrate,precip,tsea3, ...
      pcodelong,pcodelat,pcodeshipn,pcodeshipe] = GetIMET(Imet,yd2);
   Imety(:,1:averages1*LL7) = Imet(:,LL7+1:MM7);
   
   clear Imet 
   
   %************************************************************
   %  Get the YSI info   
   %************************************************************
   
   [tsea2,salt2,DO,PH] = GetYSI(YSI);
   YSIy(:,1:averages1*LL5) = YSI(:,LL5+1:MM5);
   
   clear YSI   
   
   %************************************************************
   %  Speed up the heading, speed, and scalar data
   %************************************************************
   
   yd1 = yd2(1) + (0:tot-1)*(dt/86400);
   tot1=tot-1;
   
   if imeton==1
      [lat] = faster2(pcodelat,fs,1,tot);
      [long] = faster2(pcodelong,fs,1,tot);
      [shipn] = faster2(pcodeshipn,fs,1,tot);
      [shipe] = faster2(pcodeshipe,fs,1,tot);
      [shipspd] = faster2(pcodespeed,fs,1,tot);
      ringyro = unwrap(ringyro);
      [gyro] = faster2(ringyro,fs,1,tot);
      cog = atan3(shipe,shipn)*180/pi;
      sog = sqrt(shipe.^2 + shipn.^2);
   else
      [lat] = faster2(latm,fs,1,tot);
      [long] = faster2(longm,fs,1,tot);
      [shipn] = faster2(gpsshipn,fs,1,tot);
      [shipe] = faster2(gpsshipe,fs,1,tot);
      [shipspd] = zeros(size(yd1));
      cog = atan3(shipe,shipn)*180/pi;
      sog = sqrt(shipe.^2 + shipn.^2);
   end
      
   clear latm longm gpsship gpsshipn gpsshipe yd2 pcodespeed
   clear pcodelat pcodelong pcodeshipn pcodeshipe ringyro
   
   %************************************************************
   %   Convert motion variable to physical units
   %************************************************************
   
   [Corrected,platform,deg_rate,Ubulk,Vbulk,Ustream] = ... 
      voltophys(gains,offsets,c_to_v,Sonics,deviation,rotor);
   
   clear Sonics
   
   %************************************************************
   %  Get the Ophir data
   %************************************************************
   
   [qophir,rhoair,qflag,sig_clear] = getophir(Ophir,tair,rh,pair);
   clear Ophir
   
   %************************************************************
   %  Get the Oak Ridge data
   %************************************************************
   
   Vco2 = Oakridge(1,:)*5/2048;
   Vh2o = Oakridge(2,:)*5/2048;

   Tcell = mean(Tair)+273.15;
   Pcell = pair/10;
   P0 = 101.325     % from Auble ref pressure for CO2? and H2O?
   T0 = 20+273.15   % from Auble ref temp for CO2? and H2O?

   Vh2o = Vh2o.*(P0./Pcell);                                   %Correct the volts for p
   a_h2ooak = 1.855.* Vh2o .* Vh2o + 1.606 .* Vh2o +0.96;      %g/m3  absolute hum
   q_h2ooak = a_h2ooak ./ 1.2;                                 %g/kg specific humidity
   w_h2ooak = 29 .* q_h2ooak ./ (18 + 11 .* q_h2ooak ./1000);  %mmol/mol  (PPT)

   M = 29 .* (1 - w_h2ooak ./ 1000) + 18 .* w_h2ooak ./ 1000;  % Molecular weight of moist air
   Eair = w_h2ooak .* Pcell./1000;

   H2O_corr = 1 + 0.5 .* w_h2ooak./1000;               %pressure broadening correction

   Vco2 = Vco2 .* (P0./Pcell./H2O_corr);               %correct for pressure broadening
   a_co2oak = 107 .* Vco2 + 391;                       %mg/m^3   "absolute"
   q_co2oak = a_co2oak ./ 1.2;                         %mg/m^3   "absolute"
   w_co2oak = 29 .* q_co2oak ./ 44;                    %umol/mol  (PPM)   44 is molecular weight of CO2,  29 is molecular weight of dry air

   w_co2oak = w_co2oak .* H2O_corr .* (Tcell./T0);     %umol/mol  (PPM)

   w_co2oak = w_co2oak ./ (1 - w_h2ooak./1000);        %correct for dilution (dry it)
   figure;plot(q_h2ooak);title(['H2O  = ',filename]);
   figure;plot(w_co2oak);title(['CO2  = ',filename]);

%   clear Oakridge
     
   %************************************************************
   %  Compute the Euler Angles using UCI code
   %************************************************************
   
   bf = ones(1,10*fs)/(10*fs);          %10 second filter
   af = 1;
   
   pnav = Corrected(5,:);
   cnav = filtfilt(bf,af,cos(pnav));
   snav = filtfilt(bf,af,sin(pnav));
   pnav = atan3(snav,cnav)*180/pi;
   gyro2 = pnav-12-46*sin(pnav*pi/180);
   gyro2 = gyro2*pi/180;
   clear pnav snav cnav
   
   i=find(gyro2<0);
   gyro2(i)=gyro2(i)+2*pi;
   i=find(gyro2>0);
   gyro2(i)=gyro2(i)-2*pi;   
   
   if imeton==0 
      gyro = gyro2;
   end
   
   speed = filtfilt(bf,af,shipspd);
   north = filtfilt(bf,af,shipn);
   east = filtfilt(bf,af,shipe);
   clear shipspd shipe shipn
   
   [euler,dr]=angles(ahi,bhi,fs,platform,deg_rate,gyro);    %!!!!!!!!!!

   leuler=length(euler);
   bike = std(euler(2,100:leuler-100))/std(euler(1,100:leuler-100));
   nrel = speed.*cos(gyro);
   erel = speed.*sin(gyro);
   
   clear deg_rate speed 
     
   %************************************************************
   %  Compute the platform velocity and position
   %************************************************************
   
   [uvwplat,xyzplat]=accels(bhi,ahi,fs,platform,euler);      %!!!!!!!!!!
   
   %  I AM USING THE GPS VELOCITY FOR THE TIME BEING
   
   uvwplat(1,:) = uvwplat(1,:) + north;
   uvwplat(2,:) = uvwplat(2,:) - east;
   heave = xyzplat(3,:);
   dp = -rhoair*G*heave/100;
   clear platform xyzplat
   
   %************************************************************
   %  Compute the platform velocity and position
   %************************************************************
   
   sonics = Corrected(1:3,:);
   reldir = atan3(sonics(2,:),sonics(1,:))*180/pi;
   tv = Corrected(4,:);
   
   clear Corrected
   
   [uvw] = sonic(sonics,dr,euler,uvwplat,R);       %!!!!!!!!!!
   
   clear euler dr sonics 
   
   %************************************************************
   %  Get CO2 and additional H20 data   
   %************************************************************
   
   [co2air,h2oair,pair2] = co2aircalc(Licor,dp);
   Lag=fix(1.3*fs);
   Llic = length(Licor);
   Lag2 = Llic-Lag+1;
   Licorlag = zeros(size(co2air));
   Licorlag(1:Lag2) = co2air(Lag:Llic);
   co2licor = Licorlag;
   Licorlag(1:Lag2) = h2oair(Lag:Llic);
   h2olicor = Licorlag;
   co2sea = Licorsea(2,:);
      
   Licorseay(:,1:averages1*LL9) = Licorsea(:,LL9+1:MM9);
 
   clear Licor Licorsea Licorlag pair2
   
%   if plotit==1
%      clf;figure(1);
%      mh2oair = mean(h2olicor(1:200));
%      mh2ooak = mean(h2ooak(1:200));
%      plot(yd1(1:200),h2olicor(1:200)-mh2oair,'r',yd1(1:200),h2ooak(1:200)-mh2ooak,'g');
%      pause(1)
%   end
   
   
   reldir1 = atan3(Vbulk,Ubulk)*180/pi;
   relwind = sqrt(Vbulk*Vbulk+Ubulk*Ubulk);
   unorth = mean(uvw(1,:));
   uwest = mean(uvw(2,:));
   wdir = atan3(uwest, -unorth)*180/pi;
   wspd = sqrt(uwest*uwest + unorth*unorth);
%   if plotit==1
%      clf;figure(1);
%      subplot(4,1,1);
%      plot(yd1,uvw(1,:),'r',yd1,uvwplat(1,:),'g',yd1,north,'b');
%      subplot(4,1,2);
%      plot(yd1,uvw(2,:),'r',yd1,uvwplat(2,:),'g',yd1,-east,'b');
%      subplot(4,1,3);
%      plot(yd1,uvw(3,:),'r',yd1,uvwplat(3,:),'g');
%      subplot(4,1,4);
%      plot(yd1,gyro*180/pi,'g',yd1,cog,'m');
%   end
   
   points=1024;
   aves=fix(length(uvw)/points);
   edge = fix((length(uvw)-points*aves)/2)*2;
   
   UVW(1,1:tot-edge) = uvw(1,1+edge/2:tot-edge/2);
   UVW(2,1:tot-edge) = uvw(2,1+edge/2:tot-edge/2);
   UVW(3,1:tot-edge) = uvw(3,1+edge/2:tot-edge/2);
   TsV(1:tot-edge) = tv(1+edge/2:tot-edge/2);
   CO2(1:tot-edge) = co2licor(1+edge/2:tot-edge/2);
     
   [uvwde] = align_j(UVW);
   uu = detrend(uvwde(1,:));
   vv = detrend(uvwde(2,:));
   ww = detrend(uvwde(3,:));
   cc = detrend(CO2);
   tt = detrend(TsV);
   
   clear UVW uvwde TsV
   
   overlap=0;
   [cxx1,fr] = csd2(uu,ww,points,fs,hamming(points),overlap,'none');
   cxx = real(cxx1);
   cxx(1)=0;
   j=find(fr < min(1,relwind/17.8));
   uw = sum(cxx(j))*fr(2);
   clear cxx cxx1
   
   [cxx1,fr] = csd2(vv,ww,points,fs,hamming(points),overlap,'none');
   cxx = real(cxx1);
   cxx(1)=0;
   vw = sum(cxx(j))*fr(2);
   clear cxx cxx1
   
   [cxx1,fr] = csd2(tt,ww,points,fs,hamming(points),overlap,'none');
   cxx = real(cxx1);
   cxx(1)=0;
   wt = sum(cxx(j))*fr(2);
   clear cxx cxx1
   
   [cxx1,fr] = csd2(cc,ww,points,fs,hamming(points),overlap,'none');
   cxx = real(cxx1);
   cxx(1)=0;
   wc = sum(cxx(j))*fr(2);
   [Smoothed,Frq]=specsmoo(cxx,fs);
   
 %  [cxx1,fr] = psd2(cc,512,fs,hamming(512),0,'none');
 %  [pxx1,fr] = psd2(dp,512,fs,hamming(512),0,'linear');
   plotit=0;
%   if plotit==1
%      clf;figure(1);
%      semilogx(Frq,Frq.*Smoothed,'g')
%   end
   plotit=0;
%   pause(1)
   
   clear cxx cxx1 pxx1 heave uu vv ww tt cc
   
   siggyro = std(gyro)*180/pi;
   
   stuff=sprintf('YD: %2.5f RD: %2.0f SOG: %1.1f WS: %1.1f WD: %2.0f uw: %0.4f vw: %0.4f wc: %0.4f wt: %0.4f Cair: %2.1f Csea: %2.1f Cmin: %2.0f Cmax: %2.0f PR: %1.1f SigG %1.1f', ...
      yd(1),reldir1,mean(sog)*1.946,wspd,wdir,uw,vw,wc,wt,median(CO2),median(co2sea),min(CO2),max(CO2),bike,siggyro);
   disp(stuff);
   
   clear stuff CO2
   
   filemat=setstr(filenames(nfile-averages1,:));
   
%   eval(['save c:\gasex98\' filemat ' yd1 lat long gyro Ubulk Vbulk reldir sog cog north east nrel erel uvw tv co2sea co2air h2oair qophir Vco2 Vh2o pair tair Tair tair3 Tsea tsea2 tsea3 rh RH rhoair IR1 IR2 Solar1 Solar2 Solar3 salt2 salt3 rainrate precip DO PH sig_clear qflag bike a_co2oak q_co2oak w_co2oak a_h2ooak q_h2ooak w_h2ooak Oakridge']);
   %figure;plot(Ustream)
   
   clear yd yd1 lat long gyro gyro2 reldir cog sog uvw tv qophir 
   clear uvwplat sig_clear qflag north east dp Vh2o Vco2 bike
   clear Tcamp Tsea DO pair nrel erel co2licor h2olicor h2ooak
   clear Ubulk Vbulk Ustream Acoeff h2oair co2air co2sea rhoair
   clear pair tair Tair tair3 tsea2 tsea3 rh RH rhoair IR1 IR2 
   clear Solar1 Solar2 Solar3 salt2 salt3 rainrate precip DO PH
      
end      
end
   
