clear all
close all
fclose all
%---------------------------Input values
hflnm='d:\1cva\bering\br';
nfr=1024; %-1st frequency sample is 0, 257th is 0, total 256 complex numbers
ifrhi=2047;% these are C-type indices (0..2047)
ifrlo=1792;
cgrmx=1.550%1.484; %-close estimation
xbg=30.65;
zbg=0.0;
xfn=30.65;
zfn=0.042;
%-----------start
I=sqrt(-1);
%---open highest frequency *.wdx file
i=ifrhi;
nf=int2str(i);
nz=6-length(nf);
fz='';
for j=1:nz
   fz=[fz '0'];
end
[fnm,pthnm]=uigetfile('br*.wdx','Choose the index file');
if fnm==0, return, end;
%-----------------------------
hflnm=[pthnm 'br']
fnmndx=[pthnm fnm];
[fndx,msg]=fopen(fnmndx,'rt');
if fndx==-1,
   a,
   msg,
   nrc=fndx;
   return;
end
%---get common parameters
nfl=fscanf(fndx,'%d ',1);
dz =fscanf(fndx,'%g',1);
zsurf=fscanf(fndx,'%g ',1);
xsrc=fscanf(fndx,'%g',1);
%--------------
a=fgetl(fndx)
crf=fscanf(fndx,'%g',1);
fr=fscanf(fndx,'%g',1);
ad=fscanf(fndx,'%g',1);
a=fgetl(fndx);
xrl=xfn-xbg;
nrc=0;
xr=[];
zr=[];
ez=[];
flps=[];
xs=xsrc;
xp=[];
nis=[];
nfs=[];
fpss=[];
busy=1;
%---scan *.wdx records
while 1,
   xp=xs;
   nip=nis;
   nfp=nfs;
   fpsp=fpss;
   x = fscanf(fndx,'%g ',1);
   if (feof(fndx) == 1) break, end;
   ni=fscanf(fndx,'%d ',1);
   nf=fscanf(fndx,'%d ',1);
   fps=fscanf(fndx,'%ld',1);
   xs=x;
   nis=ni;
   nfs=nf;
   fpss=fps;
   if xbg==xfn, %---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 sign(x-xbg)~=sign(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;
         if xs~=xp, busy=1; end
         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;
            flps(nrc)=(nfps+fps-1)*8;
            zr(nrc)=zrf;
            xr(nrc)=x;
         end
      end
   end
end % 1
fclose(fndx);
a=sprintf('*.wdx file scan finished')
%------------------------- of forming control data accordingly *.wdx file
%
%----------read data and form the spectrum
pr=[];
for ifr=ifrhi:-1:ifrlo
   nf=int2str(ifr);
   nz=6-length(nf);
   fz='';
   for j=1:nz
      fz=[fz '0'];
   end
   b=[hflnm fz nf '.wld'];
   [ffld,msg]=fopen(b,'r');
   if ffld==-1,
      nrc=fndx;
      return;
   end
   ifrpos=ifr-ifrlo+1;
   for irc=1:nrc,
      ist=fseek(ffld,flps(irc),-1);
      e=fread(ffld,1,'float');
      f=fread(ffld,1,'float');
      if ist~=0,
         irc,flps(nrc)
         e,f
      end
      if xbg==xfn,
         pl=e+I*f;
         pu=0.0+I*0.0;
      else
         g=fread(ffld,1,'float');
         h=fread(ffld,1,'float');
         pu=e+I*f;
         pl=g+I*h;
      end
      phase=ifr/ad*2*pi*fr*(1.0/crf-1.0/cgrmx)*(xr(irc)-xsrc);
      cyl=sqrt(xr(irc)-xsrc);
      %      p=exp(I*phase)*((1.0-ez(irc))*pu+ez(irc)*pl);
      p=((1.0-ez(irc))*pu+ez(irc)*pl);
      p=exp(I*phase)*p;
      if abs(xsrc-xr(irc))>0.001,
         p=p/cyl;
      end % if
      pr(ifrpos,irc)=p;
   end % irc
   fclose(ffld);
end % ifr
%end of arrir
a=sprintf('*.wld file scan finished')
ir=zeros(2*(ifrhi+1),nfn);
for i=1:nfn,
   for k=ifrlo:ifrhi,
      ir(k,i)=pr(k-ifrlo+1,i);
   end;
   ir(:,i)=fft(ir(:,i));%Re and Im of ir are quadrature components
end;
a=['save ' [hflnm '2047.mat'] ' pr ir xr zr nfn nfr ifrhi ifrlo fr ad cgrmx fnmndx']
eval(a);
