% mmp_pgrid_mfile.m
%  a routine to reduce and bin data from the McLane Moored Profiler
%  copyright John M. Toole
%  Woods Hole Oceanographic Institution
%
%  December, 2001; February 2002; 
%  July 2002 - change CTD filtering to avoid ringing at sharp steps and
%              implement Lueck's lag correction routine

disp([' Program mmp_pgrid version 1.2 - July 2002'])
%disp([' PROGRAM mmp_pgrid_mfile MFC  version 2 - Mar 2002'])
%disp(['    (It was mmp_pgrid       JMT  version 1 - Dec 2001)'])

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


%begin with some preliminaries 

   if ~exist('MMP_id')
    MMP_id=input('To begin, please enter the ID of the MMP you are processing: ','s');
   end

   if ~exist('CTD_id')
	CTD_id=input('    What was the serial number of the CTD that was on this Profiler? ','s');
   end

   if ~exist('ACM_id')
	ACM_id=input('    And what was the serial number of the ACM that was on it? ','s');
   end

   if ~exist('experiment_name')
    experiment_name=input('Please enter the experiment name: ','s'); 
   end

   if ~exist('mooring_position')
    mooring_position=input('    What was the position of the mooring: [lat,lon] ');
   end

   if ~exist('mag_dev')
	    [mag_dev]=input('Please enter the magnetic deviation for this site in degrees west of north ');
   end
	    mag_dev=mag_dev*pi/180.;

   disp([' ']);
   if ~exist('stas_text')
    stas_text=input('Which stations do you want to process? eg. [8:24]  ','s');
   end
   eval(['stas=',stas_text,';']);

   if ~exist('directory')
    [directory]=input('Which directory holds the raw, Matlab-format data files :   ','s');
   end

   if ~exist('outdir')
    disp([sprintf('\n'),'Into which directory do you want to put ']);
    [outdir]=input(' the output data and documentation file? :  ','s');
   end

    docname=['mmp_pgrid_',num2str(stas(1)),'_',num2str(stas(end))];
    docfilename=[outdir  slsh docname '.doc'];
 	fid=fopen(docfilename,'w');
	  count=fprintf(fid,' \n Documentation file for a run of program mmp_pgrid_mfile.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',directory);
	  count=fprintf(fid,' Output data directory = %s\n',outdir);
	  
	  
   if ~exist('gmin')
    [gmin]=input('Please enter the Pressure Grid minimum value:  ');
   end

   if ~exist('gmax')
    [gmax]=input('  Now the Pressure Grid maximum value:  ');
   end

   if ~exist('gstep')
    [gstep]=input('    Lastly, the Grid step: ');
   end

   count=fprintf(fid,' \n Grid min, max and step = %g %g %g \n',[gmin gmax gstep]);

   %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']);

   if ~exist('nscans_test')
     nscans_test=input('Enter the minimum number of scans you require to process a profile '); 
   else
     disp([' You have set the requirement of  ',num2str(nscans_test),' scans in order to process']);
   end

%now to enter the file names that hold sensor calibration information for this deployment
    
   disp([' ']);
   if ~exist('ctd_cal_filename')
    ctd_cal_filename=input('Enter the full name (with path) of the matlab file holding the CTD calibration information: ','s');
   end
   if exist(ctd_cal_filename)~=2
      disp([' File does not exist, or cannot be opened, try again ']);
      ctd_cal_filename=input('Enter the full name (with path) of the matlab file holding the CTD calibration information: ','s');
      if exist(ctd_cal_filename)~=2
	    return
      end
   end

   disp([' ']);
   if ~exist('acm_cal_filename')
    acm_cal_filename=input('Enter the full name (with path) of the matlab file holding the ACM calibration information: ','s');
   end
   if exist(acm_cal_filename)~=2
    disp([' File does not exist, or cannot be opened, try again ']);
    acm_cal_filename=input('Enter the full name (with path) of the matlab file holding the ACM calibration information: ','s');
    if exist(acm_cal_filename)~=2
 	    return
    end
   end

   eval(['load ',ctd_cal_filename]);
   eval(['load ',acm_cal_filename]);

   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=[directory,slsh,'raw',sprintf('%4.4i',stas(h)),'.mat'];
    disp([' ']);
    disp(['   working on file ',rawmatfile]);
    disp([' ']);

          if( h==1);  count=fprintf(fid,' \n'); end;
	  count=fprintf(fid,'Processing file %s ',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);
    
    
    %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);
    [bb,aa]=butter(1,.25);
    [cpres]=myfiltfilt(bb,aa,cpres);
    [ctemp]=myfiltfilt(bb,aa,ctemp);
    [ccond]=myfiltfilt(bb,aa,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 ',matfmtfile,' has a problem:  file skipped'])
	  count=fprintf(fid,'   **PROBLEM   FILE SKIPPED** \n');

    end
    
    end	%end looping on files

	fclose(fid);
