c2345678911234567892123456789312345678941234567895123456789612345678971234567898
      program pix2fits


c  g77 -O -Wuninitialized pix2fits.f -L/star/lib -lcfitsio -lsla -B/usr/lib/x86_64-linux-gnu -B/usr/lib/gcc/x86_64-linux-gnu/4.7 --no-second-underscore -o pix2fits
c
c  Takes co-ordinates and size and extracts pixel data into a FITS file along
c  with header info from the house-keeping file and a `world co-ordinate
c  system'. User has the option to keep raw transmission values, or to 
c  transform to linear intensities (insofar as dynamic range limitations will
c  allow). A FITS binary table extension is written containing the corresponding
c  IAM data, and an unformatted binary index file of pointers to those images
c  in the IAM file is written to a separate file for use by other programs
c  in this suite of software.
c
c  Compile: 
c  % f77 pix2fits.f -L/star/lib `sla_link` `fitsio_link` -o pix2fits -vms (OSF)
c                                                                    -xl 
c                                                                     (Solaris)
c
c  NB: For unswapped data on Suns, change parameter swap to true before
c      compilation.
c
c  Measurement directory specified needs to contain an IAM file in iam.srtrd
c  and a subdirectory mmdata containing the pixel data.
c
c  The maximum size of the extracted array is currently software limited.
c  Parameter INX below is set to enable extraction of 4096x4096 images. If you
c  get a run-time error of the following form:
c 
c  % ./pix2fits
c  2107:pix2fits: /sbin/loader: Fatal Error: set_program_attributes failed to
c  set heap start
c  %
c 
c  on DEC Alpha, or
c
c  % ./pix2fits
c  Killed
c  %
c
c  on SPARC stations, then memory limitation on the machine requires you to set
c  the limit lower - you should reduce parameter INX (and take note that the
c  maximum array size extractable will be smaller) until the code runs;
c  alternatively, see your system administrator to increase your quotas and/or
c  the amount of swap space on the relevant work station.
c
      implicit none
      logical swap,ok,zhflags(2),maxflg,expix,iflg
      parameter (swap=.false.)
c      parameter (swap=.true.)
      character bigim*80,line*80,path*80
      integer i,j,k,l,m,xmid,ymid,size,xs,ys,xe,ye,nstrt,adim
      integer mdim,status,strt,buff(32),length,nscans,nl,buffra,imax
      integer ihmsf(4),idmsf(4),status1
      character sign*1,units*12,objctra*20,objctdec*20
      integer inx,iny,srec,sf,ef,nsf,num,equinox,osize,idum
      parameter (adim=1152,mdim=665536,inx=4096,iny=inx)
      integer*2 i2buff(adim),image(inx,iny),lut(0:32767)
      character rfile*80,iamfile*80,name*30,output*80,idx*80
      real xmin,xmax,ymin,ymax,side,dr
      real tclear,tblack,c,gamma,t,intens
      parameter (tclear=32767.0,tblack=20.0,c=3.193922,gamma=0.3)
      integer ra,dec,points(mdim),idelra,ideldc,i2pi,nef
      double precision xi,eta,xpix,ypix,alpha,pi,msep,x0,y0,maxsep
      double precision iamra,iamdc,xye(2,mdim),xym(2,mdim),coeffs(6)
      double precision reslt,dra,ddec,ra2000,dec2000,bkwds(6),df
      double precision xz,yz,xsc,ysc,perp,orient,pc11,pc22,pc12,pc21
      double precision raref,dcref,crota2,rho1,rho2,s
      parameter (maxsep=5.0d0)
c  FITSIO declarations:
      integer unitno,blocksize,bitpix,naxis,naxes(2),dim1,group
      logical simple,extend
      maxflg=.false.
      num=0
      nef=0
      sf=0
      srec=0
      pi=4.0d0*atan(1.0d0)
      i2pi=nint(pi*2.0d8)
      msep=maxsep*pi/(60.0d0*180.0d0)

c  Calculate t-to-i LUT:
      do i=nint(tblack+1.0),nint(tclear-1.0)
       t=real(i)
       intens=10.0**(gamma*log10((tclear-tblack)/(t-tblack))+c)
       if (nint(intens).lt.32767) then
        lut(i)=intens
       else
        lut(i)=32767
       endif
      enddo
      do i=0,nint(tblack)
       lut(i)=lut(nint(tblack+1.0))
      enddo
      do i=nint(tclear),32767
       lut(i)=lut(nint(tclear-1.0))
      enddo       

      write (*,*) 'Path-name of measurement:'
      read (*,'(a)') path
      write (*,*) 'Extent of scan (xmin,xmax,ymin,ymax mm):'
      read (*,*) xmin,xmax,ymin,ymax
      write (*,*) 'Centre co-ordinates of region to be extracted:'
      read (*,'(a)') line
      write (*,*) 'Equinox of co-ordinates (0=1900, 1=B1950, 2=J2000):'
1     read (*,*) equinox
      if (equinox.lt.0.or.equinox.gt.2) then
       write (*,*) 'Invalid reply: please answer 0, 1 or 2:'
       goto 1
      endif
      write (*,*) 'Size of square image to extract (arcmin):'
2     read (*,*) side
      if (nint(side*60.0/0.6714).gt.inx) then
       write (*,*) 
     :'Error: maximum image size is ',real(inx)*0.6714/60.0
       write (*,*) 'Please try again:'
       goto 2
      endif
      write (*,*) 
     :'Output pixel values in raw transmission (t) or intensity (i)?'
3     read (*,'(a)') units
      if (units(1:1).ne.'t'.and.units(1:1).ne.'i') then
       write (*,*) 'Please answer t or i.'
       goto 3
      endif
      write (*,*) 'Output FITS filename:'
      read (*,'(a)') output

c  Test data:
c      path='/net/cosaxp1/disk1/nch/UKR444'
c      path='/disk1/nch/UKR444'
c      xmin=16.22
c      xmax=338.78
c      ymin=xmin
c      ymax=xmax
c      line='13 34 59.04 -27 30 15.05'
c      equinox=2
c      strt=1
c      side=2.86464
c      units='i'
c      output='testnowcs.fits'

c  Calculate nscans and nl:
      nscans=nint((ymax-ymin)*100.0)
      nl=nscans/adim

      size=nint(side*60.0/0.6714)
      osize=size
      dr=max(real(maxsep)*60.0,side*60.0)

c  Set up the filenames:
      i=length(path)
      bigim=path
      iamfile=path
      rfile=path
      if (path(i:i).eq.'/') then 
       bigim(i+1:i+7)='mmdata/'
       iamfile(i+1:i+9)='iam.srtrd'
       rfile(i+1:i+15)='residual.corrns'
      else
       bigim(i+1:i+8)='/mmdata/'
       iamfile(i+1:i+10)='/iam.srtrd'
       rfile(i+1:i+16)='/residual.corrns'
      endif
      if (units(1:1).eq.'t') then
       units='TRANSMISSION'
       iflg=.false.
      else
       if (units(1:1).eq.'i') units='INTENSITY'
       iflg=.true.
      endif 
      idx=output
      i=length(idx)
      idx(i+1:i+4)='.idx'

c  Try to open file residual.corrns to get x-axis magnification factor (alpha)
c  appropriate to this measurement - if file is not present, default to value
c  as at July 1998:

      open (unit=1,file=rfile,iostat=status,err=5,status='old')
      read (1,*,end=5) alpha
5     close (1)
      if (status.ne.0) then
       write (*,*) 'Warning: failed to find file '
       write (*,'(a)') rfile
       write (*,*) 
     :'Assuming default value for magnification correction:'
       alpha=2.597d-4
       write (*,*) 'alpha = ',alpha
      endif

      strt=1

      nstrt=strt
      call sla_dafin(line,nstrt,reslt,j)
      if (j.lt.0.and.j.gt.-3) then
       write (*,*) 'SLA_DAFIN error decoding RA:'
       write (*,'(a,a)') ' Offending line: '//line
       stop
      else if (j.eq.-3) then
       write (*,*) 'Warning: range check failed for seconds on line'
       write (*,'(1x,a)') line
      endif
      dra=reslt*15.0d0
      call sla_dafin(line,nstrt,reslt,j)
      if (j.lt.0.and.j.gt.-3) then
       write (*,*) 'SLA_DAFIN error decoding DEC:'
       write (*,'(a,a)') ' Offending line: '//line
       stop
      else if (j.eq.-3) then
       write (*,*) 'Warning: range check failed for arcseconds on line'
       write (*,'(1x,a)') line
      endif
      write (name,'(30x)')
      j=min(nstrt-strt+1,30)
      k=min(strt+30-1,nstrt-1)
      name(1:j)=line(strt:k)
      ddec=reslt
      if (equinox.eq.0) then
       call sla_preces('FK4',1900.0d0,1950.0d0,dra,ddec)
       call sla_fk45z(dra,ddec,1950.0D0,ra2000,dec2000)
       ra=nint(ra2000*1.0d8)
       dec=nint(dec2000*1.0d8)
      else if (equinox.eq.1) then
       call sla_fk45z(dra,ddec,1950.0D0,ra2000,dec2000)
       ra=nint(ra2000*1.0d8)
       dec=nint(dec2000*1.0d8)
      else
       ra=nint(dra*1.0d8)
       dec=nint(ddec*1.0d8)
      endif

c  Set up the 0h flags: if the requested area crosses the boundary at RA=0h,
c  then set the flag to true:
      df=dble(dr)*pi*1.0d8/(180.0d0*3600.0d0)
      idelra=nint(df/cos(dble(dec)/1.0d8))
      if ((ra+idelra).ge.nint(2.0d8*pi)) then
       zhflags(1)=.false.
       zhflags(2)=.true.
      else if ((ra-idelra).lt.0) then
       zhflags(1)=.true.
       zhflags(2)=.false.
      else 
       zhflags(1)=.false.
       zhflags(2)=.false.
      endif

c      goto 500
      status1=0
      open (unit=1,file=iamfile,status='old',access='direct',
     :form='unformatted',recl=128)
      ideldc=nint(dble(dr)*pi*1.0d8/(180.0d0*3600.0d0))
      i=0
30    i=i+1
      read (1,rec=i,iostat=status1) buff
      if (status1 .ne.0) goto 40
      if (swap) call swbyte4(buff,32)
      if (buff(29).lt.0) goto 30
c  No need to continue if less than lowest or greater than highest RA:
      idelra=nint(df/cos(dble(buff(2))/1.0d8))
      if ((buff(1)+idelra).lt.ra.and..not.zhflags(2)) goto 30
      if ((buff(1)-idelra).gt.ra.and..not.zhflags(1)) goto 30
c  Note: every time a comparison is made between centre and IAM ra, must
c  check to see if chart straddles 0h boundary and tweak IAM ra so as to 
c  include any objects on the other side of the boundary to the chart centre:
      if (zhflags(1).and.i2pi-buff(1).lt.idelra) then
       buffra=buff(1)-i2pi
      else if (zhflags(2).and.buff(1).lt.idelra) then
       buffra=buff(1)+i2pi
      else
       buffra=buff(1)
      endif
      if (abs(buff(2)-dec).le.ideldc.and.abs(buffra-ra).le.
     :idelra) then
       num=num+1
       if (num.le.mdim) then
        points(num)=i
       else if (num.eq.mdim+1) then
        write (*,*) 
     :'Warning: maximum number of matches exceeded for object ...'
        maxflg=.true.
       endif
      endif
      goto 30
40    if (maxflg) then
       imax=num
       write (*,*) 'Increase dimension MDIM and recompile -'
       write (*,*) 'maximum number of matches for object is ',imax
      endif
      write (*,*) 'Number of records read from IAM file: ',i-1

       ok=.true.

c  Now derive local transformation and extract the pixel data:
c  Check that there are enough objects found around each image to define a
c  transformation (if there are too many, then output a warning message):
  
       if (num.lt.10) then
        write (*,*) 'ERROR: not enough images found around object.'
        write (*,*) 'Number of images found: ',num
        write (*,*) 'Check that this object lies within the plate'
        write (*,*) 'measurement area; if it does, then suggest you'
        write (*,*) 'increase MAXSEP, recompile and rerun.'
        write (*,*) 'No WCS will be written for this object.'
        ok=.false.
       else if (num.gt.mdim) then
        write (*,*) 'Warning: too many images found around object.'
        write (*,*) 'No. of images found: ',num
        write (*,*) 'Co-ordinate transformation may not yield accurate'
        write (*,*) 'centroiding/WCS for this image.'
        write (*,*) 'Suggest you increase MDIM, recompile and rerun.'
        num=mdim
       endif

c  Define transformation between standard celestial co-ordinates in this
c  region and SuperCOSMOS pixel co-ords adjusted for the magnification error:

c       if (ok) then
        idelra=nint(df/cos(dble(dec)/1.0d8))
        do j=1,num
         read (1,rec=points(j)) buff
         if (swap) call swbyte4(buff,32)
         if (zhflags(1).and.(i2pi-buff(1).lt.idelra)) 
     :buff(1)=buff(1)-i2pi
         if (zhflags(2).and.(buff(1).lt.idelra)) 
     :buff(1)=buff(1)+i2pi
         iamra=dble(buff(1))/1.0d8
         iamdc=dble(buff(2))/1.0d8
         call sla_ds2tp(iamra,iamdc,dble(ra)/1.0d8,dble(dec)/1.0d8
     :,xi,eta,status)
c
c  Rather than using tangent plane projection, use simply delta RA and DEC
c  for the fit:
c
c         xi=(iamra-(dble(ra)/1.0d8))*cos(iamdc)
c         eta=iamdc-(dble(dec)/1.0d8)
c         status=0

         if (status.ne.0.and.ok) then
          write (*,*) 
     :'DS2TP error in tangent plane projection around object ',i
          write (*,*) 'Co-ordinates (radians):'
          write (*,*) iamra,iamdc,ra,dec
          write (*,*) 'Something is wrong somewhere ...'
          write (*,*) 'No WCS will be written for this object.'
          ok=.false.
         else
c  Units: magnified pixels versus degrees
          xye(1,j)=xi*180.0d0/pi
          xye(2,j)=eta*180.0d0/pi
          xpix=(dble(buff(11))/(1.0d3*(1.0d0-alpha)) - 
     :(dble(xmin)*100.0d0 - 0.5d0))
          ypix=dble(buff(12))/1.0d3 - (dble(ymin)*100.0d0 - 0.5d0)
          xym(1,j)=xpix
          xym(2,j)=ypix
c         call sla_dr2tf(3,iamra,sign,ihmsf)
c         call sla_dr2af(2,iamdc,sign,idmsf)
c         write (*,999) (ihmsf(k),k=1,4),sign,(idmsf(k),k=1,4)
c999      format (3i3,'.',i3.3,3x,a1,i2.2,2i3,'.',i2.2)
c         write (*,*) xpix,ypix
c         write (*,*) 'Press return to continue:'
c         read (*,*)
         endif
        enddo
        if (ok) then
         call newfitxy(6,num,xye,xym,coeffs,status)
         if (status.ne.0) then
          write (*,*) 'SLA_FITXY failure: '
          write (*,*) ' J = ',status
          write (*,*) 'No WCS will be written for this object.'
          ok=.false.
         endif
        endif

c        close(1)

c  Now calculate the pixels required to be extracted and extract them:
      
          call sla_invf(coeffs,bkwds,status)
          if (status.ne.0) then
           write (*,*) 'SLA_INVF failure: '
           write (*,*) ' J = ',status
           write (*,*) 'No WCS will be written for this object.'
           ok=.false.
          endif

         call sla_dcmpf(bkwds,xz,yz,xsc,ysc,perp,orient)

         x0=bkwds(1)
         y0=bkwds(4)
         xmid=int(x0)+1
         ymid=int(y0)+1
         expix=.true.
         if (xmid.lt.1.or.ymid.lt.1.or.xmid.gt.nscans.or.ymid.gt.
     :nscans) then
          write (*,*) 'ERROR: centre of requested area falls outside'
          write (*,*) 'plate bounds. No pixel data will be extracted.'
          expix=.false.
         endif
         xs=xmid-size/2
         xe=xs+size-1
         ys=ymid-size/2
         ye=ys+size-1

       if (xs.lt.1.or.xe.gt.nscans.or.ys.lt.1.or.ye.gt.nscans) then
        write (*,*) 'Warning: Part of requested area falls outside'
        write (*,*) 'plate measurement bounds. The largest square'
        write (*,*) 'image within the nearest edge will be extracted,'
        write (*,*) 'and no WCS will be written for this object.'
        xs=max(1,xs)
        ys=max(1,ys)
        xe=min(nscans,xe)
        ye=min(nscans,ye)
        size=min(xe-xs+1,ye-ys+1)
        xe=xs+size-1
        ye=ys+size-1
        ok=.false.
       endif
       call fopen(bigim,xs,ys,adim,99,srec,sf)
       ef=min(adim,sf+size-1)
       do i=ye,ys,-1
        read (99,rec=i) i2buff
        if (swap) call swbyte2(i2buff,adim)
        do j=sf,ef
         if (iflg) then
          m=i2buff(j)
          if (m.ge.0.and.m.le.32767) then
           i2buff(j)=lut(m)
          endif
         endif
         image(j-sf+1,i-ys+1)=i2buff(j)
        enddo
       enddo
       close (99)
       l=ef-sf+1
       k=1
       do while (l.lt.size)
        nsf=ef-sf+1+(k-1)*adim
         call fopen(bigim,xs+adim*k,ys,adim,99,srec,idum)
         do i=ye,ys,-1
          read (99,rec=i) i2buff
          if (swap) call swbyte2(i2buff,adim)
          if (mod(xe,adim).eq.0) then
           nef=adim
          else if (size-l.gt.adim) then
           nef=adim
          else
           nef=mod(xe,adim)
          endif
          do j=1,nef
           if (iflg) then
            m=i2buff(j)
            if (m.ge.0.and.m.le.32767) then
             i2buff(j)=lut(m)
            endif
           endif
           image(nsf+j,i-ys+1)=i2buff(j)
          enddo
         enddo
        close (99)
        l=l+nef
        k=k+1
       enddo

       if (expix) then
     
c  Open the FITS file and write in the header info:
        status=0
        call deletefile(output,status)
        unitno=2
        blocksize=1
        call ftinit(unitno,output,blocksize,status)
        simple=.true.
        bitpix=16
        naxis=2
        naxes(1)=size
        naxes(2)=size
        extend=.true.
        call ftphpr(unitno,simple,bitpix,naxis,naxes,0,1,extend,status)

c  Write the house-keeping information:
        call hkfits(unitno,path,3,status)

c  Write in the WCS information:
        if (ok) then

         call sla_dcmpf(coeffs,xz,yz,xsc,ysc,perp,orient)

c         write (*,*) (coeffs(j),j=1,6)
c         write (*,*) xz,yz
c         write (*,*) xsc,ysc
c         write (*,*) perp,orient

         pc11=cos(perp)*cos(orient)-sin(perp)*sin(orient)
         pc12=cos(perp)*sin(orient)+sin(perp)*cos(orient)
         pc21=sin(perp)*cos(orient)-cos(perp)*sin(orient)
         pc22=sin(perp)*sin(orient)+cos(perp)*cos(orient)
c         pc11=1.0d0/sqrt(2.0d0)
c         pc12=pc11
c         pc21=-pc11
c         pc22=pc11
         raref=dble(ra)*180.0d0/(pi*1.0d8)
         dcref=dble(dec)*180.0d0/(pi*1.0d8)
c  Calculate CROTA2 for older software that doesn't understand the PC rotation
c  matrix (this will only work if the skew is small):
         if 
     :(pc11.eq.0.0d0.or.pc21.eq.0.d0.or.pc12.eq.0..or.pc22.eq.0.0d0)
     : then
          crota2=0.0d0
         else
          s=ysc/xsc
          rho1=atan2(pc21*s,pc11)
          if (rho1.lt.0.0d0) rho1=2.0d0*pi+rho1
          rho2=atan2(-pc12/s,pc22)
          if (rho2.lt.0.0d0) rho2=2.0d0*pi+rho2
          crota2=180.0d0*(rho1+rho2)/(2.0d0*pi)
         endif
 
         call ftpkyd(unitno,'CRVAL1',raref,12,'Axis 1 reference value',
     :status)
         call ftpkyd(unitno,'CRPIX1',x0-dble(xs-1),12,
     :'Axis 1 pixel value',status)
         call ftpkyd(unitno,'CDELT1',xsc,12,
     :'Increment per pixel on axis 1',status)
         call ftpkys(unitno,'CTYPE1','RA---TAN',
     :'Quantity represented by axis 1',status)
         call ftpkyd(unitno,'CRVAL2',dcref,12,'Axis 2 reference value',
     :status)
         call ftpkyd(unitno,'CRPIX2',y0-dble(ys-1),12,
     :'Axis 2 pixel value',status)
         call ftpkyd(unitno,'CDELT2',ysc,12,
     :'Increment per pixel on axis 2',status)
         call ftpkys(unitno,'CTYPE2','DEC--TAN',
     :'Quantity represented by axis 2',status)
         call ftpkyd(unitno,'PC001001',pc11,12,'Axis rotation matrix',
     :status)
         call ftpkyd(unitno,'PC001002',pc12,12,'Axis rotation matrix',
     :status)
         call ftpkyd(unitno,'PC002001',pc21,12,'Axis rotation matrix',
     :status)
         call ftpkyd(unitno,'PC002002',pc22,12,'Axis rotation matrix',
     :status)
         call ftpkyd(unitno,'CROTA2',crota2,12,
     :'DEPRECATED - rotation of axis 2',status)
         call ftpkye(unitno,'EQUINOX',2000.0,6,
     :'Julian reference frame equinox',status)

        endif

        dim1=inx
        group=1
        call ftp2di(unitno,group,dim1,naxes(1),naxes(2),image,status)
        call ftpkys(unitno,'DATATYPE','INTEGER*2','Type of data',
     :status)
        call ftpkys(unitno,'DATUNITS',units,
     :'Units: transmission or "linearised" intensity',status)
        call ftpkyd(unitno,'XPIXELSZ',10.0d0*(1.0d0-alpha),12,
     :'[microns] X pixel size',status)
        call ftpkyd(unitno,'YPIXELSZ',10.0d0,12,
     :'[microns] Y pixel size',status)
        call sla_dr2tf(3,dble(ra)/1.0d8,sign,ihmsf)
        call sla_dr2af(2,dble(dec)/1.0d8,sign,idmsf)
        write (objctra,98) (ihmsf(j),j=1,4)
        write (objctdec,99) sign,(idmsf(j),j=1,4)
98      format (3i3,'.',i3.3)
99      format (a1,i2.2,2i3,'.',i2.2)
        call ftpkys(unitno,'OBJCTRA',objctra,
     :'Centre Right Ascension (J2000)',status)
        call ftpkys(unitno,'OBJCTDEC',objctdec,
     :'Centre Declination (J2000)',status)
        call ftpkyd(unitno,'OBJCTX',x0,12,'[pixels] Centre X on plate',
     :status)
        call ftpkyd(unitno,'OBJCTY',y0,12,'[pixels] Centre Y on plate',
     :status)

c  Now write the IAM parameters to a FITS binary table extension:
        write(6,*)unitno,status
        call ftcrhd(unitno,status)
        call iamtable(unitno,points,num,1,status,swap,side)
 
        call ftclos(unitno,status)
        if (status.gt.0) call printerror(status)

       endif

      close (1)

c  Finally, write the binary index of pointers to a file:

      open (unit=1,file=idx,status='new',access='direct',
     :form='unformatted',recl=1)
      do i=1,num
       write (1,rec=i) points(i)
      enddo
      close (1)

      end

      subroutine fopen(path,xs,ys,lwidth,bd,srec,sf)
      implicit none
      character path*(*),name*3,fname*80,fstr*12
      integer xs,ys,lwidth,bd,srec,sf,n,rlen,length,i

      if (mod(xs,lwidth).eq.0) then
       n=(xs/lwidth)
      else
       n=(xs/lwidth)+1
      endif
      if (n.le.9) then
       write (name,'(a1,i1,a1)') 'l',n,' '
      else
       write (name,'(a1,i2)') 'l',n
      endif
      i=length(path)
      write (fstr,100) i
100   format('(a',i2,',a3)')
      
      write (fname,fstr) path,name
      rlen=lwidth*2

      open (unit=bd,file=fname,status='old',access='direct',
     :form='unformatted',recl=rlen)
      srec=ys
      sf=mod(xs,lwidth)
      if (sf.eq.0) sf=lwidth
      return
      end
       
      include 'swbyte2.f'
      include 'swbyte4.f'
      include 'newfitxy.f'
      include 'length.f'
      include 'deletefile.f'
      include 'printerror.f'
      include 'hkfits.f'
      include 'iamtable.f'
