program geometry
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c  compute receiver coordinates along curved line 
c
c  input file should contain 3 columns with:
c     shot_point_number  source_UTM_east  source_UTM_north
c  
c  output file will contain 8 columns with:
c     shot_point_number receiver_group_number source_UTM_east \
c     source_UTM_north receiver_UTM_east receiver_UTM_north \
c     CMP_number offset
c
c  an additional file named 'cmp.dat' will be generated consisting of
c     7 columns: CMP_number shot_point_number receiver_group_number \
c                CMP_UTM_east CMP_UTM_north CMP_inline CMP_crossline
c     (all coordinates in meters)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      implicit none
      integer maxnsp,i,j,isp,lsp,nsp,ng,cmp,icmp
      parameter (maxnsp=10000)
      double precision sp(maxnsp),xc(maxnsp),yc(maxnsp),dist(maxnsp),
     +                 splx(maxnsp),sply(maxnsp),w1(maxnsp),w2(maxnsp),
     +                 x1,x2,y1,y2,ss,sg,ds,dg,dd,d1,d2,xs,ys,xr,yr,
     +                 spl3,dx,dy,x0,y0,xm,ym,dcmp,xil,yxl,offset
      character*80 infile,outfile

      write(*,*)'Name of file with UTM coordinates of shot points:'
      read(*,*)infile

      open(10,file=infile,status='old')
      lsp=0
      do 10 i=1,maxnsp
         read(10,*,end=11)isp,xc(i),yc(i)
         sp(i)=isp
         if (lsp.eq.0) write(*,*)'First shot is',isp
         if (lsp.eq.0) lsp=isp-1
         if (isp.eq.lsp) write(*,*)'Duplicate shot',isp
         if (isp.lt.lsp) write(*,*)'Duplicate shots',lsp,'-',isp
         if ((isp-lsp).eq.2) write(*,*)'Missing shot',lsp+1
         if ((isp-lsp).gt.2) write(*,*)'Missing shots',lsp+1,'-',isp-1
         lsp=isp
 10   continue
 11   nsp=i-1
      write(*,*)'Last shot is',lsp
      close(10)

      write(*,*)'Distance between antenna and source (m):'
      read(*,*) d1
c     d1=67.
      write(*,*)'Distance between source and first hydrophone group (m):
     +'
      read(*,*) d2
      write(*,*)'Distance between hydrophone groups (m):'
      read(*,*) dg
      write(*,*)'Number of hydrophone groups:'
      read(*,*) ng
      write(*,*)'Distance between CMPs (m):'
      write(*,*)'(Typically half the distance between hydrophone groups)
     +'
      read(*,*) dcmp
c     dcmp=dg/2
      write(*,*)'CMP number below first shot point:'
      read(*,*)icmp
      write(*,*)'Name of output file:'
      read(*,*)outfile

c  compute spline coefficients for x and y coordinates
c  as function of shot number

      call zspl3(nsp,sp,xc,w1,w2,splx)
      call zspl3(nsp,sp,yc,w1,w2,sply)

c  compute distance along curved line by dividing distance between
c  shots in 16, interpolate coordinates and sum line segments

      dist(1)=0.
      do 20 i=1,nsp-1
         dd=0.
         x1=xc(i)
         y1=yc(i)
         ss=sp(i)
         ds=(sp(i+1)-sp(i))/16.
         do 22 j=1,16
            ss=ss+ds
            x2=spl3(nsp,sp,xc,splx,ss)
            y2=spl3(nsp,sp,yc,sply,ss)
            dd=dd+dsqrt((x2-x1)**2+(y2-y1)**2)
            x1=x2
            y1=y2
 22      continue
         dist(i+1)=dist(i)+dd
 20   continue

c  compute spline coefficients for x and y coordinates
c  as function of distance along curved line

      call zspl3(nsp,dist,xc,w1,w2,splx)
      call zspl3(nsp,dist,yc,w1,w2,sply)

c  compute coefficients of inline-crossline coordinate system

      dx=xc(nsp)-xc(1)
      dy=yc(nsp)-yc(1)
      dd=dsqrt(dx*dx+dy*dy)
      dx=dx/dd
      dy=dy/dd
      x0=xc(1)
      y0=yc(1)

c  open output file

      open(11,file=outfile,status='unknown')
      open(12,file='cmp.dat',status='unknown')
      do 30 i=1,nsp
         isp=sp(i)
         ss=dist(i)-d1
         if (ss.gt.0) then
            xs=spl3(nsp,dist,xc,splx,ss)
            ys=spl3(nsp,dist,yc,sply,ss)
         else
            dd=dsqrt((xc(2)-xc(1))**2+(yc(2)-yc(1))**2)
            xs=xc(1)+ss*(xc(2)-xc(1))/dd
            ys=yc(1)+ss*(yc(2)-yc(1))/dd
         endif
         do 35 j=1,ng
            sg=ss-d2-(j-1)*dg
            if (sg.gt.0) then
               xr=spl3(nsp,dist,xc,splx,sg)
               yr=spl3(nsp,dist,yc,sply,sg)
            else
               dd=dsqrt((xc(2)-xc(1))**2+(yc(2)-yc(1))**2)
               xr=xc(1)+sg*(xc(2)-xc(1))/dd
               yr=yc(1)+sg*(yc(2)-yc(1))/dd
            endif
            xm=(xs+xr)/2
            ym=(ys+yr)/2
            xil=dx*(xm-x0)+dy*(ym-y0)
            yxl=-dy*(xm-x0)+dx*(ym-y0)
            cmp=icmp+nint(xil/dcmp)
            offset=dsqrt((xs-xr)**2+(ys-yr)**2)
            write(11,'(i5,i4,4i8,i6,i5)')isp,j,nint(xs),nint(ys),
     +         nint(xr),nint(yr),cmp,nint(offset)
            write(12,'(2i6,i4,2i8,i7,i5)')cmp,isp,j,nint(xm),nint(ym),
     +         nint(xil),nint(yxl)
 35      continue
 30   continue
      close(11)
      close(12)

      end