      program qscan
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c                      Q S C A N
c
c  Reads Riegl Q240 new scanner data and output ground points
c  input file must contain GPS coordinates and pitch, roll and hdg
c  Format of coor/pitch file as ".pos"-output from "gpsins" or "readegi".
c  Time assumed to be UTC in all files 
c  This version assumes GPS synchronization, no GPS synchronization
c  not tested. Built on readscan7
c
c  input: 
c
c  scanfile 
c  coorfile
c  geoidfile    (0: no geoid file)
c  outfile1     (0: no scanfile, only central point in scan)
c  outfile2
c  thr,min,sec  (0,0,0 ok for GPS synchronization)
c  dtsec, lgps_sync,
c  t1, t2 
c  lat1 lat2 lon1 lon2
c  nave, ithin, iout_type, lpr, lcalib, lgeo
c  dxl, dyl, dzl, pitch0, roll0, hdg0
c
c  where:
c
c  thr min sec  time of first scan in file (thr can be dec hr)
c               0 0 0 means GPS synchronization
c  dtsec        error in sec of start time or GPS sync time
c  lgps_sync    t=use gps sync, f=force use of thr,min,sec input 
c  t1 t2        wanted start and stop time of output data (dechr) 
c  lat1 lat2 lon1 lon2    gps position window (for aircraft),
c               data are only analyzed over this area (degrees)
c  nave         no of points averaged across scan (1 = no averaging)
c  ithin        only use every ithin scan (1 = use all)
c  iout_type    1=std text, 2=binary, 3=paraglider format  
c  lpr          t=use pitch/roll information in coorfile
c  lcalib       t=make regression for roll0, ocean or lake flight
c               (NB: the estimated roll0 value is not applied) 
c  lgeo         Normal geographic output, false for track time in hr
c  dxl dyl dzl  offset vector from GPS to scanner in body system:
c               x forward, y right, z down
c               (twin-otter back antenna:  0 0 1.66)
c  pitch0 roll0 hdg0    misorientation of laser scanner (degrees)
c
c  outfile1: thr,lat,lon,ground_h,amplitude,iscan,n_ave for all pts
c  outfile2: do + flightelev,range, 
c            for vertical (possibly averaged) scan only
c  
c  (c) Rene Forsberg, KMS/DNSC/NSI, Denmark, nov 2000
c  changed to include new laser and exact formulas, 24 jun 2001
c  updated with roll0 regression, rf 2 aug 2001
c  updated with time output
c  updated cfs alert may2002 with time along/cross par + geoid
c  updated jan 2003 with traps for timer errors (see also rieg2txt)
c  updates Alert May 2004 (time output)
c  binary output update, Moscow airport, Oct 2004
c  (version can now be compiled with lahey or ming - lahey is faster)
c  changed to Q240 and lahey f95 compiler (only)
c  Updates for time sync, averaging and iscandir, oct 2008, rf 
c  NB: iscandir = 1 for DTU-Space 2008 (-1 for paraglider iout_type=3)
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      implicit double precision(a-h,o-z)
      character*36 ifile,pfile,gfile,ofile1,ofile2
      logical lbino,ltxt,lpr,lcalib,lgeo,lgeoid,lgps,lgps_sync,
     .lscan,lspike
      dimension gpst(48801),rlat(48801),rlon(48801),rh(48801),
     .pitch(48801),roll(48801),hdg(48801)
      real*4 ggrid(250000),rvel(48801)
c      
      integer*2 syncseq,bytespl,measid,linecount,i2
      integer*1 bampl,polarid,i1,iampo,i1o,iflags
      integer*1 brange(4),bang(4),btime(4),bscount(4),bstimer(4)
      integer*4 irange,iang,itime,iscount,istimer
      equivalence(brange(1),irange),(bang(1),iang),(btime(1),itime),
     .(bscount(1),iscount),(bstimer(1),istimer),(bampl,iampl),
     .(linecount,ilinecount) 
      dimension range(1000),timer(1000),vgon(1000),iamp(1000)
c
      nmax = 48801
      igdim = 250000
      iscandim = 1000
c
c  reset equivalence integers 
c
      irange = 0
      iang = 0
      itime = 0
      iscount = 0
      istimer = 0   
c
      radeg = 180.d0/3.1415926536d0
      radgon = radeg/90*100
      degm = 111120
c

      read(*,'(a)') ifile
      read(*,'(a)') pfile
      read(*,'(a)') gfile
      read(*,'(a)') ofile1
      read(*,'(a)') ofile2
c
      read(*,*) thr,min,sec
      read(*,*) dtsec, lgps_sync
      read(*,*) t1,t2
      read(*,*) rlat1,rlat2,rlon1,rlon2
      read(*,*) nave,ithin,iout_type,lpr,lcalib,lgeo
      read(*,*) dxl,dyl,dzl
      read(*,*) pitch0,roll0,hdg0
c
c  open files - special byte handling, can only compile in lahey fortran
c
      open(10,file=ifile,form='binary',
     .access='transparent',status='old')
c
      lbino = (iout_type.eq.2)
      ltxt = (iout_type.eq.1)
      if (iout_type.gt.3) stop '*** iout_type wrong, check input'
      lscan = (ofile1.ne.'0')   
      if (.not.lscan) lcalib = .false.
c 
      open(12,file=pfile,status='old')
      if (lscan) then
        if (lbino) then
          open(20,file=ofile1,form='binary',access='transparent',
     .    status='unknown')
        else          
          open(20,file=ofile1,status='unknown')
        endif
      endif
      open(21,file=ofile2,status='unknown')
c  
      if (lcalib) open(22,file='scan.cal',status='unknown')
c
      write(*,4) 
4     format(' === QSCAN === output Q240 scanner surface points')
      if (.not.lscan) write(*,*) '*** only vertical output'
      if (ithin.gt.1) 
     .write(*,*) '- Thinned output, every:',ithin,' scan'
      write(*,5) dxl,dyl,dzl,pitch0,roll0,hdg0
5     format(' offset/misalignment parameters used: ',6f6.2)
c
c  read geoid file
c  ---------------
c
      lgeoid = (gfile.ne.'0')
      if (lgeoid) then
        open(14,file=gfile,status='old')
        read(14,*) rfi1,rfi2,rla1,rla2,dfi,dla
        nng = (rfi2-rfi1)/dfi + 1.5
        neg = (rla2-rla1)/dla + 1.5
        if (nng*neg.gt.igdim) then
          write(*,*) nng*neg,igdim
          stop '*** geoid grid too large, sorry'
        endif
        do 7 i = nng,1,-1
          j = (i-1)*neg
          read(14,*) (ggrid(j+k),k=1,neg)
7       continue
        write(*,*) '- geoid file read OK -',nng,neg
        close(14)
        endif
c
c  read file with gps and pitch/roll data
c  first number in file is UTC time in dechr 
c  -----------------------------------------
c
      i = 1
10    if (lpr) then
        read(12,*,end=12)
     .  gpst(i),rlat(i),rlon(i),rh(i),pitch(i),roll(i),hdg(i)
        if (hdg(i).gt.180) hdg(i) = hdg(i)-360
      else
        read(12,*,end=12) gpst(i),rlat(i),rlon(i),rh(i)
        pitch(i) = 0
        roll(i) = 0
        hdg(i) = 0
      endif
      igps = 0
c
      if (gpst(i).lt.t1-1.d0/3600) goto 10
      if (gpst(i).gt.t2+1.d0/3600) goto 12
      i = i+1
      if (i.ge.nmax) stop 'too many gps points - increase nmax'
      goto 10
12    ngps = i-1
      dtmin = 9999.9
      dtmax = 0
      rvelmin = 9999.9
      rvelmax = -9999.9
      do i = 2, ngps
        dt = (gpst(i)-gpst(i-1))*3600
        dr = sqrt((rlat(i)-rlat(i-1))**2 + 
     .  ((rlon(i)-rlon(i-1))*cos(rlat(i)/radeg))**2)*111194.d0
        rvel(i) = dr/dt*3.6
        if (dt.gt.dtmax) dtmax = dt
        if (dt.lt.dtmin) dtmin = dt
        if (rvel(i).gt.rvelmax) rvelmax = rvel(i)
        if (rvel(i).lt.rvelmin) rvelmin = rvel(i)
      enddo
      rvel(1) = rvel(2)
      write(*,13) ngps,dtmin,dtmax,rvelmin,rvelmax
13    format(/' GPS read:',i6,
     .' gps/attitude data, min max spacing (sec):',2f9.2,
     ./' minimum and maximum velocity (km/h): ',2f8.1)
c       
c  main laserscanner reading code
c  ------------------------------
c
c  the reading is restricted to Q240, for old scanner see readscan
c  the Q240 scan are read scan by scan because timing is in trailer record
c
c  scanner installation orientation:
c  iscandir = 1 for scan from left to right, -1 right to left 
c
c  NSI2008: iscandir = 1
c  Ling Vinghof paraglider: iscandir = -1
C
c  dx dy dz is offset vector from GPS antenna to laser in acft body system
c  x forward  y right  z down
c
      if (iout_type.le.2) then
        iscandir = 1
      else
        iscandir = -1
      endif
c
      timer0 = 0
      igpssec = 0
c
      write(*,*)
      write(*,*) 'Reading scanner data from: ',ifile
c
c  loop initialization parameters
c  iscan is counting scan no
c
      throp = 0
      angp = 999
      nlas = 0
c
      rmin = 9999
      rmax = -9999
      ahmin = 9999
      ahmax = -9999
      rfimin = 9999
      rfimax = -9999
      rlamin = 9999
      rlamax = -9999
      isum = 0
      sumal = 0
      sumal2 = 0
      sumy = 0
      sumaly = 0
c
      iscan = 1
      thrp = 0
      thro = 0
c
c  read riegl header
c  thr0 is start time
c  ------------------
c
      call readhdr(npscan,rangefak,anglefak,timerfak,measid,
     .  isync,polarid,thr0)
      if (measid.ne.77) stop '*** only measis 77 implemented'
      if (npscan.gt.iscandim) stop '*** increase iscandim'
      nmod = 400.d0/anglefak/polarid
c
      if (thr0.eq.0) then
        thr0 = thr+min/60.d0+sec/3600
        lgps = .false.
        write(*,*) '- Riegl header record: no gps sync found'
      else
        if (lgps_sync) then
          lgps = .true.
        else
          write(*,*) '- GPS sync found, but not used'
          lgps = .false.
        endif  
      endif
c
      ihr = thr0
      imin = (thr0-ihr)*60
      write(*,11) thr0,ihr,imin,(thr0-ihr-imin/60.d0)*3600,
     .dtsec
11    format(' file start epoch:',f9.5,' ',i5,i3,f6.2,', dtsec:',f7.2)
c
c  scan record loop - store in array to get timing from trailer
c  entry point for loop over scans
c  ------------------------------------------------------------
c
20    read(10,err=90) syncseq
      if (syncseq.ne.isync) 
     .write(*,*) '*** warning - syncseq not ok',syncseq,isync
      do 30 i = 1,npscan
        read(10,err=90) brange(1),brange(2),brange(3),bampl,
     .  bang(1),bang(2),bang(3),btime(1),btime(2),btime(3)
c
        vgon(i) = 2*mod(iang,nmod)*anglefak
        range(i) = irange*rangefak
        timer(i) = itime*timerfak
        iamp(i) = iampl
30    continue
c
c  read scan trailer and complete timer
c  linecount er riegls "dataportlinecounter" 
c  2 bytes = 65535 max
c  -----------------------------------------
c
      read(10) i1,linecount,iflags,bscount(1),bscount(2),bscount(3),
     .bstimer(1),bstimer(2),bstimer(3)
      if (iscan.eq.1) then 
        write(*,31) iscount + istimer*timerfak 
31      format(/' First scan epoch (seconds from GPS init):',f12.6)
        istimerp = istimer
        lc0 = ilinecount
      elseif (.not.lgps) then
c
c  check for timer overruns - see trailer structure in Riegl manual
c
        if (istimer-istimerp.le.-16700000) 
     .  istimer = istimer+16777215
      endif
c
      ilinecount = ilinecount-lc0+1
      if (ilinecount.lt.ilinecountp.and.iscan.gt.1) 
     .ilinecount = ilinecount+65336
      ilinecountp = ilinecount
c
c  trap if counter increments but scount is still high (riegl error?)
c
      if (iscount.gt.iscountp.and.iscan.gt.1) then
        if (istimer*timerfak.gt.0.999) iscount = iscount-1
      endif
      iscountp = iscount
c
      thrscan = (iscount+istimer*timerfak + dtsec)/3600+thr0
c
      if (iscan.le.3) then
        write(*,32) iscan,thrscan,iflags
32      format(' scan no:',i4,', time (hr)', f10.6,
     .  ', GPS sync flags: ',b7.7)
        if (iscan.eq.3) then
          write(*,321) 
321       format(' Centre of scan data every 30 sec: ',
     .    /' scan no   thr       lat        lon        h   ampl',
     .    '  flags')
        endif
      endif
c
c  skip scan if thinning mode or time before start or beyond end
c  --------------------------------------------------------------     
c
      if (thrscan.lt.t1) then
        if (mod(iscan,10000).eq.0) write(*,33) iscan,thrscan
33      format(i7,f9.5,' skipped, before t1')
        thrp = thrscan
        goto 49
      endif
      if (thrscan.gt.t2) goto 90
      if (ithin.gt.1.and.mod(iscan,ithin).ne.1) goto 49
c
c  loop over points in scan 
c  correct times with dtsec
c  ------------------------
c
      thrsum = 0
      glatsum = 0
      glonsum = 0
      ahsum = 0
      ghsum = 0
      angsum = 0
      ampsum = 0
      nsum = 0
c
      do 40 i = 1,npscan
        thr = (timer(i)+iscount+istimer*timerfak + dtsec)/3600 + thr0
        if (i.eq.npscan.and.iscan.eq.1) thr1 = thr
c
        if (range(i).lt.1.5) goto 40
c
c  interpolate gps pos at time thr
c  -------------------------------
c  find which index igps,igps+1 to use in GPS file
c
42      if (gpst(igps+1).lt.thr) then 
          igps = igps+1
          if (igps.gt.ngps-1) then 
            igps = ngps-1
          else
            goto 42
          endif
        endif  
        r = (thr-gpst(igps))/(gpst(igps+1)-gpst(igps))
        if (r.lt.0) r = 0
        if (r.gt.1) r = 1
c
        alat = (1-r)*rlat(igps) + r*rlat(igps+1)
        alon = (1-r)*rlon(igps) + r*rlon(igps+1)
        ah = (1-r)*rh(igps) + r*rh(igps+1)
c
c  only use data inside geographical limits
c
        if (alat.lt.rlat1.or.alat.gt.rlat2) goto 40
        if (alon.lt.rlon1.or.alon.gt.rlon2) goto 40
c
        if (ah.lt.ahmin) ahmin = ah
        if (ah.gt.ahmax) ahmax = ah
        if (range(i).lt.rmin) rmin = range(i)
        if (range(i).gt.rmax) rmax = range(i) 
        if (alat.lt.rfimin) rfimin = alat 
        if (alat.gt.rfimax) rfimax = alat 
        if (alon.lt.rlamin) rlamin = alon 
        if (alon.gt.rlamax) rlamax = alon
c
        aroll = (1-r)*roll(igps) + r*roll(igps+1) + roll0
        apitch = (1-r)*pitch(igps) + r*pitch(igps+1) + pitch0
c
c  check heading for interpolation at 180
c
        hdg1 = hdg(igps)
        hdg2 = hdg(igps+1)
        if (hdg1.gt.90.and.hdg2.lt.-90) hdg2 = hdg2+360
        if (hdg1.lt.-90.and.hdg2.gt.90) hdg1 = hdg1+360
        ahdg = (1-r)*hdg1 + r*hdg2 + hdg0
c
c  transfer coordinate with interpolated position and attitude values
c  see honeywell manual 6.5.3
c  local geodetic: x north, y west, z up
c
        ang = iscandir*(vgon(i)-100)/radgon
c
        if (lgeo) then
          cosa = cos(ahdg/radeg)
          sina = sin(ahdg/radeg)
          cosp = cos(apitch/radeg)
          sinp = sin(apitch/radeg)
          cosr = cos(aroll/radeg)
          sinr = sin(aroll/radeg)
          coslat = cos(alat/radeg)
c
          dy = -sin(ang)*range(i) + dyl
          dz = cos(ang)*range(i) + dzl
          dxg = cosa*cosp*dxl + (cosa*sinp*sinr-sina*cosr)*dy
     .                                + (cosa*sinp*cosr+sina*sinr)*dz
          dyg = -sina*cosp*dxl - (sina*sinp*sinr+cosa*cosr)*dy
     .                                + (cosa*sinr-sina*sinp*cosr)*dz
          dzg = sinp*dxl - cosp*sinr*dy - cosp*cosr*dz
c
          glat = alat + dxg/degm
          glon = alon - dyg/(degm*coslat)
          gh = ah + dzg 
          if (lgeoid) gh = gh - 
     .    bilin((glat-rfi1)/dfi+1,(glon-rla1)/dla+1,ggrid,nng,neg,igdim)
        else
c
c  approximative along-track/across-track time in hr
c  add cross-track bias for every 0.02 hr = 1.2 min
c
          vel = ((1-r)*rvel(igps) + r*rvel(igps+1))*1000
          glat = thr + sin(apitch/radeg)*range(i)/vel
          glon = - sin(ang+aroll/radeg)*range(i)/vel
          gh = ah - cos(apitch/radeg)*cos(ang+aroll/radeg)*range(i) - dzl
          if (lgeoid) gh = gh - 
     .    bilin((alat-rfi1)/dfi+1,(alon-rla1)/dla+1,ggrid,nng,neg,igdim)
        endif
c
c  average data if nave > 1
c  ------------------------
c
        thro = thr
        if (nave.gt.1) then
c
c  same scan, accumulate values within nave
c
          if (mod(i,nave).ne.0) then
            thrsum = thrsum + thr
            glatsum = glatsum + glat
            glonsum = glonsum + glon
            ahsum = ahsum + ah
            ghsum = ghsum + gh
            angsum = angsum + ang
            ampsum = ampsum + iamp(i)
            nsum = nsum+1
            goto 40
          else
c
c  write sum of previous values and start accumulate new value    
c
            if (nsum.ge.1) then
              thro = thrsum/nsum
              glato = glatsum/nsum
              glono = glonsum/nsum
              gho = ghsum/nsum
              aho = ahsum/nsum
              ango = angsum/nsum
              iampo = ampsum/nsum+0.5
              nsumo = nsum
c
              thrsum = thr
              glatsum = glat
              glonsum = glon
              ghsum = gh
              ahsum = ah
              angsum = ang
              ampsum = iamp(i) 
              nsum = 1
            else
              thrsum = 0
              glatsum = 0
              glonsum = 0
              ahsum = 0
              ghsum = 0
              angsum = 0
              ampsum = 0
              nsum = 0
              goto 40
            endif
c
            glat = glato
            glon = glono
            gh = gho
            ah = aho
            ang = ango
            iamp(i) = iampo
          endif
        endif
c           
c  output of results
c  -----------------
c
        iampo = iamp(i)
c
c  swath data - binary output do not contain nsumo
c  binary record 18 bytes time(4),lat(4),lon(4),h(4),amp(1),i(1)
c
        if (lscan) then
          if (lbino) then
            ithro = nint(thro*10000000)
            iglat = nint(glat*10000000)
            iglon = nint(glon*10000000)
            igh = nint(gh*1000)
            i1o = i
            write(20) ithro,iglat,iglon,igh,iampo,i1o
          elseif (ltxt) then
            if (nave.le.1) then
              write(20,44) thro,glat,glon,gh,iampo,i
            else
              write(20,44) thro,glat,glon,gh,iampo,i,nsumo
            endif
          else                             
            if (nave.le.1) then
              write(20,44) thro,glat,glon,gh,iampo,i
            else
              write(20,44) thro,glat,glon,gh,iampo,i,nsumo
            endif
          endif
        endif
44      format(f9.6,f10.6,f11.6,f8.2,1x,i3,i4,i3)
c
c  vertical data
c
        if (angp.lt.0.and.ang.ge.0) then
          if (iscan.eq.1.or.thro-throp.gt.0.0083333) then
            write(*,45) iscan,thro,glat,glon,gh,iampo,iflags
            throp = thro
          endif
          write(21,46) thro,glat,glon,gh,iampo,ah,range(i)
        endif
45      format(' ',i6,f11.7,f10.6,f11.6,f8.2,1x,i3,1x,b7.7)
46      format(' ',f11.7,f10.6,f11.6,f8.2,1x,i3,2f8.2)
        nlas = nlas+1
        angp = ang
c
c  linear regression parameters gh versus cross-offset
c
        if (lcalib) then
          alfa = tan(ang)*ah
          isum = isum+1
          sumal = sumal + alfa
          sumal2 = sumal2 + alfa**2
          sumy = sumy + gh
          sumaly = sumaly + alfa*gh
          write(22,47) alfa,gh
47       format(2f10.2)
        endif
c
c  end of scan loop, spike check
c  -----------------------------
c
        if (abs(thr-thrp).gt.0.0001.and.(i.gt.1.or.iscan.gt.1)) 
     .  write(*,48) iscan,i,thr,abs(thr-thrp)*3600
48      format(' ** time gap scan',i7,i4,f11.7,', seconds',f10.3)
        thrp = thr
40    continue
c
49    iscan = iscan+1
      goto 20
c
c  output information on scans
c  ----------------------------
c
90    write(*,*)
      write(*,*) 'Total number of scans read: ', iscan
      write(*,*) 'Total number of points output: ',nlas
      write(*,91) rmin,rmax
91    format(' min and max laser range measured: ',2f9.2)
      write(*,92) ahmin,ahmax
92    format(' min and max flight elevation: ',2f9.2)
      write(*,93) rfimin,rfimax,rlamin,rlamax
93    format(' aircraft positions within area: ',4f9.4)
      write(*,*)
      write(*,*) 
     .'scanpt data (thr,lat,lon,h,ampl,i) output on: ',ofile1
      write(*,*) 
     .'center pts (thr,lat,lon,h,ampl,i,airh,range): ',ofile2
c
      if (lcalib) then
        write(*,*)
        write(*,*) 'Regression of surface relative to ellipsoid/geoid:'
        write(*,*) 'no of points: ',isum
        r = isum*sumal2 - sumal**2
        if (r.eq.0) r = 1.d-99
        bias = (sumal2*sumy - sumal*sumaly)/r
        alfa = (-sumal*sumy + isum*sumaly)/r
        write(*,95) bias,-alfa*radeg
95      format(' bias:',f7.2,' m, slope (roll0):',f6.2,' degrees')
c  
        write(*,*) 'cross profiles written on file "scan.cal"'
      endif
c
      close(20)
      end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      real*8 function sqrtx(r)
      implicit double precision(a-h,o-z)
      if (r.lt.0) then
        sqrtx = -sqrt(-r)
      else
        sqrtx = sqrt(r)
      endif
      return
      end
c
      subroutine readhdr(npscan,rangefak,anglefak,timerfak,measid,
     .isync,polarid,dechr0)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c                     r e a d h d r
c
c  reads header in Riegl 2dd-data from unit 10
c  supports both new Q120/Q240 header and old Q140 headers
c  old scanner headers given by measid = 140
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      implicit real*8(a-h,o-z) 
      integer*1 hdr(47),protocolid,headerid,polarid,syncflag,i1
      integer*4 wsize,i4
      integer*2 wsize1,wsize2,i2
      integer*2 bytespl,meassize,meascount,measid
      character*8 char8,ctime
      character*32 cepoch
      real*4 rangeunit,angleunit,timerunit
      equivalence (hdr(1),wsize),
     .(hdr(5),bytespl),(hdr(11),meassize),
     .(hdr(13),meascount),(hdr(35),rangeunit),(hdr(39),angleunit),
     .(hdr(43),timerunit) 
c
      read(10,err=90) (hdr(j),j=1,47)
      if (wsize.le.47) goto 80
c
c   Q240 new header
c   ---------------
c
      rewind(10)
      read(10,err=90) wsize,bytespl,protocolid,headerid
      write(*,15) wsize,bytespl,protocolid
15    format(/' RIEGL Q120/220 HEADER bytes:',i4,', bytes pr scan:',i5,
     .', protocol id:',i2) 
      if (headerid.ne.10) then
        write(*,*) '*** error in header: headerid is not 10:',headerid
        stop
      endif
      read(10,err=90) i2,meassize,meascount,i1,i2,i1,measid,i1,i2,i1,i2
      write(*,*) 'scanner output: measid ',measid
      if (mod(ishft(measid,0),2).eq.1) write(*,*) '- range' 
      if (mod(ishft(measid,-2),2).eq.1) write(*,*) '- intensity' 
      if (mod(ishft(measid,-3),2).eq.1) write(*,*) '- angle' 
      if (mod(ishft(measid,-4),2).eq.1) write(*,*) '- meas qual' 
      if (mod(ishft(measid,-6),2).eq.1) write(*,*) '- synctimer' 
      if (measid.ne.77) stop '*** only measid 77 is implemented'
      if (i1.ne.8) then
        write(*,*) '*** error in header: parIDmain not 8:',headerid
        stop
      endif
      read(10,err=90) char8,rangeunit,angleunit,timerunit,polarid,
     .(i1,j=1,122),cepoch,ctime,syncflag
c
      if (ctime(1:3).eq.'GPS') then
        read(cepoch,16) ihr,imin,isec
16      format(11x,i2,1x,i2,1x,i2)      
      else
        ihr = 0
        imin = 0
        isec = 0
      endif
      dechr0 = ihr+imin/60.d0+isec/3600.d0
      isync = bytespl
      npscan = meascount
      rangefak = rangeunit
      anglefak = angleunit
      timerfak = timerunit
      write(*,17) meassize,meascount,rangefak,anglefak,timerfak
17    format(' bytes pr meas:',i5,', measurements per scan:',i5/
     .' range,angle and timer units:',4f12.6)
      write(*,*) 'instrument id: ',char8,' polarangle id: ',polarid
      write(*,*) 'time epoch: ',cepoch,' hr,min,sec =',ihr,imin,isec
      write(*,18) ctime,syncflag
18    format(' time source: ',a8,' syncflag: ',b7.7)
      return
c 
c  Q140 old header
c  ---------------
c
80    isync = bytespl
      npscan = meascount
      rangefak = rangeunit
      anglefak = angleunit
      timerfak = timerunit
      write(*,85) wsize,bytespl,meassize,meascount,
     .rangeunit,angleunit,timerunit
85    format(/' RIEGL Q140 HEADER bytes:',i4,', bytes pr scan:',i5,
     .' bytes pr meas:',i5,/' measurements per scan:',i5/
     .' range,angle and timer units:',4f12.6)
      write(*,*)
c
      if (hdr(47).ne.4) 
     .write(*,*) '*** warning: header angleid is not 4 !!'
      measid = 140
      return
90    stop '*** end of file met in header'   
      end
c
      double precision function bilin(ri,rj,a,imax,jmax,iadim)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c                                                                  c
c                           b i l i n                              c
c                                                                  c
c  interpolates values in an array a using bilinear                c
c  (parabolic hyperboloid) interpolation.                          c
c                                                                  c
c  parameters:                                                     c
c                                                                  c
c  bilin       interpolated value                                  c
c                                                                  c
c  ri, rj      interpolation argument, (1,1) in lower left corner, c
c              (imax, jmax) in upper right.                        c
c                                                                  c
c  a           integer*2 array with arguments                      c
c                                                                  c
c  imax, jmax  number of points in grid                            c
c                                                                  c
c  iadim       declared dimension of 'a'                           c
c                                                                  c
c  outside area covered by 'a', the function returns the value of  c
c  the nearest boundary point. if unknown (9999) values are met    c
c  the subroutine returns 9999.                                    c
c                                                                  c
c  programmer:                                                     c
c  rene forsberg, july 1983, modified june 1989                    c
c                                                                  c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      implicit double precision(a-h, o-z)
      real*4 a
      dimension a(iadim)
c
      in = ifrac(ri)
      ie = ifrac(rj)
      rn = ri - in
      re = rj - ie
c
      if (in.lt.1) go to 1
      go to 2
 1      in = 1
        rn = 0.0
      go to 5
 2    if (in.ge.imax) go to 3
      go to 5
 3      in = imax-1
        rn = 1.0
 5    continue
      if (ie.lt.1) go to 4
      go to 7
 4      ie = 1
        re = 0.0
       go to 6
 7    if (ie.ge.jmax) go to 8
      go to 6
 8      ie = jmax-1
        re = 1.0
 6    continue
c
c  let possible corner value give value if close
c
      k = (in-1)*jmax + ie
      a1 = a(k)
      a2 = a(k+jmax)
      a3 = a(k+1)
      a4 = a(k+jmax+1)
      if (rn.gt.0.0001.or.re.gt.0.0001) goto 10
        if (int(a1).eq.9999) bilin = 9999
        if (int(a1).ne.9999) bilin = a1
        return
10    if (int(a1).ne.9999.and.int(a2).ne.9999.and.
     .int(a3).ne.9999.and.int(a4).ne.9999) goto 11
        bilin = 9999
        return
11    bilin = (1-rn)*(1-re)*a1 +
     .rn*(1-re)*a2 + (1-rn)*re*a3 +
     .rn*re*a4
      return
      end
c
      integer function ifrac(r)
      implicit double precision(a-h,o-z)
      if (r.lt.0) goto 10
      ifrac = r
      return
10    i = r
      if (i.eq.r) goto 20
      ifrac = i-1
      return
20    ifrac = i
      return
      end
