C+
C NAME:
C	A_mod_mainf
C PURPOSE:
C	This is the main Fortran program to index TMO/SMEI data using utilizing
C	Spatial Index functions from the c++ file A_mod_L12_all.cpp
C CATEGORY:
C CALLING SEQUENCE:
C INPUTS:
C OUTPUTS:
C CALLS:
C SEE ALSO:
C INCLUDE:
C COMMON BLOCKS:
C PROCEDURE:
C	The data frame array, data(600,1280), has its indices reversed
C	from the usual oreder in Fortran because C accesses columns before rows
C	when dealing with arrays.
C MODIFICATION HISTORY:
C	???-2002, Aaron Smith (UCSD/CASS)
C		Original version
C	DEC-2002, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
C		Defined LEVEL and KMAX as parameters. KMAX is determined from LEVEL.
C		KMAX is then used to dimension arrays node, nodeid and hits
C	DEC-2002, Andy Buffington (UCSD/CASS; abuffington@ucsd.edu)
C		Changed bunch of scalars from single to double precision.
C-
	PROGRAM A_mod_mainf

	! Dimensions for node, nodeid, hits above must be at least 2*(4**(level+1))
	! For LEVEL = 12, KMAX = 134217728
	! For LEVEL =  8, KMAX =    524288

	integer*2	  TOP_LEVEL
	parameter	( TOP_LEVEL = 12 )

	integer*8	  TOP_KMAX
	parameter	( TOP_KMAX = 2*(4**(TOP_LEVEL+1)) )

	! The Intel compiler doesn't like this initialization, apparently because
	! TOP_KMAX is an integer*8. Initialization is done explicitly later.

	real*4		node  ( TOP_KMAX )	!/TOP_KMAX*0.0/
	integer*8	nodeid( TOP_KMAX )	!/TOP_KMAX*0  /
	integer*2	hits  ( TOP_KMAX )	!/TOP_KMAX*0  /

	integer*2	level			/TOP_LEVEL/
	integer*8	kmax			/TOP_KMAX/

	real*4		data  (600,1280)
	real*8		euler (3)
	real*8		affine(6)		/6*0.0D0/		! Affine parameters zero for now
	real*8		cra
	real*8		cdec
	real*8		alfa
	real*8		beta
	real*8		gamma
	real*8		camlat
	real*8		degperminute
	real*8		ra0
	real*8		arg
	real*8		arg1
	real*8		dusk

	integer*4	iframe
	integer*4	ifmin                   /2/
	integer*4	ifmax                   /600/
	integer*4	i
	integer*4	j
	integer*4	night			/40/
	integer*4	image(3200,1600)	/5120000*0/
	integer*8	k

	character	infile*18
	character	outfile*10

	print*, 'Fortran IN ...'
	print*, 'Level: ', level, ' (array size: ', kmax,')'
	
	if(night.eq.40)ifmin = 12	!remove 10 minute early start for night 40
	if(night.eq.30)ifmin = 322	!sky reversed
	if(night.eq.30)ifmax = 5808	!sky reversed
	ra0 = 279.425D0

	!cc  adjustments relative to Andy's ra0 for selected stars...night 42 start 19:40
	!    ra0 = ra0 - 0.985626D0           !by-the-clock converts from night 42 to night 41, start 19:40
	!    ra0 = ra0 - 0.985626D0 - 2.5079D0  !by-the-clock converts from night 41 to night 40, start 19:30
	!    note that gap interval setting was 56 ms larger on night 40: see degperminute adjustments below

	if(night.eq.41) ra0 = 279.425D0 + 0.005D0       		!Andy's value with small ad hoc adjustment
	if(night.eq.42) ra0 = 279.425D0 + 0.97D0        		!Aaron's value with ad hoc adjustment (of -0.015626)
	if(night.eq.40) ra0 = 279.425D0 - 0.985626D0 - 2.5079D0 - 0.007D0	!Andy's 41 value, with by-the-clock & adhoc adjustment

	degperminute = 0.2506845D0
	camlat	= 34.9471D0
	if(night.eq.42) degperminute=degperminute*1.0001D0	!Ad hoc adjustment
	if(night.eq.41) degperminute=degperminute*1.00044444D0	!Uses observed end time
	if(night.eq.40) degperminute=degperminute*1.001389D0	!56-millisecond adjustment (+ observed end time)
	if(night.eq.40) camlat = camlat -0.61D0			!Southward looking ad hoc adjustment...
	alfa	= -0.5D0
	beta	= camlat-90.0D0

	do k=1,kmax
	   node(k)=0.0
	   nodeid(k)=0
	   hits(k)=0
	enddo

	do iframe=2,600

		write(infile,'(A8,I2.2,I4.4,A4)') 'DATA/TMO',night,iframe,'.nic'
		print*, infile
		
		arg =dexp(dfloat(1-iframe  )/4.9D0)		!calculate correction timeconstant for dusk
		arg1=dexp(dfloat(iframe-601)/4.9D0)		!and dawn
		if(night.eq.42)  dusk= 16100.D0*arg + 3300.D0*arg1
		if(night.eq.41)  dusk= 24000.D0*arg + 4000.D0*arg1
		if(night.eq.40)  dusk=235000.D0*arg + 1500.D0*arg1

		gamma=(-dfloat(iframe)*degperminute)-ra0
		euler(1)=alfa
		euler(2)=beta
		euler(3)=gamma

		CALL READFILE(infile,data)

		CALL PEDESTAL(data)

		CALL INDEXFRAME(level,node,nodeid,hits,data,night,iframe,dusk,euler,affine)

	enddo

	cra=355.D0
	cdec=40.D0
	CALL AVERAGENIGHT(level,node,nodeid,hits,cra,cdec,image)


	! Write imagecnt(sky atlas) to file

	write (outfile,'(I2.2,A2,I2.2,A4)') night,'_L',level,'.grd'

	open (10,file = outfile, form = 'formatted')
	write(10,10)
10	format('DSAA'/'3198 1598'/'-75 85'/'0 80'/'0 20000')
	do j=2,1599
	  write(10,'(20i6)')(image(i,j),i=2,3199)
	enddo
	close(10)

	!CALL WRITEDATA(outfile,image)

	print*, 'Fortran OUT...'

      END