      program readscan
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c                      R E A D S C A N
c
c  Reads Riegl 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
c  input: 
c
c  scanfile 
c  coorfile
c  geoidfile   (0: no geoid file)
c  outfile1    (0: no scanfile, only vert.)
c  outfile2
c  thr,min,sec,
c  tshift, tscale 
c  t1, t2 
c  lat1 lat2 lon1 lon2
c  nave, ithin, lbin, lbino, 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  tshift       empirical time shift of lidar data (dechr)
c  tscale       empirical time stretch factor (should be 1 except for ".log"-files)
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  lbin         t=binary format (.2dd)    
c  lbino        t=output in binary format, f=txt
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  minamp       amplitude threshold, the minimum signal strength 
c               (amplitude) for the received return pulse
c               std value min. 10, maybe more over ice
c  rangemin     min. range std. 10 m 
c               data of laser range below rangemin are excluded,
c
c  outfile1: thr,lat,lon,ground_h,amplitude,iscan,n_ave for all pts
c  outfile2: do + flightelev,range, for vertical (possibly averaged) scan only
c
c  Note: Riegl data has occasional irregularities which means time can slip 
c  typically 1 sec. Program repairs some errors, but data should always
c  be checked. For a quick view of contents of file use "rieg2txt".
c  
c  (c) Rene Forsberg, DTU-Space/Denmark, Nov 2000
c
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  updated: vertical output, arrays extended
c           added thresholds: min range, scan angle interval, min aplitude
c           time spike test deleted  l361-366                krk 18 nov 04
c  updated: corrected interpolation of heading when interpolating
c           between -179 and 179 degrees             las 19 nov 04 
c  updated: timer routine changed (same routine in all places) las 22-24 nov 04
c  updated: scan angle interval corrected las 26 nov 04
c  updated: scans missing when thin > 1 las 17 jan 05
c  updated for gaps in gps/ins files (skipped if more than 5 sec), rf feb 08
c  updated for ".log" 5-byte format, rf oct 2011 - WARNING: no clock info in 5-byte format,
c  - use only small batches as time scale will slip with lost scans in logging
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      implicit double precision(a-h,o-z)
      character*36 ifile,pfile,gfile,ofile1,ofile2
      logical lnew,lbin,lbino,lpr,lcalib,lnewscan,lgeo,lgeoid,
     .lscan,lspike,lwarn
      dimension gpst(488010),rlat(488010),rlon(488010),rh(488010),
     .pitch(488010),roll(488010),hdg(488010)
      real*4 ggrid(500000),rvel(488010)
      integer*1 sync(2)
      integer*2 syncseq
      equivalence (sync(1),syncseq)
      integer*1 iamp,iampo,iiscano
      integer*1 brange(4),bang(4),btime(4)
      integer*4 irange,iang,itime
      equivalence(brange(1),irange),(bang(1),iang),(btime(1),itime) 
c
      nmax = 488010
      igdim = 500000
c
      irange = 0
      iang = 0
      itime = 0
      lnew = .true.
      nsum = 0
      numamp = 0
      numsc = 0
      numrange = 0
      lwarn = .false.
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(*,*) tshift, tscale
      read(*,*) t1,t2
      read(*,*) rlat1,rlat2,rlon1,rlon2
      read(*,*) nave,ithin,lbin,lbino,lpr,lcalib,lgeo
      read(*,*) dxl,dyl,dzl
      read(*,*) pitch0,roll0,hdg0
      read(*,*) minamp, rangemin
c
      if (lbin) then
        open(10,file=ifile,form='unformatted',
     .  access='direct',recl=1)
      else
        open(10,file=ifile,status='old')
      endif
      lscan = (ofile1.ne.'0')
      if (.not.lscan) lcalib = .false. 
      open(12,file=pfile,status='old')
      if (lscan) then
        if (lbino) then
          open(20,file=ofile1)
          close(20,status='delete')
          open(20,file=ofile1,form='unformatted',access='direct',
     .    recl=18)
          ireco = 1
        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')
      tstart = thr + min/60.d0 + sec/3600.d0 + tshift
      if (t1.lt.tstart) t1 = tstart
      i0 = 1
c
      write(*,4) 
4     format(' READSCAN - output laser scanner surface points')
      if (.not.lscan) write(*,*) '*** only vertical output'
      if (lbin) write(*,*) 'Riegl binary format (.2dd) assumed'
      if (.not.lnew) write(*,*) '*** old 2000-scanner assumed ***'
      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) stop '*** geoid grid too large, sorry'
        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)
      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  three cases: 2000 text input
c               new text input - with GPS time tag
c               binary input - do
c
c  scanner installation orientation:
c  iscandir = 1 for scan from left to right, -1 right to left 
c  dx dy dz is offset vector from GPS antenna to laser in acft body system
c  x forward  y right  z down
c
      iscandir = -1
c
      timer0 = 0
      igpssec = 0
c
      write(*,*)
      write(*,*) 'Reading scanner data from: ',ifile
c
c  read riegl .2dd binary file header
c  ----------------------------------
c
      if (lbin) then
        call readhdr(npscan,rangefak,anglefak,timerfak,nbytes,irec)
        iscan = 0
c
c  40 hz assumed for spike checks
c
        ihz = 40
        dt = 1.d0/ihz/(4*npscan)
        dtscan = 1.d0/ihz
        dtscan1 = dtscan - npscan*dt
c 
c time for first scan added by las 24 nov 2004
c
        read(10,rec=irec,err=90) sync(1)
        read(10,rec=irec+1,err=90) sync(2)
c        
        irec = irec+2
        if (nbytes.eq.8) then
          read(10,rec=irec+5,err=90) btime(1)
          read(10,rec=irec+6,err=90) btime(2)
          read(10,rec=irec+7,err=90) btime(3)
          timer = itime*timerfak
        else
c
c  ".log"-files: no timing information, set time counter to zero 
c
          timer = 0
          timer0 = 0
          igpssec = 0
        endif
        thr = tstart + (igpssec+timer)/3600
        irec = irec + npscan*nbytes     
c
        if (syncseq.ne.npscan*nbytes) then
        write(*,*) 
     .    '*** warning - syncseq not ok',syncseq,npscan*nbytes,nscan
        else
          write(*,15) thr
15        format(' --- First data sync ok,', 
     .    ' time of first scan in file:',f11.5)
          write(*,*) 
        endif
c
        do while (thr.lt.t1)
          read(10,rec=irec,err=90) sync(1)
          read(10,rec=irec+1,err=90) sync(2)
          irec = irec+2
          if (syncseq.ne.npscan*nbytes) 
     .    write(*,*) '*** warning - syncseq recs not ok',
     .    syncseq,npscan*nbytes,' scan no:',iscan
c         
          if (nbytes.eq.8) then 
            read(10,rec=irec+5,err=90) btime(1)
            read(10,rec=irec+6,err=90) btime(2)
            read(10,rec=irec+7,err=90) btime(3)
            timer = itime*timerfak
          else
c  time for scan: 40 Hz, calibrated by 2002 data
            timer = timer + 0.025002d0*tscale
          endif
c
          iscan = iscan+1
c
c  check for timer errors in initially skipped scans
c
          if (timer.lt.timer0) then
            igpssec = igpssec+1
            if (timer0+dtscan-timer.lt.0.97)
     .      write(*,16) iscan,timer,timer0+dtscan
16          format(' *** timer error in scan start:',i7,' ',2f11.5)
          endif
c
          thr = tstart + (igpssec+timer)/3600
          timer0 = timer
          irec = irec + npscan*nbytes
        enddo
      endif
      write(*,17) iscan
17    format(' Number of scans skipped in start of file:',i7)
c
c  parameters for old scanner (2000 flights)
c  data with 1 sec break every one minute
c  10 scans/second, 396 points/scan in 0.033 sec
c
      if (.not.lnew) then 
        dt = 0.033d0/396 / 3600
        dtscan = 0.067d0 / 3600
        tbreak = 1
      endif
c
c  loop initialization parameters
c  iscan is counting scan no
c  iiscan is point no within scan
c
      thr = tstart - dtscan
      throp = 0
      angp = 999
      nscan = 0
      nlas = 0
      iiscan = 1
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
      timer1 = timer
c
c     main scanner loop entry 
c     =======================
c
20    if (lbin) then
c
c  read binary .2dd data
c  ---------------------
c  byte order is dos .. modify for unix/sun use
c
        if (iiscan.eq.1) then
          read(10,rec=irec,err=90) sync(1)
          read(10,rec=irec+1,err=90) sync(2)
          irec=irec+2
          if (syncseq.ne.npscan*nbytes) 
     .    write(*,*) 
     .    '*** warning - syncseq not ok',syncseq,npscan*nbytes,nscan
          nscan = nscan+1
          iscan = iscan+1
          lnewscan = .true.
        endif
        if (nbytes.eq.8) then 
          read(10,rec=irec,err=90)   brange(1)
          read(10,rec=irec+1,err=90) brange(2)
          read(10,rec=irec+2,err=90) iamp
          read(10,rec=irec+3,err=90) bang(1)
          read(10,rec=irec+4,err=90) bang(2)
          read(10,rec=irec+5,err=90) btime(1)
          read(10,rec=irec+6,err=90) btime(2)
          read(10,rec=irec+7,err=90) btime(3)
          irec = irec+8
c        
          timer = itime*timerfak
        elseif (nbytes.eq.5) then
          read(10,rec=irec,err=90)   brange(1)
          read(10,rec=irec+1,err=90) brange(2)
          read(10,rec=irec+2,err=90) iamp
          read(10,rec=irec+3,err=90) bang(1)
          read(10,rec=irec+4,err=90) bang(2)
          irec = irec+5
c
c  time for scan: 40 Hz, calibrated by 2002 data
c
          if (iiscan.eq.1) timer1 = timer1 + 0.025002d0*tscale
          timer = timer1 + (iiscan-1)/207*0.00828d0*tscale
        else
          stop '*** header nbytes wrong'
        endif
c
        range = irange*rangefak
        vgon = 2*mod(iang,10000)*anglefak
        if (timer.lt.timer0) then
          igpssec = igpssec+1
          if (timer0+dtscan-timer.lt.0.97) 
     .    write(*,25) iscan,timer,timer0+dtscan
25        format(' *** timer error in scan start:',i7,2f9.5)
        endif
c	    
        timer0 = timer
        thr = tstart + (igpssec+timer)/3600
c
c  check output first scan timer info
c  skip scans if ithin > 1
c
        if (iiscan.eq.1) then
          if (nscan.eq.1) then
            write(*,26) iscan,igpssec,timer,thr
26          format
     .      (' First scan used:',i6,', igpssec timer thr:',
     .      i6,' ',f9.2,' ',f11.7)
            write(*,*) 'Centre of scan data every 30 sec: '
            write(*,*)
          endif
c
          if (ithin.gt.1.and.mod(nscan,ithin).ne.1) then
            irec = irec + (npscan-1)*nbytes
            iscan = iscan+1
	      iiscan = 1
            goto 20
          endif 
        endif
c
        iiscan = iiscan+1
        if (iiscan.gt.npscan) iiscan = 1
      elseif (lnew) then
c
c  read text format scan input
c  NB: not checked recently!
c  ---------------------------
c
        read(10,*,err=90) iiscan,range,iamp,vgon,timer 
        if (iiscan.eq.1) then
          nscan = nscan+1
          iscan = iscan+1
          lnewscan = .true.
        endif
        if (timer.lt.timer0) igpssec = igpssec+1
        timer0 = timer
        thr = tstart + (igpssec+timer)/3600
c
        if (ithin.gt.1) then
          if (mod(nscan,ithin).ne.1) goto 20
        endif 
      else
c
c  read old 2000 text format
c  -------------------------
c
        read(10,*,err=90) iiscan,range,iamp,vgon 
        if (iscan.eq.1) then
          thr = thr + dtscan
          nscan = nscan+1
          iscan = iscan+1
          lnewscan = .true.
c  make empirical 1 sec break every 1 min
          if (mod(nscan,600).eq.0) thr = thr + tbreak/3600
        else  
          thr = thr + dt
        endif
c
        if (ithin.gt.1) then
          if (mod(nscan,ithin).ne.1) goto 20
        endif 
      endif
c
c  exit loop for t > t2, skip laser point for too small ranges
c  -----------------------------------------------------------
c
      if (thr.gt.t2) goto 90
      if (range.lt.1.5) goto 20
c
c  only use data with accepted amplitude
c
      if (iamp.lt.minamp) then
        numamp=numamp+1      
        goto 20
      endif
c
c  only use data with laser range more than rangemin
c
      if (range.le.rangemin) then
        numrange=numrange+1
        goto 20
      endif
c
c  interpolate gps pos at time thr
c  -------------------------------
c  find which index igps,igps+1 to use in GPS file
c
30    if (gpst(igps+1).lt.thr) then 
        igps = igps+1
        if (igps.gt.ngps-1) then 
          igps = ngps-1
        else
          goto 30
        endif
      endif
      ddt = gpst(igps+1)-gpst(igps)
      if (ddt.le.0.or.ddt*3600.gt.5.1) then
        if (.not.lwarn) write(*,31) gpst(igps), thr
31      format(' *** scan skipped, time gap in gps at time: ',2f9.4)
        lwarn = .true.
        goto 20
      endif
      lwarn = .false.    
      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 20
      if (alon.lt.rlon1.or.alon.gt.rlon2) goto 20
c
      if (ah.lt.ahmin) ahmin = ah
      if (ah.gt.ahmax) ahmax = ah
      if (range.lt.rmin) rmin = range
      if (range.gt.rmax) rmax = range 
      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 correct interpolation between -179 and 179
c 
      if ( hdg(igps)*hdg(igps+1)>-8100.0 ) then
        ahdg = (1-r)*hdg(igps) + r*hdg(igps+1) + hdg0
      else
        if ( hdg(igps)<-90.0 ) then 
          ahdg = (1-r)*(hdg(igps)+360.0) + r*hdg(igps+1) + hdg0
        else 
          ahdg = (1-r)*hdg(igps) + r*(hdg(igps+1)+360.0) + hdg0
        endif
        if(ahdg>180.0) ahdg=ahdg-360.0        
      endif
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-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 + dyl
        dz = cos(ang)*range + 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/vel
        glon = - sin(ang+aroll/radeg)*range/vel
        gh = ah - cos(apitch/radeg)*cos(ang+aroll/radeg)*range - 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 (lnewscan) then
	  thrsum = thr
          glatsum = glat
          glonsum = glon
          ghsum = gh
          ahsum = ah
          angsum = ang
          ampsum = iamp 
          nsum = 1
	    i0 = 1 
	    lnewscan = .false.
	  goto 20
	endif
c
c next 1 line changed from if (iiscan-i0.lt.nave) then 17 jan 05 by las
c
        if ((iiscan-i0).lt.nave) then
          thrsum = thrsum + thr
          glatsum = glatsum + glat
          glonsum = glonsum + glon
          ahsum = ahsum + ah
          ghsum = ghsum + gh
          angsum = angsum + ang
          ampsum = ampsum + iamp
          nsum = nsum+1
          goto 20
        endif
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
        endif
        nsumo = nsum
c
        thrsum = thr
        glatsum = glat
        glonsum = glon
        ghsum = gh
        ahsum = ah
        angsum = ang
        ampsum = iamp 
        nsum = 1
        i0 = iiscan
        if (nsumo.lt.1) goto 20
c
        glat = glato
        glon = glono
        gh = gho
        ah = aho
        ang = ango
        iamp = iampo
      endif
c           
c  output of results
c  -----------------
c
      iampo = iamp
c
c  swath data - binary output do not contain nsumo
c
      if (lscan) then
        if (lbino) then
          ithro = nint(thro*10000000)
          iglat = nint(glat*10000000)
          iglon = nint(glon*10000000)
          igh = nint(gh*1000)
          iiscano = iiscan
          write(20,rec=ireco) ithro,iglat,iglon,igh,iampo,iiscano
          ireco = ireco+1
        else
          if (nave.le.1) then
            write(20,39) thro,glat,glon,gh,iampo,iiscan
          else
            write(20,39) thro,glat,glon,gh,iampo,iiscan,nsumo
          endif
        endif
      endif
39    format(f9.6,f10.6,f12.6,f8.2,1x,i3,i4,i3)
c
c  vertical data
c
      if (angp.ge.0.and.ang.lt.0) then
        if (nscan.eq.1.or.thro-throp.gt.0.0083333) then
          write(*,40) iscan,thro,glat,glon,gh,iampo
          throp = thro
        endif
        write(21,41) thro,glat,glon,gh,iampo,ah,range
      endif
40    format(' ',i6,f11.7,f10.6,f11.6,f8.2,1x,i3,i3,3i5)
41    format(' ',f11.7,f10.6,f12.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,411) alfa,gh
411     format(2f10.2)
      endif
c
c  end of scanner loop
c  -------------------
      goto 20
c
90    write(*,*)
      write(*,*) 'Total number of scans read: ', nscan
      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(*,*) 
     .'scanpt data (thr,lat,lon,h,ampl,i) output on: ',ofile1
      write(*,*) 
     .'center pts (thr,lat,lon,h,ampl,i,airh,range): ',ofile2
c
c
      write(*,*) 'no of rejected data due to amplitude:', numamp
      write(*,*) 'no of rejected data due to range:', numrange
      if (nbytes.eq.5) then 
        write(*,*)
        write(*,*)
     .  '** WARNING: 5-byte Riegl records, time scale may slip **'
      endif
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,nbytes,
     .irec)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c                     r e a d h d r
c
c  reads header in Riegl 2dd-data from unit 10
c  will accept unknown headers from greenland 2002 flights 
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      implicit real*8(a-h,o-z) 
      integer*1 hdr(47)
c      character*4 csize
      integer*4 wsize
      integer*2 wsize1,wsize2
      integer*2 bytespl,meassize,meascount
      real*4 rangeunit,angleunit,timerunit
      equivalence (hdr(1),wsize1),(hdr(3),wsize2),
     .(hdr(5),bytespl),(hdr(11),meassize),
     .(hdr(13),meascount),(hdr(35),rangeunit),(hdr(39),angleunit),
     .(hdr(43),timerunit) 
c
      irec = 1
      do j = 1,47
        read(10,rec=irec,err=90) hdr(j)
        irec = irec+1
      enddo
c
      wsize = wsize2*256**2 + wsize1
      if (wsize.eq.538976288) goto 88
      npscan = meascount
      nbytes = meassize
      rangefak = rangeunit
      anglefak = angleunit
      timerfak = timerunit
c
      write(*,15) wsize,bytespl,meassize,meascount,
     .rangeunit,angleunit,timerunit
15    format(/' RIEGL 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 !!'
      irec = 48
      return
c
c  year 2002 unknown headers
c
88    write(*,*) '- no 2dd-header, default 2002 params assumed -'
      npscan = 208
      rangefak = 0.008
      anglefak = 0.01
      timerfak = 0.00001
      irec = 1675
      return
c
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
