function trcpvzr(path,tracename,nfr,ifrlo,ifrhi,ifrndx,cgrmx,xbg,zbg,xfn,zfn)
%Construction of the grid of the spectra in 
%pressure, vertical and radial velocity
%for the given source linear trajectory 
%in layered medium (reciprocity used)
%function
%  trcpv
%   (path - string matrix,
%        first string - path to the directory
%        second  string - first two chars of *.wdx filename
%                          and first *.wld filename
%        third string - first two chars of next *.wld filename and so on
%   ,tracename - path to the file with results
%   ,nfr   - number of frequencies, power of 2 (256 or 1024)
%   ,ifrlo - lower nonzero spectrum sample number ( 0..nfr-1 )
%   ,ifrhi - higher nonzero spectrum sample number ( 0..nfr-1 )
%   ,ifrndx- frequency number with index file
%   ,cgrmx - estimation of the maximum group speed along the propagation path
%   ,xbg   - range coordinate of the initial point of the rectilinear trajectory
%   ,zbg   - depth coordinate of the initial point of the rectilinear trajectory
%   ,xfn   - range coordinate of the final point of the rectilinear trajectory
%   ,zfn   - depth coordinate of the final point of the rectilinear trajectory
%   )
%All *.wld files are assumed to have the identical structure
%--------------------------------------------------------
%  Results are stored in the MAT-file with the pathname [hflnm intrpl.mat] as
%  save [path tracename intrpl.mat] pr xr np rp , where
%  pr(ifrhi-ifrlo+1,nrc) - pressure frequency spectrum
%  vz(ifrhi-ifrlo+1,nrc) - vertical velocity frequency spectrum
%  xr - array of range coordinates of the grid points
%  np=[nar nfr ifrhi ifrlo nrc];
%  nar -
%  nrc - number of intrpolated points (receivers)
%  rp=[ztr xsrc fr ad crf cgrmx];
%  ztr - depth of the trajectory - assumed constant and equal to zbg
%  xsrc - range coordinate of the source
%  fr - frequency corresponding to the (nfr-1)th spectrum sample
%  ad - nfr-1
%  crf - reference sound speed on the propagation path

%-----------start
dbstop if error
disp(' ------- TrcPVzr running -------');
a=sprintf( ...
   'nfr=%d, ifrlo=%i, ifrhi=%d,\ncgrmx=%e, xbg=%f, zbg=%f' ...
   ,nfr,ifrlo,ifrhi,cgrmx,xbg,zbg);
disp(a);
disp(path);
I=sqrt(-1);
%---open ifrndx frequency *.wdx file
i=ifrndx;
nf=int2str(i);
nz=6-length(nf);
fz='';
for j=1:nz
   fz=[fz '0'];
end
hflnm=deblank(path(1,:))
clear a;
a=[hflnm deblank(path(2,:)) fz nf '.wdx'];
disp(['Index file is ' a]);
[fndx,msg]=fopen(a,'rt');
if fndx==-1,
   disp(msg);
   nrc=fndx;
   return;
end
%---get common parameters
wdx=fread(fndx);
fclose(fndx);
lwdx=size(wdx);
wdx=setstr(wdx');
b=sprintf('%s %i',a,lwdx);
disp(b);
%  [A,COUNT,ERRMSG,NEXTINDEX] = SSCANF(S,FORMAT,SIZE)
[nfl,cnt,errmsg,nxt]=sscanf(wdx,'%d',1);
wdx=wdx(1,nxt:lwdx); lwdx=lwdx-nxt+1;
[dz ,cnt,errmsg,nxt]=sscanf(wdx,'%g ',1);
wdx=wdx(1,nxt:lwdx); lwdx=lwdx-nxt+1;
[zsurf,cnt,errmsg,nxt]=sscanf(wdx,'%g ',1);
wdx=wdx(1,nxt:lwdx); lwdx=lwdx-nxt+1;
[xsrc,cnt,errmsg,nxt]=sscanf(wdx,'%g ',1);
wdx=wdx(1,nxt:lwdx); lwdx=lwdx-nxt+1;
c=real(wdx(1:70)); nxt=1; while c(nxt)~=10, nxt=nxt+1; end
wdx=wdx(1,nxt:lwdx); lwdx=lwdx-nxt+1;
[crf,cnt,errmsg,nxt]=sscanf(wdx,'%g ',1);
wdx=wdx(1,nxt:lwdx); lwdx=lwdx-nxt+1;
[fr,cnt,errmsg,nxt]=sscanf(wdx,'%g ',1);
wdx=wdx(1,nxt:lwdx); lwdx=lwdx-nxt+1;
[ad,cnt,errmsg,nxt]=sscanf(wdx,'%g ',1);
wdx=wdx(1,nxt:lwdx); lwdx=lwdx-nxt+1;
c=real(wdx(1:70)); nxt=1; while c(nxt)~=10, nxt=nxt+1; end
wdx=wdx(1,nxt:lwdx); lwdx=lwdx-nxt+1;
[c,cnt,errmsg,nxt]=sscanf(wdx,'%g ',[4,inf]);
nvx=cnt/4;
%-------------attention, please
xbg=c(1,1);
xfn=c(1,nvx);
%-------------------------------
xrl=xfn-xbg;
nrc=0;
xr=[];
zr=[];
ez=[];
flps=[];
xs=xsrc;
xp=[];
nis=[];
nfs=[];
fpss=[];
busy=1;
ivx=1;
%---scan *.wdx records
while 1,
   xp=xs;
   nip=nis;
   nfp=nfs;
   fpsp=fpss;
   if (ivx>nvx) break, end;
   x = c(1,ivx);
   ni= c(2,ivx);
   nf= c(3,ivx);
   fps=c(4,ivx);
   ivx=ivx+1;
   xs=x;
   nis=ni;
   nfs=nf;
   fpss=fps;
   if xbg==xfn, 
      %---this is vertical slice,only one output region assumed
      if sign(xp-xbg)~=sign(xs-xbg)&xp~=[],
         if abs(xp-xbg)>abs(xs-xbg) %---choose closer
            x =xs
            ni=nis;
            nf=nfs;
            fps=fpss
         else
            x=xp
            ni=nip;
            nf=nfp;
            fps=fpsp
         end
         nbg=max(ni,fix((zbg-zsurf)/dz+1.0))
         nfn=min(nf,fix((zfn-zsurf)/dz+1.0))
         for j=nbg:nfn
            nrc=nrc+1;
            ez(nrc)=1.0;
            xr(nrc)=xbg;
            zr(nrc)=dz*j-dz;
            flps(nrc)=(fps+j-ni-1)*8;
         end
         a='v finish'
         break,
      end
   else
      
      if x>=xbg & x<=xfn,
         zsi=zsurf+(ni-1)*dz;
         zsf=zsurf+(nf-1)*dz;
         ex=(x-xbg)/xrl;
         fx=(xfn-x)/xrl;
         zrf=fx*zbg+ex*zfn; %current depth
         if xs~=xp, busy=1; end
         disp(['x=' num2str(x) ', zrf=' num2str(zrf) ', busy=' int2str(busy)]);
         disp(['zsi=' num2str(zsi) ', zsf=' num2str(zsf) ', busy=' int2str(busy)]);
         if (zrf>zsi & zrf<=zsf) & busy,
            busy=0;
            nrc=nrc+1;
            ez(nrc)=rem((zrf-zsurf),dz);
            nfps=fix((zrf-zsurf)/dz+1.0)-ni;
         disp(['nrc=' num2str(nrc) ', nfps=' num2str(nfps) ', busy=' int2str(busy)]);
            flps(nrc)=(nfps+fps)*2-1;
            zr(nrc)=zrf;
            xr(nrc)=x;

         end
      end
   end
end % 1
a=sprintf('*.wdx file scan finished')
%------------- of forming control data accordingly *.wdx file
%
[nar,a]=size(path);
nar=nar-1;
nar=1; %this allows only one source-receiver pair
%----------read data and form the spectrum
for iar=1:nar
   pr=[];vz=[];apf=[];
   for ifr=ifrhi:-1:ifrlo
      nf=int2str(ifr);
      nz=6-length(nf);
      fz='';
      for j=1:nz
         fz=[fz '0'];
      end
      b=[hflnm deblank(path(iar+1,:)) fz nf '.wld']; disp(b);
      [ffld,msg]=fopen(b,'r');
      if ffld==-1,
         disp(msg)
         nrc=ffld;
         return;
      end
      fld=fread(ffld,inf,'float');
      fclose(ffld);
      ifrpos=ifr-ifrlo+1;
      ap=0.0;
      for irc=1:nrc,
         e=fld(flps(irc)  );
         f=fld(flps(irc)+1);
         if xbg==xfn,
            pl=e+I*f;
            pu=0.0+I*0.0;
         else
            g=fld(flps(irc)+2);
            h=fld(flps(irc)+3);
            pu=e+I*f;
            %pu=e-I*f;
            pl=g+I*h;
            %pl=g-I*h;
         end
         cyl=sqrt(xr(irc)-xsrc);
         p=((1.0-ez(irc))*pu+ez(irc)*pl);
         v=(pl-pu)/dz/(I*2*pi*fr/ad*ifr);
         if abs(xsrc-xr(irc))>0.001,
            p=p/cyl;v=v/cyl;
         end % if
%----full phase is added in pvzvxvy         
         pr(ifrpos,irc)=p;
         vz(ifrpos,irc)=v;
         ap=ap+abs(p);
      end % irc
      apf(ifr)=ap;
      disp(['ifr=' int2str(ifrpos) ', ap=' num2str(ap)]);
      vx(ifrpos,:)=([0 diff(pr(ifrpos,:))./diff(xr)]...
         +(I*2*pi*fr/ad*ifr)/crf*pr(ifrpos,:))/(I*2*pi*fr/ad*ifr);
   end % ifr
%------------start   
nrrf=20;
   figure;
   sgt=zeros(2*nfr,1);
   sgt(ifrlo:ifrhi)=exp(xr(nrrf)*(cgrmx-crf)/(crf*cgrmx)...
                       *I*2*pi*fr/ad*(ifrlo:ifrhi))'.*pr(:,nrrf);
   sgt=fft(sgt);
   ymax=max(real(sgt));
   subplot(3,3,1);
   plot(real(sgt));axis([1 size(sgt,1) -ymax ymax]);
   title(['range=' num2str(xr(nrrf))]);
   ylabel('p');
   subplot(3,3,2);
   plot(abs(sgt));axis([1 size(sgt,1) 0 ymax]);
   subplot(3,3,3);
   plot(unwrap(angle(sgt)));v=axis;axis([1 size(sgt,1) v(3) v(4)]);
%------vz   
   sgt=zeros(2*nfr,1);
   sgt(ifrlo:ifrhi)=exp(xr(nrrf)*(cgrmx-crf)/(crf*cgrmx)...
                       *I*2*pi*fr/ad*(ifrlo:ifrhi))'.*vz(:,nrrf);
   sgt=fft(sgt);
   subplot(3,3,4);
   plot(real(sgt));axis([1 size(sgt,1) -ymax ymax]);
   ylabel('vz');
   subplot(3,3,5);
   plot(abs(sgt));axis([1 size(sgt,1) 0 ymax]);
   subplot(3,3,6);
   plot(unwrap(angle(sgt)));v=axis;axis([1 size(sgt,1) v(3) v(4)]);
%------vx   
   sgt=zeros(2*nfr,1);
   sgt(ifrlo:ifrhi)=exp(xr(nrrf)*(cgrmx-crf)/(crf*cgrmx)...
                       *I*2*pi*fr/ad*(ifrlo:ifrhi))'.*vx(:,nrrf);
   sgt=fft(sgt);
   subplot(3,3,7);
   plot(real(sgt));axis([1 size(sgt,1) -ymax ymax]);
   ylabel('vx');
   subplot(3,3,8);
   plot(abs(sgt));axis([1 size(sgt,1) 0 ymax]);
   subplot(3,3,9);
   plot(unwrap(angle(sgt)));v=axis;axis([1 size(sgt,1) v(3) v(4)]);
   xlabel(['Source depth=100m']);
   zoom on;
%-----------------checking time gate at the middle of the trace
nrrf=fix(nrc/2);
   figure;
   sgt=zeros(2*nfr,1);
   sgt(ifrlo:ifrhi)=exp(xr(nrrf)*(cgrmx-crf)/(crf*cgrmx)...
                       *I*2*pi*fr/ad*(ifrlo:ifrhi))'.*pr(:,nrrf);
   sgt=fft(sgt);
   ymax=max(real(sgt));
   subplot(3,3,1);
   plot(real(sgt));axis([1 size(sgt,1) -ymax ymax]);
   title(['range=' num2str(xr(nrrf))]);
   ylabel('p');
   subplot(3,3,2);
   plot(abs(sgt));axis([1 size(sgt,1) 0 ymax]);
   subplot(3,3,3);
   plot(unwrap(angle(sgt)));v=axis;axis([1 size(sgt,1) v(3) v(4)]);
%------vz   
   sgt=zeros(2*nfr,1);
   sgt(ifrlo:ifrhi)=exp(xr(nrrf)*(cgrmx-crf)/(crf*cgrmx)...
                       *I*2*pi*fr/ad*(ifrlo:ifrhi))'.*vz(:,nrrf);
   sgt=fft(sgt);
   subplot(3,3,4);
   plot(real(sgt));axis([1 size(sgt,1) -ymax ymax]);
   ylabel('vz');
   subplot(3,3,5);
   plot(abs(sgt));axis([1 size(sgt,1) 0 ymax]);
   subplot(3,3,6);
   plot(unwrap(angle(sgt)));v=axis;axis([1 size(sgt,1) v(3) v(4)]);
%------vx   
   sgt=zeros(2*nfr,1);
   sgt(ifrlo:ifrhi)=exp(xr(nrrf)*(cgrmx-crf)/(crf*cgrmx)...
                       *I*2*pi*fr/ad*(ifrlo:ifrhi))'.*vx(:,nrrf);
   sgt=fft(sgt);
   subplot(3,3,7);
   plot(real(sgt));axis([1 size(sgt,1) -ymax ymax]);
   ylabel('vx');
   subplot(3,3,8);
   plot(abs(sgt));axis([1 size(sgt,1) 0 ymax]);
   subplot(3,3,9);
   plot(unwrap(angle(sgt)));v=axis;axis([1 size(sgt,1) v(3) v(4)]);
   xlabel(['Source depth=100m']);
   zoom on;
%-----------------checking time gate at the end of the trace
nrrf=nrc;
   figure;
   sgt=zeros(2*nfr,1);
   sgt(ifrlo:ifrhi)=exp(xr(nrrf)*(cgrmx-crf)/(crf*cgrmx)...
                       *I*2*pi*fr/ad*(ifrlo:ifrhi))'.*pr(:,nrrf);
   sgt=fft(sgt);
   ymax=max(real(sgt));
   subplot(3,3,1);
   plot(real(sgt));axis([1 size(sgt,1) -ymax ymax]);
   title(['range=' num2str(xr(nrrf))]);
   ylabel('p');
   subplot(3,3,2);
   plot(abs(sgt));axis([1 size(sgt,1) 0 ymax]);
   subplot(3,3,3);
   plot(unwrap(angle(sgt)));v=axis;axis([1 size(sgt,1) v(3) v(4)]);
%------vz   
   sgt=zeros(2*nfr,1);
   sgt(ifrlo:ifrhi)=exp(xr(nrrf)*(cgrmx-crf)/(crf*cgrmx)...
                       *I*2*pi*fr/ad*(ifrlo:ifrhi))'.*vz(:,nrrf);
   sgt=fft(sgt);
   subplot(3,3,4);
   plot(real(sgt));axis([1 size(sgt,1) -ymax ymax]);
   ylabel('vz');
   subplot(3,3,5);
   plot(abs(sgt));axis([1 size(sgt,1) 0 ymax]);
   subplot(3,3,6);
   plot(unwrap(angle(sgt)));v=axis;axis([1 size(sgt,1) v(3) v(4)]);
%------vx   
   sgt=zeros(2*nfr,1);
   sgt(ifrlo:ifrhi)=exp(xr(nrrf)*(cgrmx-crf)/(crf*cgrmx)...
                       *I*2*pi*fr/ad*(ifrlo:ifrhi))'.*vx(:,nrrf);
   sgt=fft(sgt);
   subplot(3,3,7);
   plot(real(sgt));axis([1 size(sgt,1) -ymax ymax]);
   ylabel('vx');
   subplot(3,3,8);
   plot(abs(sgt));axis([1 size(sgt,1) 0 ymax]);
   subplot(3,3,9);
   plot(unwrap(angle(sgt)));v=axis;axis([1 size(sgt,1) v(3) v(4)]);
   xlabel(['Source depth=100m']);
   zoom on;
   ztr=zr(1);
   np=[nar nfr ifrhi ifrlo nrc];
   rp=[ztr xsrc fr ad crf cgrmx];
   save([hflnm tracename 'trcfir1.mat']...
      ,'np','rp','xr','pr','vz','vx');
   figure;
   plot(apf);
end; %iar
disp(' -------- TrcPVzr finished -------- ');
