%scale_rotate_hilo.m
%
% Now we will read in each profile that has dynamic height and
% Do scaling and rotating a la John 
% Do high and low pass filtering of uave,vave,dynh 
% using john's new mirrorfiltfilt
% overwrite dyn files with new ones
%
%2 May 2001
%
% basic input information

clear all

[direc]=['/j2/data/mvp/weddell_9900/grdfmt/dynh/'];
stas=input('Which stations do you want to process?   ');
disp([sprintf('\n'),'input stations will be :  ',direc,'dyn*[',num2str(min(stas)),':',num2str(max(stas)),']','.mat',sprintf('\n')]);

for i=1:length(stas)

% load input file

  matfile=([direc,'dyn',sprintf('%4.4i',stas(i)),'.mat']);

  if(exist(matfile)==2)
    disp([sprintf('\n'),'   working on file ',matfile]);

    eval(['load ',matfile]);
    clear matfmtfile

% scale and rotate a la John

    speed=1.39*sqrt(uave.*uave+vave.*vave);
    dir=atan2(vave,uave)-0.384;
    uave=speed.*cos(dir);
    vave=speed.*sin(dir);
    wave=wave*1.39;

%  edit last non nan value of uave, vave, wave to be equal to previous

    jj=find(isnan(wave));
    uave(jj(1)-1)=uave(jj(1)-2);
    vave(jj(1)-1)=vave(jj(1)-2);
    wave(jj(1)-1)=wave(jj(1)-2);

    clear jj;

% do highpass filters of u,v,w and dynh
% this is actually hipass2.m encorporated here.

    xx=find(~isnan(uave));

    ulow=uave;vlow=vave;dlow=dynh;
    [b,a]=mybutter(4,1/1000,.5);
    ulow(xx)=(mirrorfiltfilt(b,a,uave(xx)'));
    uhigh=(uave-ulow);
    vlow(xx)=(mirrorfiltfilt(b,a,vave(xx)'));
    vhigh=(vave-vlow);
    dlow(xx)=(mirrorfiltfilt(b,a,dynh(xx)))';
    dhigh=(dynh-dlow)';

 end   % if exist

%  save output

   eval(['save ',matfile]);
   clear speed dir uave vave wave ulow uhigh vlow vhigh dlow dhigh;

end  %main loop
