program spectext

c Program is a refined version of the program ttfspect which itself was a
c more sophisticated version of bandspectra which originally grew out of
c bandsort. The program here:
c
c 	- groups individual observations of an object together
c	- sorts object observations in order of increasing wavelength
c 	- discards objects with dud (negative) fluxes
c	- combines the SExtractor flags for each observation of an object
c	  into a single number
c	- searches for and flags calibration stars in the field and uses
c	  their raw fluxes to produce mean normalisation factors at each
c	  wavelength
c	- assigns fresh ID numbers to each object
c	- calculates wavelength correction based on distance of object
c	  from optical axis 
c	- chooses fixed aperture size that optimizes S/N for each object
c	- normalises individual observations to calibrated stellar spectra
c	- converts measured counts to physical flux units
c	- performs a weighted least-squares fit to each spectrum
c	- calculates an emission-line quality factor based on the maximum
c	  change in flux between two adjacent slices
c
c
c Input file name is 'spectext.in'
c
c Input format requires the following 16 columns:
c     slice no., SExt ID no., X, Y, flux(ap1), err(ap1), flux (ap2), err(ap2)
c     flux(ap3), err(ap3), flux(ap4), err(ap4), ellipticity, FWHM,
c     star/gal classification, SExt flag
c
c Output format is one object per line, with each wavelength
c observation in a separate column. Other fields such as continuum slope
c and intercept, FWHM, star/galaxy classification, ellipticity and
c aperture size, are also printed. Missing flux measurements are padded
c out with the mean detection threshold for that particular aperture.
c
c Assumes stellar calibration data are sorted L to R in order of stars
c with increasing X and X,Y being the first two entries in each column.
c
c DHJ 9/9/97


c **********************************************************************
c *** N.B.: *** Must change format statement in fortran lines 80/100 *** 
c ************* to accommodate number of flux slices, manually !!    ***
c **********************************************************************


c Max number of observations = 60000 (must be > N)
c Max number of objects = 10000
c Max number of slices = 12 (must be > Nsl)
c Max number of calibration stars = 3 
c Number of observation fields (input) = 16


      implicit none
      character*50 inFile, outFile1, outFile2, outFile3, inStarFile
      character*50 inWaveFile, inApertFile
      integer col, row, inpRow, N, numObj, numObsv, outRow, objRow, Nsl
      integer inpStart, bestCol, Nstars, nextstar, report, numReport
      integer refRow, rowDepth, distflag1, distflag2, maxRowDepth, k

      real input(60000,16), objbuffer(12,16), meanRow(16)
      real apDiam(4), bestSN, currSN, scaleFactors(12,4)
      real meanFactors(12), centwaveL(12), sigma(4,(12+1)), meanSigma(4)
      real output1(60000,16), output2(10000,25), stardata((12+2),(3+1))
      real corrWaveL(12), flux(12), delF(12), t(12)
      real Xprev(3*12), Yprev(3*12)
      real dist, maxDist, Xref, Yref, X, Y, Xcen, Ycen, realRow
      real Dpix, beta, Rpix2, placeholder, lowSN, starflag, starDist
      real gn, Area, sig2, siglevel, S, Sx, Sy, Stt, a, b, numer, denom
      real maxChange, currChange, uncert, f1, f2, corrWave1, corrWave2
      real delf1, delf2, fitflux1, fitflux2, contSub1, contSub2
      real distXY, distXYrefs, freq, singleFactor, relErrLimit


	maxRowDepth = 100
c --- maximum search distance allowed (in rows) from reference obsvn

	k = 1
c --- maximum comparison distance to go when checking current obsvn
c --- against those previous (in units of Nsl; i.e. k=3 => dist=3*Nsl)

	freq = 500.0
c --- frequency of input lines to report on progress through input
c --- catalogue (freq=500.0 means progress reported every 500th line)


c      inFile = 'tsex.out'
c      inStarFile = 'test_stars2.txt'
c      inWaveFile = 'ac106waveL.txt'
c      inApertFile = 'ac106aps.txt'
c      N = 190
c      Nsl = 7
c      maxDist = 5.0
c      starDist = 3.0
c      gn = 1.36
c      Dpix = 1024
c      beta = 0.14
c      Xcen = 520
c      Ycen = 308
c      Nstars = 3
c      siglevel = 3
c      report = 8
c      outFile1 = 'jnk1'
c      outFile2 = 'jnk2'


      open(unit=1, file='spectext.in', status='old')

      read(1,*,end=10,err=10) inFile
      read(1,*,end=10,err=10) inStarFile
      read(1,*,end=10,err=10) inWaveFile
      read(1,*,end=10,err=10) inApertFile
      read(1,*,end=10,err=10) N
      read(1,*,end=10,err=10) Nsl
      read(1,*,end=10,err=10) maxDist
      read(1,*,end=10,err=10) starDist
      read(1,*,end=10,err=10) gn
      read(1,*,end=10,err=10) Dpix
      read(1,*,end=10,err=10) beta
      read(1,*,end=10,err=10) Xcen
      read(1,*,end=10,err=10) Ycen
      read(1,*,end=10,err=10) Nstars
      read(1,*,end=10,err=10) siglevel
      read(1,*,end=10,err=10) report
      read(1,*,end=10,err=10) relErrLimit
      read(1,*,end=10,err=10) outFile1
      read(1,*,end=10,err=10) outFile2
      read(1,*,end=10,err=10) outFile3




      write (6,*) '  '
      write (6,*) 'Input (near-sorted) catalogue: ',inFile
      write (6,*) 'Input stellar calibration file name: ',inStarFile
      write (6,*) 'Input central wavelengths file: ',inWaveFile
      write (6,*) 'Input aperture sizes/noise file: ',inApertFile
      write (6,*) 'Number of input lines (observations) to be read: ',N
      write (6,*) 'Number of wavelength slices: ',Nsl
      write (6,*) 'Maximum pix-dist for an object match: ',maxDist
      write (6,*) 'Maximum pix-dist for calibn star match: ',starDist
      write (6,*) 'Gain of the CCD (e-/ADU): ',gn
      write (6,*) 'Side length (pix) of CCD: ',Dpix
      write (6,*) 'Angle subtended by CCD: ',beta
      write (6,*) 'X-coord of optical axis centre: ',Xcen
      write (6,*) 'Y-coord of optical axis centre: ',Ycen
      write (6,*) 'Number of calibration stars in this field: ',Nstars
      write (6,*) 'Nominal detection threshold (in sigma): ',siglevel
      write (6,*) 'Threshold line flux (cts); object count: ',report
      write (6,*) 'Limiting precision of scaling: ',relErrLimit
      write (6,*) 'Output catalogue 1 name: ',outFile1
      write (6,*) 'Output catalogue 2 name: ',outFile2

 10    continue


c --------------------------------------------------------------------
c     Initialising input array:
c --------------------------------------------------------------------

      write (6,*) ' '
      write (6,*) 'Initialising large object arrays...'
      do row=1,N
	do col=1,16
		input(row,col)=0.0
	enddo
      enddo



c --------------------------------------------------------------------
c     Initialising output arrays:
c --------------------------------------------------------------------


      do row=1,N
	do col=1,16
		output1(row,col)=0.0
	enddo	
      enddo

      do row=1,N
	do col=1,(Nsl+13)
		output2(row,col)=0.0
	enddo	
      enddo




c --------------------------------------------------------------------
c     Reading input file and filling input array:
c --------------------------------------------------------------------


      write (6,*) 'Reading input data...'
      write (6,*) ' '

      open(unit=11, file=inFile, status='old')
      do row=1,N
	read(11,*,end=20,err=20) (input(row,col),col=1,16)
      enddo
20    continue

c --- N.B: Any left over space is padded with 0.0s.




c --------------------------------------------------------------------
c     Reading stellar calibration data:
c --------------------------------------------------------------------


      open(unit=14, file=inStarFile, status='old')
      do row=1,(Nsl+2)
	read(14,*,end=30,err=30) (stardata(row,col),col=1,Nstars)
      enddo
30    continue



c --------------------------------------------------------------------
c     Reading central wavelength data:
c --------------------------------------------------------------------


      open(unit=15, file=inWaveFile, status='old')
      do row=1,Nsl
	read(15,*,end=40,err=40) centwaveL(row)
      enddo
40    continue




c --------------------------------------------------------------------
c     Reading aperture sizes/noise data:
c --------------------------------------------------------------------


      open(unit=16, file=inApertFile, status='old')
      do row=1,4
	read(16,*,end=50,err=50) apDiam(row),(sigma(row,col+1),col=1,Nsl)
      enddo
50    continue




c --------------------------------------------------------------------
c     Calculating mean sigma values for each aperture:
c --------------------------------------------------------------------

      do row=1,4
        meanSigma(row) = 0.0
      	do col=1,Nsl
	  meanSigma(row) = meanSigma(row)+sigma(row,col+1)
        enddo
        meanSigma(row) = meanSigma(row)/Nsl
      enddo



c =====================================================================
c     Object Loop:

      inpRow = 1
      inpStart = 1
      numObj = 0
      outRow = 1
      numObsv = 0
      nextstar = 1
      starflag = 0
	write (6,*) ' ------------------------------------------------'


c --- Copy X and Y positions and input row number for reference; also
c --- initialise previous values of (Xref,Yref).
	
	do row=1,Nsl
		Xprev(row)=0.0
		Yprev(row)=0.0
	enddo
	Xref = input(inpRow,3)
	Yref = input(inpRow,4)
	refRow = inpRow


c --- New object: initialise object buffer array:

      	do row=1,Nsl
		objbuffer(row,1)=row
		do col=2,16
			objbuffer(row,col)=0.0
		enddo
      	enddo

c --- Copy first entry to object buffer:

	objRow = int(input(inpRow,1))
	do col=1,16
		objbuffer(objRow,col)=input(inpRow,col)
	enddo
	inpRow=inpRow+1
	numObsv = numObsv+1




c -------------------------------------------------------------
      do while (inpRow.le.(N+1))



	X = input(inpRow,3)
	Y = input(inpRow,4)
	rowDepth = inpRow - refRow
	dist = sqrt((X-Xref)**2+(Y-Yref)**2)
c --- i.e. calculating distances between current observation and
c --- current (Xref,Yref), as well the distance from the reference row.



c --- Checking current observation against previous ones:
	distflag1 = 0
	distflag2 = 0
c --- means that neither the current (X,Y) or (Xref,Yref) are nearby
c --- any of the Nsl previous (Xref,Yref) points.
	row=1
	do while (row.le.(k*Nsl))
		distXY = sqrt((X-Xprev(row))**2+(Y-Yprev(row))**2)
		distXYrefs = sqrt((Xref-Xprev(row))**2+(Yref-Yprev(row))**2)
		if (distXY.lt.maxDist) then
			distflag1 = 1
		endif
c --- i.e. current (X,Y) is no good for adding to object buffer since it
c --- is nearby one of the Nsl previous (Xref,Yref).
		if (distXYrefs.lt.maxDist) then
			distflag2 = 1
		endif
c --- i.e. current (Xref,Yref) is no good to write to either output file
c --- since it is nearby one of the Nsl previous (Xref,Yref).
		if ((distflag1.eq.1).and.(distflag2.eq.1)) then
			row = (k*Nsl)+1
		endif
c --- i.e. break out of loop.
	row=row+1
	enddo




	if (inpRow.eq.(N+1)) then
		rowDepth = maxRowDepth+1
        end if
c --- i.e. special detour straight to writing out object buffer if
c --- last time through (final input line)



	if (distflag1.eq.1) then
		dist = 1.1*maxDist
        end if
c --- i.e. also special detour straight past distance comparison 
c --- if current matches previous two obsvns. 



c	write (6,*) inpRow, Xref, Yref, X, Y, int(input(inpRow,1)), dist

	if (dist.le.maxDist) then
c		write (6,*) '	---> obsvn copied across !!' 
		
		objRow = int(input(inpRow,1))
		do col=1,16
			objbuffer(objRow,col)=input(inpRow,col)
		enddo
		numObsv = numObsv+1
c --- i.e. starting position for next round moves down one since
c --- object was copied out and across to the object buffer.
c --- ***N.B: This is the only place where inpStart+1 occurs.

	endif


c =====================================================================
c === Only enters writing-out phase if obsvns thought to be complete:

	if ((rowDepth.gt.maxRowDepth).or.(numObsv.eq.Nsl)) then


	       inpStart=inpStart+1


	       realRow = inpStart*1.0
 	       if ( (inpStart/freq).eq.(int(inpStart/freq)) ) then
  	             write (6,*) '  Passing input line: ',inpStart,
     $                         '    Objects found: ',numObj
   	       endif
c --- i.e. reporting on progress made on working through input catalogue


	if (distflag2.eq.0) then
c --- i.e. only permit newly encountered (Xref,Yref) obsvns past this
c --- point for writing of object buffer to output files.

c	write (6,*) '	** in the assessing/output stage ...'




c --- Copy contents of object buffer to output1 array: --
		do row=1,Nsl
		    do col=1,16
			output1(outRow,col)=objbuffer(row,col)
		    enddo
		    outRow = outRow+1
          	enddo



c --- Create contents for output2 array: ----------------

		do col=1,16
		    meanRow(col) = 0.0
		enddo
		do row=1,Nsl
		    do col=1,16
		      meanRow(col) = meanRow(col)+objbuffer(row,col)
		    enddo
		enddo
		do col=1,16
		    meanRow(col) = meanRow(col)/numObsv
		enddo
c --- calculating mean value of each column in object buffer


	
		bestCol = 0
		bestSN = 0.0
		do col=1,4
		    currSN = meanRow(2*col+3)/meanRow(2*col+4)
		    if (currSN.gt.bestSN) then
			bestSN = currSN
			bestCol = col
		    endif
		enddo
c --- selecting aperture that maximises S/N



		
		if (nextstar.le.Nstars) then
      		dist =(meanRow(3)-stardata(1,nextstar))**2
		dist = dist+(meanRow(4)-stardata(2,nextstar))**2
		dist = sqrt(dist)
		   if (dist.le.starDist) then
			write (6,*) '    Obj ',(numObj+1),' is calibn star '
     $                                         ,nextstar,'.'
c                        write (6,*) '   --> scale factors:'
		   	starflag = 9*(10**Nsl)
			do row=1,Nsl
			   if (objbuffer(row,2*bestCol+3).le.0.0) then
			      write (6,*) '*** star:',nextstar,'    slice:',
     $                                        row,' ***'
			      write (6,*) '*** --> flux = ', 
     $                                        objbuffer(row,2*bestCol+3)
			      write (6,*) '*** set to mean over all',
     $                                          meanRow(2*bestCol+3)
			   objbuffer(row,2*bestCol+3)=meanRow(2*bestCol+3)
                           endif
			   scaleFactors(row,nextstar) = 
     $                       stardata(row+2,nextstar)/
     $			     objbuffer(row,2*bestCol+3)
c			   write (6,*) '       ',scaleFactors(row,nextstar)
			enddo
		   	nextstar = nextstar+1

		   endif
		endif
c --- Checking if current object matches one of the calibration stars;
c --- if it does, determine scaling factors pertaining to it.
c --- Prints and error message if star flux is in error; uses mean value 





		lowSN = 0
		row = 1
		do while (row.le.Nsl)
		  currSN = objbuffer(row,2*bestCol+3)/
     $                       (objbuffer(row,2*bestCol+4)+0.00001)
		  if (currSN.lt.0) then
			lowSN = 1
			row = Nsl
		  endif
		  row = row+1		    
		enddo
c --- checking for neg flux (i.e. S/N<0); flagging this (0=no,1=yes)



c --- Deciding if current object is legit: --------------
c --- Criteria: At least 2 obsvns; SN +ve; SExt flag<8 (neighbs only)

		if ((numObsv.gt.1).and.(lowSN.eq.0).and.
     $			(meanRow(16).lt.8)) then

c		write (6,*) '	==> real object; written to catalogue !!'

		numObj=numObj + 1
		output2(numObj,4) = apDiam(bestCol)
		output2(numObj,13) = meanSigma(bestCol)

c		write (6,*) 'Object',numObj,' aper: ',bestCol

  		do row=1,Nsl
		   output2(numObj,row+13) = objbuffer(row,2*bestCol+3)

		   if (output2(numObj,row+13).le.0.0) then
                      Area = ( (output2(numObj,4)/2)**2 )*3.1415927
c ----------------    numer = 1+Area*( (2*siglevel*gn/meanSigma(bestCol))**2 )
		      numer = 1+Area*( (2*meanSigma(bestCol)*gn/siglevel)**2 )
		      numer = 1+sqrt(numer)
		      denom = 2*gn/(siglevel*siglevel)
		      output2(numObj,row+13) = numer/denom
		   endif
c --- i.e. assigns the threshold flux if none measured by SExtractor


		end do
c --- Copying best aperture data from object buffer to row of output2;
c --- also writing chosen aperture size and sigma for that aperture.
c --- ** Note: ** The sigma are in units of ADU in inApertFile; here we
c --- use same.
c		output2(numObj,13) = meanRow(2*bestCol+4)
c --- above line used to write the mean error in the flux to output.





		Rpix2 = (meanRow(3)-Xcen)**2 + (meanRow(4)-Ycen)**2
		output2(numObj,5) = ((beta/Dpix)**2)*Rpix2*0.5
c --- copying wavelength correction to output2 row
		


		output2(numObj,1) = numObj
		do col=2,3
		   output2(numObj,col) = meanRow(col+1)
		enddo
		output2(numObj,10) = meanRow(15)
c		output2(numObj,6) = meanRow(13)
		output2(numObj,7) = meanRow(14)
c --- i.e. writing ID no., X, Y, star/gal class, ellipticity, FWHM
c --- to output2 (in that order of the program lines above).
c ---
c --- All columns of current row of output2 have been written except
c --- for cols 8-9 (continuum fit coeff; calculated after flux 
c --- normalisation), col 11 (SExtractor flag; calculated below) and
c --- col 12 (emission-line quality factor; after flux normalisation
c --- and continuum fit)



		output2(numObj,11) = 0.0
		do row=1,Nsl
		   placeholder = 10**(row-1)
		   output2(numObj,11) = output2(numObj,11) +
     $			(placeholder*objbuffer(row,16))
		   output2(numObj,11) = output2(numObj,11)+starflag
		   starflag = 0
		enddo
c --- Combining SExtractor flags into a single number; a 9 in leftmost
c --- position means that object is one of the calibration stars

		end if

		end if






c --- Set-up for Next Object: ---------------------------
		inpRow=inpStart

		do row=0,(k*Nsl-2)
			Xprev(k*Nsl-row) = Xprev(k*Nsl-row-1)
			Yprev(k*Nsl-row) = Yprev(k*Nsl-row-1)
		enddo	
		Xprev(1) = Xref
		Yprev(1) = Yref		
		Xref = input(inpRow,3)
		Yref = input(inpRow,4)
		refRow = inpRow
c --- copy X and Y positions and input row number for reference

      	do row=1,Nsl
		objbuffer(row,1)=row
		do col=2,16
			objbuffer(row,col)=0.0
c			write (6,*) row,col,objbuffer(row,col)
		enddo
      	enddo
c --- initialise object buffer array


		objRow = int(input(inpRow,1))
		do col=1,16
		    objbuffer(objRow,col)=input(inpRow,col)
		enddo
c --- copy first entry to object buffer

		numObsv = 1


	endif


	inpRow = inpRow+1
c --- ***N.B: This is the only place where inpRow+1 occurs.

      enddo		

      numObj = numObj-1
 


c =====================================================================



c --------------------------------------------------------------------
c     Determining mean scaling factor over all calibration stars;
c     scaling all entries by this (wavelength by wavelength).
c
c     Then calculating the least-squares fit to the continuum.
c 
c     Then using this to determine emission-line quality factor.
c --------------------------------------------------------------------

	write (6,*) ' ------------------------------------------------'
      write (6,*) '  '
c      write (6,*) 'Mean scaling factors:'

      singleFactor = 0.0
c --- Mean scaling factor over all wavelengths
      do row=1,Nsl
	meanFactors(row) = 0.0
        do col=1,Nstars
	    meanFactors(row) = meanFactors(row)+
     $          scaleFactors(row,col)
	enddo
	meanFactors(row) = meanFactors(row)/Nstars
c	write (6,*) meanFactors(row)	
	singleFactor = singleFactor + meanFactors(row)
      enddo



      singleFactor = singleFactor/Nsl
      write (6,*) 'Single Mean Scaling Factor =',singleFactor,
     $ ' erg/s/cm2/band'


      write (6,*)
      write (6,*) 'Mean scaling factors (after normalisation):'
      do row=1,Nsl
	meanFactors(row) = meanFactors(row)/singleFactor
	write (6,*) meanFactors(row)
      enddo

      write (6,*) ' '
      write (6,*) 'Scaling fluxes and fitting continua ...'
      write (6,*) '  '
c      write (6,*) 'Emission-line candidates above reporting threshold:'
      write (6,*) 'Max. Residual from Cont. Fit for Calibration Stars:'

      do row=1,numObj

	Area = ( (output2(row,4)/2)**2 )*3.1415927
	sig2 = output2(row,13)*output2(row,13)
	S = 0.0
	Sx = 0.0
	Sy = 0.0
	Stt = 0.0
	b = 0.0
	a = 0.0	

	do col=1,Nsl

            output2(row,col+13) = output2(row,col+13)*meanFactors(col)
		

	    t(col) = 0.0
	    corrWaveL(col) = (1-output2(row,5))*centWaveL(col)
	    flux(col) = output2(row,col+13)
	    delF(col) = sqrt( (Area*sig2)+(flux(col)/gn) )

c	    write (6,*) row, corrWaveL(col), flux(col), delF(col)

	    S = S+( 1/(delF(col)*delF(col)) )
	    Sx = Sx+( corrWaveL(col)/(delF(col)*delF(col)) )
	    Sy = Sy+( flux(col)/(delF(col)*delF(col)) )

	enddo 


	do col=1,Nsl
	   t(col) = (corrWaveL(col)-(Sx/S))/delF(col)
	   Stt = Stt+(t(col)*t(col))
	   b = b+( (t(col)*flux(col))/delF(col) )		
	enddo

	b = b/Stt
	a = (Sy-Sx*b)/S
	output2(row,8) = a
	output2(row,9) = b
c --- Continuum-Fitting Routine: corrected set of wavelength values
c --- generated and least-squares fitting done to fluxes.


      enddo


      numReport = 0.0

      do row=1,numObj
	maxChange = 0.0
	uncert = 0.0001
	Area = ( (output2(row,4)/2)**2 )*3.1415927
	sig2 = output2(row,13)*output2(row,13)
	a = output2(row,8)
	b = output2(row,9)


	do col=1,(Nsl-1)

	    f1 = output2(row,col+13)
	    f2 = output2(row,col+1+13)

	    corrWave1 = (1-output2(row,5))*centWaveL(col)
	    corrWave2 = (1-output2(row,5))*centWaveL(col+1)

	    delf1 = sqrt( (Area*sig2)+(f1/gn) )
	    delf2 = sqrt( (Area*sig2)+(f2/gn) )

	    fitflux1 = b*corrWave1+a
	    fitflux2 = b*corrWave2+a

	    contSub1 = f1-fitflux1
	    contSub2 = f2-fitflux2

	    currChange = abs(contSub1-contSub2)


	    if (currChange.gt.maxChange) then
		  maxChange = currChange

		  uncert = relErrLimit*fitflux1
		  if ((delf1/fitflux1).gt.relErrLimit) then
		  	uncert = delf1
		  endif
	    endif

 	enddo

	output2(row,12) = int(maxChange/uncert)
	output2(row,6) = maxChange
c --- calculation of ELQF/writing-out of max. flux change above contin.


	    if (output2(row,11).ge.(8*(10**Nsl))) then
		write (6,*) output2(row,1),' --- Max. Residual 
     $ (after scaling): ', (maxChange/fitflux1*100.0),' percent'
	    endif
c --- i.e. if object row is one of the calibration stars, print
c --- out its maximum flux deviation from the continuum line
c --- (now that scaling has been done) ==> limiting relative error.



	if ((output2(row,6).ge.report).and.(output2(row,12).ge.2)
     $.and.(output2(row,11).eq.0)) then
c	    write (6,*) '    Object: ',row,
c     $                  '        ELQF = ',int(output2(row,12))
            numReport = numReport+1
	endif
c --- ============================================================
c --- SELECTION CRITERIA:
c ---		(1) Max. flux change must be > (report threshold)
c --- 	    and (2) ELQF must be >= 2
c ---       and (3) All SExtractor flags must be 0.
c --- ============================================================


      enddo

      write (6,*) ' '
      write (6,*) 'Number of emission-line objects (> ',report,
     $ '): ',numReport
      write (6,*) 'Total number of Objects =', numObj


c --- Calculating emission-line quality factor and writing this
c --- to output.



      if (nextstar.le.Nstars) then
	write (6,*) '   '
	write (6,*) ' *** NOT ALL CALIBRATION STARS DETECTED ***'
	write (6,*) '     --> ONLY ',(nextstar-1),'/',Nstars,' FOUND.'
	write (6,*) ' ******************************************'
      endif


c --------------------------------------------------------------------
c     Writing output files:
c --------------------------------------------------------------------

      write (6,*) '   '
      write (6,*) 'Writing to files...'

      open(unit=12, file=outFile1, status='unknown')
      do row=1,outRow
		write (12,60) (output1(row,col),col=1,16)
      enddo
c --- i.e. long version of all obsvns, vertical sort style

60    format (f4.0, f7.0, f8.2, f8.2, 4(f12.3, f12.4), f7.3, f6.2,
     $f6.2, f4.0)



      open(unit=13, file=outFile2, status='unknown')
      do row=1,numObj
		write (13,80) (output2(row,col),col=1,(Nsl+13))
      enddo
c --- proper catalogue (full version, all columns)

80    format (f5.0, f8.2, f8.2, f6.1, f12.8, f12.3, f8.2, f12.3, f9.3,
     $f6.2, f14.0, f6.0, f8.4, 7(f12.3))


c 80    format (f5.0, f8.2, f8.2, f6.1, f12.8, f7.3, f8.2, e12.3, 
c     $e12.3, f6.2, f14.0, f6.0, f8.4, 7(e12.3))




      open(unit=10, file=outFile3, status='unknown')
      do row=1,numObj
		write (10,100) output2(row,1), output2(row,2),
     $ output2(row,3),output2(row,4),output2(row,5),output2(row,8),
     $ output2(row,9),output2(row,12),
     $ (output2(row,col),col=14,(Nsl+13))
      enddo
c --- proper catalogue (condensed version, (7+Nsl) columns)

100   format (f5.0, f8.2, f8.2, f6.1, f12.8, f12.3, f9.3, f6.0,
     $  7(f12.4))


c 100   format (f5.0, f8.2, f8.2, f6.1, f12.8, e12.3, e12.3, f6.0,
c     $  7(e12.3))



      write (6,*) 'Done.'
      write (6,*) ' '


      stop
      end


c ===========================================================================