close('all');
fclose('all');
clear;

datasetv='_v1';
read_parameters_HiWinGS_2013;

expp='HiWinGS_2013';
drive='C:\data\CRUISES\HiWinGS_2013\Knorr\flux\Raw';
% drive_w='C:\Data\CRUISES\DYNAMO_2011\flux_data';  
drive_w='C:\data\CRUISES\HiWinGS_2013\flux_data';  

graphformat='.png';  %select graphics format files
graphdevice='-dpng'; %select graphic device


clearit=1;
printit=1;
year='2013';

yyyy='2013';

decorr=1;  %to decorrelate winds from motion...
dcname=[];
if decorr==1
    dcname='_vdecorr';
end;

rge=310:318; %283:301;

for ddd=rge
    [m,d]=yd2md(str2num(year), ddd);
    Vdate=[str2num(year) m d];
       
    for hhh=00:23;   %
      read_hp_HiWinGS_2013;      %
      read_shp_HiWinGS_2013;     %
      read_motion_HiWinGS_2013;  %
      read_son_HiWinGS_2013;     %
      read_gps_HiWinGS_2013;     %

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%      
      %uses the SCS heading when we lost one of the antenna
      if ddd+hhh/24>=(298+15/24) & ddd+hhh/24<=(300+16/24)
          jd_mot=ddd+(hhh+tmotion/3600)/24;
          lrgc=cos(lrg'*rdcon);        %convert to N
          lrgs=sin(lrg'*rdcon);        %convert to E
            lrgcm=interp1(jd_scs(~isnan(jd_scs)),lrgc(~isnan(jd_scs)),jd_mot); 
            lrgsm=interp1(jd_scs(~isnan(jd_scs)),lrgs(~isnan(jd_scs)),jd_mot); 
          heading=mod(atan2(lrgsm,lrgcm)/rdcon+360,360);         %laser ring gyro
      else
          read_head_pitch_HiWinGS_2013; % 
          heading=headx_pitch;
      end
      
      ii=(find(isnan(heading)));
      try
        heading(find(isnan(heading)))=heading(find(isnan(heading))+1);
        heading(find(isnan(heading)))=heading(find(isnan(heading))+1);
      catch
        if ii(1) == 1        
            ll=1;
            while   ii(1)==1
                heading(1)=heading(ii(ll)+1);
                ii=(find(isnan(heading)));
                ll=ll+1;
            end
        end    
        if ii(end) == 36000 & ~isempty(ii)
            ll=0;
            while   ii(end) == 36000
                heading(end)=heading(ii(end-ll)-1);
                ii=(find(isnan(heading)));
                if isempty(ii);    break;  end;                 
                ll=ll+1;
            end
        end  
        ii=(find(isnan(heading)));
        for ll=1:length(ii)
            heading(ii(ll))=heading(ii(ll)-1);  
        end    
%           disp('use meddian')        
      end
      shead=sin(heading./180.*pi);     %can use also unwrap function to avoid heading problem between 360 and 0 deg.
      chead=cos(heading./180.*pi);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

      licor_prefix='lic';
      read_co2_HiWinGS_2013;    %
      
%       %gets the surface current speed and direction from ADCP
%       ii=find(MMF(:,1)>=ddd+hhh/24-1/60/24 & MMF(:,1)<ddd+(hhh+1)/24);
%       Cspd=MMF(ii,7); %Current speed (m/s)
% %       Cdir=MMF(ii,8)+180; %Current direction (deg). +180 is to convert 'from' to 'to'
%        Cdir=MMF(ii,8); %Current direction (deg) 'from' -->needed that way in nau_correct.!
     
      nau_correct_HiWinGS_2013;
      nau_process_HiWinGS_2013;  
    end;
end;

% 
% for ddd=rge
%     jd=sprintf('%04i',ddd);
%     filecat([drive_w '\procturb\da_' jd '*' dcname '.prn'],[drive_w '\procturb\' cruise '_da' jd dcname '.txt']);
%     filecat([drive_w '\procturb\cr_' jd '*' dcname '.prn'],[drive_w '\procturb\' cruise '_cr' jd dcname '.txt']);
%     filecat([drive_w '\procturb\co_' jd '*' dcname '.prn'],[drive_w '\procturb\' cruise '_co' jd dcname '.txt']);
%     filecat([drive_w '\procturb\sp_' jd '*' dcname '.prn'],[drive_w '\procturb\' cruise '_sp' jd dcname '.txt']);
% end
% filecat([drive_w '\procturb\' cruise '_da*' dcname '.txt'],[drive_w '\procturb\' cruise '_da' datasetv dcname '.prn']);
% filecat([drive_w '\procturb\' cruise '_cr*' dcname '.txt'],[drive_w '\procturb\' cruise '_cr' datasetv dcname '.prn']);
% filecat([drive_w '\procturb\' cruise '_co*' dcname '.txt'],[drive_w '\procturb\' cruise '_co' datasetv dcname '.prn']);
% filecat([drive_w '\procturb\' cruise '_sp*' dcname '.txt'],[drive_w '\procturb\' cruise '_sp' datasetv dcname '.prn']);


% if ~exist([way_proc_data_flux 'procturb\' cruise '_da.txt'],'file')
%     filecat([way_proc_data_flux 'procturb\da_*.prn'],[way_proc_data_flux 'procturb\' cruise '_da_v1.txt']);
% end
% if ~exist([way_proc_data_flux 'procturb\' cruise '_cr.txt'],'file')
%     filecat([way_proc_data_flux 'procturb\cr_*.prn'],[way_proc_data_flux 'procturb\' cruise '_cr_v1.txt']);
% end
% if ~exist([way_proc_data_flux 'procturb\' cruise '_co.txt'],'file')
%     filecat([way_proc_data_flux 'procturb\co_*.prn'],[way_proc_data_flux 'procturb\' cruise '_co_v1.txt']);
% end
% if ~exist([way_proc_data_flux 'procturb\' cruise '_sp.txt'],'file')
%     filecat([way_proc_data_flux 'procturb\sp_*.prn'],[way_proc_data_flux 'procturb\' cruise '_sp_v1.txt']);
% end


% % calnex2010_da_combine;
% % calnex2010_cr_combine;
% % calnex2010_co_combine;
% % calnex2010_sp_combine;
% % 


% [SU,Fson_U]=psd2(detrend(U),36000,10,hamming(36000));[SU,Fson_U]=specsmoo(SU,10);
% [SV,Fson_V]=psd2(detrend(V),36000,10,hamming(36000));[SV,Fson_V]=specsmoo(SV,10);
% [SW,Fson_W]=psd2(detrend(W),36000,10,hamming(36000));[SW,Fson_W]=specsmoo(SW,10);
% figure;loglog(Fson_U,Fson_U.*SU,'b--',Fson_V,Fson_V.*SV,'r--',Fson_W,Fson_W.*SW,'g--')
% title('Normalized Spectra of raw U (blue), V (red), W (green)');


%check angles spectra
% % sf=sampleRate;accm=Plat;ratem=deg_rate;gyro=heading;its=5;
% % 
%***phi / roll
% % phiHF=filtfilt(bhi,ahi,atan2(accm(2,:),gravity));
% % [SphiHF,FphiHF]=psd2(detrend(phiHF),36000,10,hamming(36000));[SphiHF,FphiHF]=specsmoo(SphiHF,10);
% % 
% % phiraw=atan2(accm(2,:),gravity);
% % [Sphiraw,Fphiraw]=psd2(detrend(phiraw),36000,10,hamming(36000));[Sphiraw,Fphiraw]=specsmoo(Sphiraw,10);
% % 
% % [Sphislow,Fphislow]=psd2(detrend(phislow),36000,10,hamming(36000));[Sphislow,Fphislow]=specsmoo(Sphislow,10);
% % 
% % figure;loglog(Fphiraw,Fphiraw.*Sphiraw,'b',FphiHF,FphiHF.*SphiHF,'r--',Fphiraw,Fphiraw.*(Sphiraw-SphiHF),'g--');
% % title('Normalized Spectra of phi slow. Blue= raw, red= HF, green= BF')
% % 
% % phiint=filtfilt(bhi,ahi,((cumsum(ratem(1,:)) -.5*ratem(1,:) - .5*ratem(1,1))/sf));
% % [Sphiint,Fphiint]=psd2(detrend(phiint),36000,10,hamming(36000));[Sphiint,Fphiint]=specsmoo(Sphiint,10);
% % [Sphi,Fphi]=psd2(detrend(phi),36000,10,hamming(36000));[Sphi,Fphi]=specsmoo(Sphi,10);
% % figure;loglog(Fphi,Fphi.*abs(Sphi),'b',Fphiint,Fphiint.*abs(Sphiint),'r--',Fphislow,Fphislow.*abs(Sphislow),'g--');
% % title('Normalized Spectra of phi slow. Blue= euler angle computed, red= rate integration, green= BF')
% % 
% % figure;loglog(Fphiint,Fphiint.*abs(Sphiint),'r--',FphiHF,FphiHF.*abs(SphiHF),'r');
% % title('Normalized Spectra of phi HF. broken line= by integration of ratex, line=with gravity')
% % 
% % phiintnf=(cumsum(ratem(1,:)) -.5*ratem(1,:) - .5*ratem(1,1))/sf;
% % [Sphiintnf,Fphiintnf]=psd2(detrend(phiintnf),36000,10,hamming(36000));[Sphiintnf,Fphiintnf]=specsmoo(Sphiintnf,10);
% % figure;loglog(Fphiraw,Fphiraw.*Sphiraw,'b',Fphiintnf,Fphiintnf.*abs(Sphiintnf),'r--',Fphi,Fphi.*abs(Sphi),'g')
% % legend('accy./gravity', 'integral of angular rate x','Final euler roll angle computed')
% % 
% % 
% % 
% % % % %***theta / pitch
% % thetaHF=filtfilt(bhi,ahi,atan2(-accm(1,:),gravity));
% % [SthetaHF,FthetaHF]=psd2(detrend(thetaHF),36000,10,hamming(36000));[SthetaHF,FthetaHF]=specsmoo(SthetaHF,10);
% % 
% % thetaraw=atan2(-accm(1,:),gravity);
% % [Sthetaraw,Fthetaraw]=psd2(detrend(thetaraw),36000,10,hamming(36000));[Sthetaraw,Fthetaraw]=specsmoo(Sthetaraw,10);
% % 
% % [Sthetaslow,Fthetaslow]=psd2(detrend(thetaslow),36000,10,hamming(36000));[Sthetaslow,Fthetaslow]=specsmoo(Sthetaslow,10);
% % 
% % figure;loglog(Fthetaraw,Fthetaraw.*Sthetaraw,'b',FthetaHF,FthetaHF.*SthetaHF,'r--',Fthetaraw,Fthetaraw.*(Sthetaraw-SthetaHF),'g--');
% % title('Normalized Spectra of theta slow. Blue= raw, red= HF, green= BF')
% % 
% % thetaint=filtfilt(bhi,ahi,((cumsum(ratem(2,:)) -.5*ratem(2,:) - .5*ratem(2,1))/sf));
% % [Sthetaint,Fthetaint]=psd2(detrend(thetaint),36000,10,hamming(36000));[Sthetaint,Fthetaint]=specsmoo(Sthetaint,10);
% % [Stheta,Ftheta]=psd2(detrend(theta),36000,10,hamming(36000));[Stheta,Ftheta]=specsmoo(Stheta,10);
% % figure;loglog(Ftheta,Ftheta.*abs(Stheta),'b',Fthetaint,Fthetaint.*abs(Sthetaint),'r--',Fthetaslow,Fthetaslow.*abs(Sthetaslow),'g--');
% % title('Normalized Spectra of theta slow. Blue= euler angle computed, red= rate integration, green= BF')
% % 
% % figure;loglog(Fthetaint,Fthetaint.*abs(Sthetaint),'r--',FthetaHF,FthetaHF.*abs(SthetaHF),'r');
% % title('Normalized Spectra of theta HF. broken line= by integration of ratey, line=with gravity')
% % 
% % thetaintnf=(cumsum(ratem(2,:)) -.5*ratem(2,:) - .5*ratem(2,1))/sf;
% % [Sthetaintnf,Fthetaintnf]=psd2(detrend(thetaintnf),36000,10,hamming(36000));[Sthetaintnf,Fthetaintnf]=specsmoo(Sthetaintnf,10);
% % figure;loglog(Fthetaraw,Fthetaraw.*Sthetaraw,'b',Fthetaintnf,Fthetaintnf.*abs(Sthetaintnf),'r--',Ftheta,Ftheta.*abs(Stheta),'g')
% % legend('-accx./gravity', 'integral of angular rate y','Final euler pitch angle computed')

% % %****psi  / yaw
% % psiHF=filtfilt(bhi,ahi,-gyro);
% % [SpsiHF,FpsiHF]=psd2(detrend(psiHF),36000,10,hamming(36000));[SpsiHF,FpsiHF]=specsmoo(SpsiHF,10);
% % 
% % psiraw=-gyro;
% % [Spsiraw,Fpsiraw]=psd2(detrend(psiraw),36000,10,hamming(36000));[Spsiraw,Fpsiraw]=specsmoo(Spsiraw,10);
% % 
% % [Spsislow,Fpsislow]=psd2(detrend(psislow),36000,10,hamming(36000));[Spsislow,Fpsislow]=specsmoo(Spsislow,10);
% % 
% % figure;loglog(Fpsiraw,Fpsiraw.*Spsiraw,'b',FpsiHF,FpsiHF.*SpsiHF,'r--',Fpsiraw,Fpsiraw.*(Spsiraw-SpsiHF),'g--');
% % title('Normalized Spectra of psi slow. Blue= raw, red= HF, green= BF')
% % 
% % psiint=filtfilt(bhi,ahi,((cumsum(ratem(3,:)) -.5*ratem(3,:) - .5*ratem(3,1))/sf));
% % [Spsiint,Fpsiint]=psd2(detrend(psiint),36000,10,hamming(36000));[Spsiint,Fpsiint]=specsmoo(Spsiint,10);
% % [Spsi,Fpsi]=psd2(detrend(psi),36000,10,hamming(36000));[Spsi,Fpsi]=specsmoo(Spsi,10);
% % figure;loglog(Fpsi,Fpsi.*abs(Spsi),'b',Fpsiint,Fpsiint.*abs(Spsiint),'r--',Fpsislow,Fpsislow.*abs(Spsislow),'g--');
% % title('Normalized Spectra of psi slow. Blue= euler angle computed, red= rate integration, green= BF')
% % 
% % figure;loglog(Fpsiint,Fpsiint.*abs(Spsiint),'r--',FpsiHF,FpsiHF.*abs(SpsiHF),'r');
% % title('Normalized Spectra of psi HF. broken line= by integration of ratez, line=with gravity')
% % 
% % psiintnf=(cumsum(ratem(3,:)) -.5*ratem(3,:) - .5*ratem(3,1))/sf;
% % [Spsiintnf,Fpsiintnf]=psd2(detrend(psiintnf),36000,10,hamming(36000));[Spsiintnf,Fpsiintnf]=specsmoo(Spsiintnf,10);
% % figure;loglog(Fpsiraw,Fpsiraw.*Spsiraw,'b',Fpsiintnf,Fpsiintnf.*abs(Spsiintnf),'r--',Fpsi,Fpsi.*abs(Spsi),'g')
% % legend('Gyrocompass', 'integral of angular rate z','Final euler yaw angle computed')
% % 


% % % % % % %more plots
% % wp=1/240/(sf/2);
% % ws=1/120/(sf/2);
% % [n,wn]=buttord(wp,ws,1,14);
% % [b,a]=butter(n,wn,'high');
% % 
% % wp=1/120/(sf/2);
% % ws=.8*wp;
% % [n,wn]=buttord(wp,ws,1,14);
% % % % %disp([num2str(ws),'  ',num2str(wn),'  ',num2str(wp)])
% % [b,a]=butter(n,wn,'high');
% % 
% % wp=1/25/(sf/2);
% % ws=0.8*wp;
% % [n,wn]=buttord(wp,ws,3,7);
% % [b,a]=butter(n,wn,'high');
% % 
% % [h,f]=freqz(b,a,36000,10);
% % figure;semilogx(f,h,'k')
% % 
% % % % 
% [Sx,Fx]=psd2(detrend(uvwplat(1,:)),36000,10,hamming(36000));[Sx,Fx]=specsmoo(Sx,10);
% [Sy,Fy]=psd2(detrend(uvwplat(2,:)),36000,10,hamming(36000));[Sy,Fy]=specsmoo(Sy,10);
% [Sz,Fz]=psd2(detrend(uvwplat(3,:)),36000,10,hamming(36000));[Sz,Fz]=specsmoo(Sz,10);
% uvwplatf(1,:) = filtfilt(b,a,uvwplat(1,:));
% uvwplatf(2,:) = filtfilt(b,a,uvwplat(2,:));
% uvwplatf(3,:) = filtfilt(b,a,uvwplat(3,:));
% [Sxf,Fxf]=psd2(detrend(uvwplatf(1,:)),36000,10,hamming(36000));[Sxf,Fxf]=specsmoo(Sxf,10);
% [Syf,Fyf]=psd2(detrend(uvwplatf(2,:)),36000,10,hamming(36000));[Syf,Fyf]=specsmoo(Syf,10);
% [Szf,Fzf]=psd2(detrend(uvwplatf(3,:)),36000,10,hamming(36000));[Szf,Fzf]=specsmoo(Szf,10);
% figure;loglog(Fx,Fx.*Sx,'b',Fxf,Fxf.*Sxf,'b--')
% figure;loglog(Fy,Fy.*Sy,'b',Fyf,Fyf.*Syf,'b--')
% figure;loglog(Fz,Fz.*Sz,'b',Fzf,Fzf.*Szf,'b--')
% 
% 
% 
% % [Sx,Fx]=psd2(detrend(AAAA(1,:)),36000,10,hamming(36000));[Sx,Fx]=specsmoo(Sx,10);
% % [Sy,Fy]=psd2(detrend(AAAA(2,:)),36000,10,hamming(36000));[Sy,Fy]=specsmoo(Sy,10);
% % [Sz,Fz]=psd2(detrend(AAAA(3,:)),36000,10,hamming(36000));[Sz,Fz]=specsmoo(Sz,10);
% % [Sxc,Fxc]=psd2(detrend(uvwplat(1,:)),36000,10,hamming(36000));[Sxc,Fxc]=specsmoo(Sxc,10);
% % [Syc,Fyc]=psd2(detrend(uvwplat(2,:)),36000,10,hamming(36000));[Syc,Fyc]=specsmoo(Syc,10);
% % [Szc,Fzc]=psd2(detrend(uvwplat(3,:)),36000,10,hamming(36000));[Szc,Fzc]=specsmoo(Szc,10);
% % [Sxf,Fxf]=psd2(detrend(DDDD(:,1)),36000,10,hamming(36000));[Sxf,Fxf]=specsmoo(Sxf,10);
% % [Syf,Fyf]=psd2(detrend(DDDD(:,2)),36000,10,hamming(36000));[Syf,Fyf]=specsmoo(Syf,10);
% % [Szf,Fzf]=psd2(detrend(DDDD(:,3)),36000,10,hamming(36000));[Szf,Fzf]=specsmoo(Szf,10);
% % figure;loglog(Fx,Fx.*Sx,'b',Fxc,Fxc.*Sxc,'r',Fxf,Fxf.*Sxf,'g')
% % figure;loglog(Fy,Fy.*Sy,'b',Fyc,Fyc.*Syc,'r',Fyf,Fyf.*Syf,'g')
% % figure;loglog(Fz,Fz.*Sz,'b',Fzc,Fzc.*Szc,'r',Fzf,Fzf.*Szf,'g')
% 
% lu=find(Fzc.*Szc==max(Fzc.*Szc))
% figure;loglog(Fz(lu-2:lu+2),Szc(lu-2:lu+2)./Sz(lu-2:lu+2),'b')
% correctionratio=nanmean(Szc(lu-2:lu+2)./Sz(lu-2:lu+2));
% 
% % figure;loglog(Fz,Fz.*Sz,'b',Fzc,Fzc.*Szc,'r',Fzf,Fzf.*Szf,'g',Fzc,Fzc.*Szc/correctionratio,'m')
% figure;loglog(Fx,Fx.*Sx,'b',Fxc,Fxc.*Sxc,'r',Fxf,Fxf.*Sxf,'g',Fxc,Fxc.*Sxc/correctionratiox,'m')
% 
% lu=find(Fxc.*Sxc==max(Fxc.*Sxc))
% correctionratiox=nanmean(Sxc(lu-2:lu+2)./Sx(lu-2:lu+2));
% lu=find(Fyc.*Syc==max(Fyc.*Syc))
% correctionratioy=nanmean(Syc(lu-2:lu+2)./Sy(lu-2:lu+2));
% lu=find(Fzc.*Szc==max(Fzc.*Szc))
% correctionratioz=nanmean(Szc(lu-2:lu+2)./Sz(lu-2:lu+2))
% 
% uvwplatc=ones(3,36000);
% uvwplatc(1,:) = ifft(Sxc/correctionratiox);
% 
% uvwplatc=ones(3,36000);
% Fs = 10;                    % Sampling frequency
% ft = Fs*linspace(0,1,NFFT);
% NFFT = 36000; % Next power of 2 from length of y
% YU = fft(AAAA(1,:),NFFT);YUc = fft(CCCC(1,:),NFFT);
% YV = fft(AAAA(2,:),NFFT);YVc = fft(CCCC(2,:),NFFT);
% YW = fft(AAAA(3,:),NFFT);YWc = fft(CCCC(3,:),NFFT);
% 
% titi=find(ft>5e-2 & ft<3e-1);
% lu=find(abs(YUc(titi))>abs(YU(titi)));
% correctionratiox=nanmean(abs(YUc(titi(lu)))./abs(YU(titi(lu))));
% lu=find(abs(YVc(titi))>abs(YV(titi)));
% correctionratioy=nanmean(abs(YVc(titi(lu)))./abs(YV(titi(lu))));
% lu=find(abs(YWc(titi))>abs(YW(titi)));
% correctionratioz=nanmean(abs(YWc(titi(lu)))./abs(YW(titi(lu))));
% 
% %does not allow to overcorrect spectra...
% YU2=YUc;
% if correctionratiox>1
%     YU2(titi(lu))=(YUc(titi(lu)))./correctionratiox;
% end
% YV2=YVc;
% if correctionratioy>1
%     YV2(titi(lu))=(YVc(titi(lu)))./correctionratioy;
% end;
% YW2=YWc;
% if correctionratioz>1
%     YW2(titi(lu))=(YWc(titi(lu)))./correctionratioz;
% end;
% uvwplatc(1,:) = real(ifft(YU2));
% uvwplatc(2,:) = real(ifft(YV2));
% uvwplatc(3,:) = real(ifft(YW2));


% % 
% % 
% % fs=10;         %Sonic are output at 20 Hz
% % fc=1/25;       %Cutoff frequency for angle calculation
% % nfreq = fs/2.0;                   %Nyquist frequency   
% % wp = fc/nfreq;          
% % ws = .7*wp;   
% % [n,wn] = buttord(wp,ws,10,25);        %Filter coefficients
% % [bhi,ahi] = butter(n,wn,'high');