% mag_pgrid.m
% a routine to reduce and bin data from the McLane Moored Profiler
%  copyright John M. Toole
%  Woods Hole Oceanographic Institution
%  December, 2001

disp([' PROGRAM mmp_pgrid  MFC  version 2 - Mar 2002'])
disp(['    (It was mmp_pgrid       JMT  version 1 - Dec 2001)'])
disp([' *******************************************************'])

%begin with some preliminaries 
% get the input info from the gui window 
%   each variable has a discrete TAG, and everything that's input from an edit 
% window has ED as the suffix.

tmp=findobj(gcf,'Tag','MMPidEd');
  MMP_id=char(get(tmp,'String'));
tmp=findobj(gcf,'Tag','CTDidEd');
  CTD_id=char(get(tmp,'String'));
tmp=findobj(gcf,'Tag','ACMidEd');
  ACM_id=char(get(tmp,'String'));
tmp=findobj(gcf,'Tag','exnamEd');
  experiment_name=char(get(tmp,'String'));
tmp=findobj(gcf,'Tag','inEd');
  direc=char(get(tmp,'String'));
tmp=findobj(gcf,'Tag','outEd');
  outdir=char(get(tmp,'String'));
tmp=findobj(gcf,'Tag','posEd');
  mooring_position=char(get(tmp,'String'));
tmp=findobj(gcf,'Tag','declEd');
  mag_dev=str2num(char(get(tmp,'String')));
    mag_dev=mag_dev*pi/180.;
tmp=findobj(gcf,'Tag','drngEd');
 dat=str2num(char(get(tmp,'String')));
 stas_text=char(get(tmp,'String'));
    eval(['stas=[' num2str(dat(1)),':1:',num2str(dat(2)),'];']);
tmp=findobj(gcf,'Tag','pgrEd');
  pgs=str2num(char(get(tmp,'String')));
    gmin=pgs(1);
    gstep=pgs(2);
    gmax=pgs(3);
% now the cal files
tmp=findobj(gcf,'Tag','CcalEd');
  ctd_cal_filename=char(get(tmp,'String'));
tmp=findobj(gcf,'Tag','AcalEd');
  acm_cal_filename=char(get(tmp,'String'));
% now the bounds
tmp=findobj(gcf,'Tag','pbndEd');
  p_bounds=str2num(char(get(tmp,'String')));
tmp=findobj(gcf,'Tag','tbndEd');
  t_bounds=str2num(char(get(tmp,'String')));
tmp=findobj(gcf,'Tag','cbndEd');
  c_bounds=str2num(char(get(tmp,'String')));
tmp=findobj(gcf,'Tag','vbndEd');
  vel_bounds=str2num(char(get(tmp,'String')));
  compass_bounds=[-1.1 1.1];
tmp=findobj(gcf,'Tag','nscnEd');
   nscans_test=str2num(char(get(tmp,'String')));

% assign the proper sort of slash, for PC & Unix
arch=computer;
if (strcmp(arch,'PCWIN')) 
   slsh='\';
else
   slsh='/';
end

%  write a header on to the output file
    docname=['mmp_pgrid_',num2str(stas(1)),'_',num2str(stas(end))];
      if (exist(outdir,'dir'))
        docfilename=[outdir slsh docname '.doc'];
      else
        msg='Output path does not exist.  Re-enter the correct path';
        msg_window(msg);
         while (findobj('type','figure','tag','msgwin'))
          pause(2)
         end
	 break;
      end

     docfilename=[outdir slsh docname '.doc'];
 	fid=fopen(docfilename,'w');

	  count=fprintf(fid,' \n Documentation file for a run of program mmp_pgrid.m version %g.\n',2.0);
          count=fprintf(fid,' Date of Pressure Gridding:  %s\n',date);
 	  count=fprintf(fid,' \n Experiment name %s\n',experiment_name);
	  count=fprintf(fid,'\n Mooring position = %s\n',mooring_position);
	  count=fprintf(fid,'\n MMP identification number = %s\n',MMP_id);
	  count=fprintf(fid,'     with CTD identification number %s\n',CTD_id);
	  count=fprintf(fid,'     and  ACM identification number %s\n',ACM_id);
	  count=fprintf(fid,'\n Magnetic deviation used in this run is %g.\n',mag_dev*180/pi);
	  count=fprintf(fid,'\n This run is for stations %s\n',stas_text);
	  count=fprintf(fid,'\n Input data directory = %s\n',direc);
	  count=fprintf(fid,' Output data directory = %s\n',outdir);
	  
	  

   %use the grid bounds to estimate how many CTD and ACM data scans you expect from a full profile
   nscans_expected=fix((gmax-gmin)/.25*1.5)+400;
   disp([' Based on the gridding, we expect ',num2str(nscans_expected),' scans in a data file']);

%now to enter the file names that hold sensor calibration information for this deployment
 % load the cal files     
      if exist(ctd_cal_filename,'file')
	eval(['load ',ctd_cal_filename]);
      else
        msg='Could not find CTD cal file.  Re-enter the correct path and filename';
        msg_window(msg);
         while (findobj('type','figure','tag','msgwin'))
          pause(2)
         end
	 break;
       end

      if exist(acm_cal_filename,'file')
	eval(['load ',acm_cal_filename]);
      else
        msg='Could not find ACM cal file.  Re-enter the correct path and filename';
        msg_window(msg);
         while (findobj('type','figure','tag','msgwin'))
          pause(2)
         end
	 break;
       end

   count=fprintf(fid,'\n CTD calibration filename = %s\n',ctd_cal_filename);
   count=fprintf(fid,'     Pres coeffs = %g  %g \n',p_coefs);
   count=fprintf(fid,'     Temp coeffs = %g %g \n',t_coefs);
   count=fprintf(fid,'     Cond coeffs = %g  %g \n',c_coefs);
   count=fprintf(fid,'     Temp Lag    = %g \n',tlag);

   count=fprintf(fid,'\n ACM calibration filename = %s\n',acm_cal_filename);
   count=fprintf(fid,'     Hx_bias  = %g \n',Hx_bias);
   count=fprintf(fid,'     Hx_range = %g \n',Hx_range);
   count=fprintf(fid,'     Hy_bias  = %g \n',Hy_bias);
   count=fprintf(fid,'     Hy_range = %g \n',Hy_range);
   count=fprintf(fid,'     Vab_bias = %g \n',Vab_bias);
   count=fprintf(fid,'     Vcd_bias = %g \n',Vcd_bias);
   count=fprintf(fid,'     Vef_bias = %g \n',Vef_bias);
   count=fprintf(fid,'     Vgh_bias = %g \n',Vgh_bias);
   count=fprintf(fid,'\n     Compass bias   = %g \n',compass_bias);
   count=fprintf(fid,'     Direction sign = %g \n',dir_sign);
   count=fprintf(fid,'     Velocity Scale = %g \n',velocity_scale);



   if ~exist('p_bounds') | ~exist('t_bounds') | ~exist('c_bounds') | ~exist('vel_bounds')

      disp([' ']);
      disp([' Now you must enter sensor edit parameters: ']);
      disp(['     For each raw parameter, enter the minimum and maximum acceptable value as a two-element vector, e.g. [2 24] '])
      disp([' ']);

	p_bounds=input(' Pressure bounds: ');
	t_bounds=input(' Temperature bounds: ');
	c_bounds=input(' Conductivity bounds: ');
	vel_bounds=input(' Velocity bounds: ');
   end
   compass_bounds=[-1.1 1.1];

   count=fprintf(fid,'\n Pressure edit bounds     = %g %g \n',p_bounds);
   count=fprintf(fid,' Temperature edit bounds  = %g %g \n',t_bounds);
   count=fprintf(fid,' Conductivity edit bounds = %g %g \n',c_bounds);
   count=fprintf(fid,' Velocity edit bounds     = %g %g \n',vel_bounds);
   count=fprintf(fid,' Compass edit bounds      = %g %g \n',compass_bounds);

%now start the processing

    for h=1:length(stas)
    
    rawmatfile=[direc slsh,'raw',sprintf('%4.4i',stas(h)),'.mat'];
    disp([' ']);
    disp(['   working on file ',rawmatfile]);

	  count=fprintf(fid,' \n');
	  count=fprintf(fid,'Processing file %s \n',rawmatfile);

   if exist(rawmatfile)
    [okay Vab Vcd Vef Vgh aHx aHy aHz aTx aTy cpres ctemp ccond ...
	 psdate pstart pedate pstop epres ecurr edpdt evolt engtime]...
	    =load_MMP_rawmatfmt(rawmatfile,nscans_test);
      else
        msg=['Could not find ' rawmatfile '.  Re-enter the correct path and filename'];
        msg_window(msg);
         while (findobj('type','figure','tag','msgwin'))
          pause(2)
         end
	 break;
       end
    
    
    %check that this file is okay, if not skip this profile
    if (okay==1);


%extract the profile start and stop times 

    j=find(psdate=='/');
    psdate(j)=' ';
    date=str2num(psdate);
    
    mon=date(1);
    day=date(2);   
    year=2000+date(3);
        
    j=find(pstart==':');
    pstart(j)=' ';
    start=str2num(pstart);
    
    hour=start(1);
    minute=start(2);
    second=start(3);
    
    startdaytime=datenum(year,mon,day,hour,minute,second);
    
    j=find(pedate=='/');
    pedate(j)=' ';
    date=str2num(pedate);
    
    mon=date(1);
    day=date(2);
    year=2000+date(3);
        
    j=find(pstop==':');
    pstop(j)=' ';
    stop=str2num(pstop);
    
    hour=stop(1);
    minute=stop(2);
    second=start(3);
    
    stopdaytime=datenum(year,mon,day,hour,minute,second);

%next we create a time variable for each CTD data scan using the engineering data epres and etime.

    jj=find(epres~=0);
    nep=jj(1);
    mep=jj(end);
    diffp=abs(epres(nep)-cpres);
    mindifp=min(diffp);
    jj=find(diffp==mindifp);
    ncp=jj(1);
    diffp=abs(epres(mep)-cpres);
    mindifp=min(diffp);
    jj=find(diffp==mindifp);
    mcp=jj(end);
    
    ctime=[1:length(cpres)];
    ctime=ctime-ctime(ncp);
    ctdsamplerate=(engtime(mep)-engtime(nep))/(mcp-ncp);
    ctime=(engtime(nep)+ctime*ctdsamplerate)';



% now we edit the CTD and ACM data

    %first, we tag ACM path current data very near zero as bad.
	%this is done to eliminate some data spikes occasionally seen in the data
	
	j=find(abs(Vab)<0.05);
	if(length(j)>0)
	Vab(j)=-1000;
	end
	
	j=find(abs(Vcd)<0.05);
	if(length(j)>0)
	Vcd(j)=-1000;
	end
	
	j=find(abs(Vef)<0.05);
	if(length(j)>0)
	Vef(j)=-1000;
	end
	
	j=find(abs(Vgh)<0.05);
	if(length(j)>0)
	Vgh(j)=-1000;
	end
	 
    [cpres]=MMP_edit(p_bounds,3,cpres);
    [ctemp]=MMP_edit(t_bounds,3,ctemp);
    [ccond]=MMP_edit(c_bounds,3,ccond);
    [Vab]=MMP_edit(vel_bounds,11,Vab);
    [Vcd]=MMP_edit(vel_bounds,11,Vcd);
    [Vef]=MMP_edit(vel_bounds,11,Vef);
    [Vgh]=MMP_edit(vel_bounds,11,Vgh);
    [aHx]=MMP_edit(compass_bounds,11,aHx);
    [aHy]=MMP_edit(compass_bounds,11,aHy);   
%    [epres]=MMP_edit(p_bounds,1,epres);


%here we apply calibrations to the CTD and ACM sensor data
    
    cpres=polyval(p_coefs,cpres);
    ctemp=polyval(t_coefs,ctemp);
    ccond=polyval(c_coefs,ccond);
    
    Vab=(Vab-Vab_bias)*velocity_scale;
    Vcd=(Vcd-Vcd_bias)*velocity_scale;
    Vef=(Vef-Vef_bias)*velocity_scale;
    Vgh=(Vgh-Vgh_bias)*velocity_scale;
    
    aHx=(aHx-Hx_bias)/Hx_range;
    aHy=(aHy-Hy_bias)/Hy_range;

% and recursively filter the conductivity and pressure to match the temperature

    [ccond]=MMP_recur_filt(ccond,tlag);
    [cpres]=MMP_recur_filt(cpres,tlag);


% here we lowpass filter the pres, temp & cond, Vab, Vcd, Vef, Vgh, aHx, and aHy data

    [b,a]=butter(8,.08);
    [cpres]=myfiltfilt(b,a,cpres);
    [ctemp]=myfiltfilt(b,a,ctemp);
    [ccond]=myfiltfilt(b,a,ccond);
    [Vab]=myfiltfilt(b,a,Vab);
    [Vcd]=myfiltfilt(b,a,Vcd);
    [Vef]=myfiltfilt(b,a,Vef);
    [Vgh]=myfiltfilt(b,a,Vgh);
    [aHx]=myfiltfilt(b,a,aHx);
    [aHy]=myfiltfilt(b,a,aHy);
%    [epres]=myfiltfilt(b,a,epres);

% Now calculate the velocity in instrument and geographic coordinates

   % first, derive the velocity data in Cartesian instrument coordinates
	
	Vx=-(Vab+Vef)/(2.*.707);
	Vy=(Vab-Vef)/(2.*.707);
	Vz1=Vx-Vgh/.707;
	Vz2=-Vx+Vcd/.707;

	
   % now select the vertical velocity channel that is "upstream"
	
	avevz=mean(Vz1);
	if(avevz<0);
	Vz=Vz1;
	else
	Vz=Vz2;
	end
    

    % here we normalize the compass data onto a unit circle, and derive the compass heading of the ACM sting

	Hmag=sqrt(aHx.*aHx+aHy.*aHy);
	aHx=aHx./Hmag;
	aHy=aHy./Hmag;
    	mvpdir=dir_sign*atan2(aHx,aHy)+compass_bias+mag_dev;

%note to me: here is where I might try a correction to the velocity data for the wagging

    %now derive East and North velocity components

	speed=sqrt(Vx.*Vx+Vy.*Vy);
	reldir=atan2(Vy,Vx);
	trudir=mvpdir+reldir;
	Veast=speed.*cos(trudir);
	Vnorth=speed.*sin(trudir);

%now we start the pressure binning 

    %first find the scan numbers for the ctd and acm at profile start and end
    
    	[startc,endc,starta,enda,dpdt]=MMP_align_ctdacm(cpres(:),Vz,ctdsamplerate);
    
    %find the CTD scans that correspond to each pressure interval
    
	cscan1=[];cscan2=[];
	[pgrid,cscan1,cscan2]=pbins(gmin,gmax,gstep,cpres);
    
    %now average everything inside each pressure bin
    
    pave=[];tave=[];cave=[];uave=[];vave=[];wave=[];dpdtave=[];ascan1=[];ascan2=[];
    ctimave=[];
        
    for j=1:length(pgrid);
    
    %first check if we have any data in a given bin,  if no, set all = nan;
 
   if isnan(cscan1(j))==1 | isnan(cscan2(j))==1;
    
    pave(j)=nan;
    tave(j)=nan;
    cave(j)=nan;
    uave(j)=nan;
    vave(j)=nan;
    wave(j)=nan;
    dpdtave(j)=nan;
    ascan1(j)=nan;
    ascan2(j)=nan;
    ctimave(j)=nan;
    
    else;
    
    pave(j)=mean(cpres(cscan1(j):cscan2(j)));
    tave(j)=mean(ctemp(cscan1(j):cscan2(j)));
    cave(j)=mean(ccond(cscan1(j):cscan2(j)));
    dpdtave(j)=mean(dpdt(cscan1(j):cscan2(j)));
    ctimave(j)=mean(ctime(cscan1(j):cscan2(j)));
    
    %next derive the acm scans corresponding to the ctd scan interval
    
    ascan1(j)=fix(starta+(cscan1(j)-startc)*(enda-starta)/(endc-startc));
    ascan2(j)=fix(starta+(cscan2(j)-startc)*(enda-starta)/(endc-startc));
    
    if(ascan1(j)<1);
    ascan1(j)=1;
    end
    if ascan1(j)>length(Vz);
    ascan1(j)=length(Vz);
    end
    if(ascan2(j)<1);
    ascan2(j)=1;
    end
    if ascan2(j)>length(Vz);
    ascan2(j)=length(Vz);
    end
    
    uave(j)=mean(Veast(ascan1(j):ascan2(j)));
    vave(j)=mean(Vnorth(ascan1(j):ascan2(j)));
    wave(j)=mean(Vz(ascan1(j):ascan2(j)));
    
    
  end;  %if ;
    
  end;  %j loop;

    
    % calculate derived CTD quantities 
    
    CR=cave/sw_c3515.*(1-(6.5e-6).*(tave-2.) + (1.5e-8).*(pave-2000.));
    s_ave=sw_salt(CR,tave,pave);
    thetave=sw_ptmp(s_ave,tave,pave,0);
    sigthave=sw_pden(s_ave,tave,pave,0)-1000.;
        
%lastly write the binned data to file

    eval(['save ', outdir slsh,'grd',sprintf('%4.4i',stas(h)),'.mat'...
     ' pgrid pave  tave cave s_ave thetave sigthave uave vave wave dpdtave'...
     ' cscan1 cscan2 ascan1 ascan2 startdaytime stopdaytime psdate pstart' ...
     ' pedate pstop ctimave'])

	  count=fprintf(fid,'   OKAY \n');
    
    else	%if data file okay loop
    disp(['File ',rawmatfile,' skipped- try with a smaller min scans to process value'])
	  count=fprintf(fid,'   **PROBLEM   FILE SKIPPED** \n');

    end
    
    end	%end looping on files

	fclose(fid);
