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	SMEI data analysis
C CALLING SEQUENCE:
C INPUTS:
C	Various parameters are presently specified by data statements at the top
C	This program presently combines data from a single SMEI camera
C OUTPUTS:
C	A data file listing pedestal and dark current values versus frame number
C	A grid file containing the resulting sky map
C CALLS:
C	set_flatfields		!sets the appropriate flatfields for the selected camera
C	read_frame		!reads a frame of data
C	pedestalr		!removes pedestal from empty pixels
C	darkr			!removes dark current from covered pixels
C	cosmicra		!cosmic ray removal
C	lsff			!large scale flatfield
C	ssff			!small scale flatfield
C	rotate8			!euler angle rotation
C	indexframe		!frame indexing routine
C SEE ALSO:
C INCLUDE:
C COMMON BLOCKS:
C PROCEDURE:
C	This program uses the hierarchical triangle coordinate system from Johns Hopkins,
C	as a base upon which the photometric measurements from a single SMEI camera are
C	combined and then averaged to produce a surface-brightness sky map. The "level"
C	parameter governs the granularityof the triangular coordinate frame. The
C	coordinates of the four corners of a pixel (or bin) are transformed down into
C	this frame, and the photometric surface brightness as measured by that pixel is
C	added to all triangles whose centers lie within the rectangle. 
C	This surface brightness is corrected by 1/cosine of the pixel's angle to the camera
C	axis, and by r/r0 the pixel's1 fractional distance to the FOV's rotational center
C	When the data sequence is finished, all triangle response sums are
C	averaged and the resulting sky map binned into an output surface brightness map.
C MODIFICATION HISTORY:
C	2002,2003 Aaron Smith (UCSD/CASS)
C		2002->Original version
C		2003->removed most C++ functions and re-wrote in fortran as subroutines
C		    ->add quaternion capabilities
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	JAN-2003, Andrew Buffington (UCSD/CASS)
C		Integrated standard pedestal and dark current routines, FORTRAN
C		Rendered names compliant with those of B.V. Jackson
C		Incorporated 1/cosine surface-brightness correction for aperture perspective (currently in indexframe routine)
C		Incorporated large and small scale flatfield corrections
C		Incorporated FORTRAN read subroutine
C-
	PROGRAM A_mod_mainf
	implicit none

	! 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)) )

	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		frame      (1272,256)
	real*4		framessff  (1272,256)
	real*4		framelsff  (1272,256)
	real*4		framelsffu (1272,256)
	real*4  	ped			/0./
	real*4		dark			/0./
	real*4		avg        (1272,256)    /0./
	real*4		hroom			/49400./		!'headroom'=new onboard value of 65k ADUs
	real*4		adu2e
	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		euler(3)
	real*8		camlat
	real*8		degperminute
	real*8		ra0
	real*8		arg
	real*8		arg1
	real*8		dusk

	integer*4	mode			/1/			!onboard binning mode
	integer*4	ic			/3/			!Camera #
	integer*4	ifm			/0/			!ROI mask (0=no,1=yes)
	integer*4	id			/1272/			!1280 for TMO
	integer*4	jd			/256/			!600 for TMO
	integer*4	icrcnt			/0/
	integer*4	neighbors		/0/
	integer*4	iframe
	integer*4	iframemask (1272,256)	/325632*1/
	integer*4	ihit(1272,256)
	integer*4	ihist(250)
	integer*4	ifmin                   /193/
	integer*4	ifmax                   /469/
	integer*4	i
	integer*4	j
	integer*4	night			/99/			!night=99 for SMEI data
	integer*4	image1(1200,1200)	/1440000*0/
	integer*4	image2(1200,1200)	/1440000*0/
	integer*4	image3(1200,1200)	/1440000*0/
	integer*4	npolarimage(1600,1600)	/2560000*0/
	integer*4	spolarimage(1600,1600)	/2560000*0/
	integer*8	k

	character	infile*47	!not compatable with TMO data
	character	listfile*27,framename*28, outfile1*21,outfile2*21,outfile3*21, outfile4*17, outfile5*17
c
c	AARON'S DECLARATIONS
c
	integer*4 iairplane, airplane, idarkcr, ii, jj, icra, fcnt
	real*8 rq, dr, thetax, thetaxsq, thetay, cx, cy, cz, crsq
	real*8 rccd, tempx
	real*8 response, arg, arg1, phi, rlat
	real*8 x(4), y(4), ra(4), dec(4)
	real*8 f0, f, xf, yf, zf, rfsq, rf, raf, decf
	real*8 qcam(3,4), qc(4), qsc(4), qandy(4), qdum(4), qfinal(4)
	
C	write(listfile,'(A27)') 'DATA_SMEI/20030227/list.txt'
C	write(listfile,'(A27)') 'DATA_SMEI/20030228/list.txt'
C	write(listfile,'(A27)') 'DATA_SMEI/20030301/list.txt'
	write(listfile,'(A27)') 'DATA_SMEI/20039999/list.txt'

	!!quaternion for cameras from SMEIPoint v1.0
	qcam(1,1) = -0.24387742D0
	qcam(1,2) = 0.46807739D0
	qcam(1,3) = 0.84666176D0
	qcam(1,4) = 0.067758569D0
	qcam(2,1) = -0.27454405D0
	qcam(2,2) = 0.017651734D0
	qcam(2,3) = 0.96032546D0
	qcam(2,4) = 0.045705368D0
C	qcam(3,1) = -0.21298450D0
C	qcam(3,2) = -0.43870382D0
C	qcam(3,3) = 0.86674816D0
C	qcam(3,4) = 0.10451884D0
	!! (1) updated test quaternions
c	qcam(2,1) = -0.265175D0
c	qcam(2,2) = 0.0219207D0
c	qcam(2,3) = 0.963390D0
c	qcam(2,4) = 0.0328890D0
	!! Quaternion from Andy's Euler angles
	qcam(3,1) = 0.216504D0
	qcam(3,2) = 0.432944D0
	qcam(3,3) = -0.868351D0
	qcam(3,4) = -0.107945D0


	!!quaternion to rotate Andy's coordinates to Don's
	qandy(1)=0.5D0
	qandy(2)=0.5D0
	qandy(3)=0.5D0
	qandy(4)=0.5D0
	data ihist /250*0/
	data ihit /81408*0/

	print*, 'Fortran In ...'
c	print*, 'Level: ', level, ' (array size: ', kmax,')'
	adu2e=4.5*65536./hroom
c	call set_flatfields(ic,framessff,framelsff,framelsffu)	!note dimensions should be checked before ff is done...
c	call set_cr_backgnd(mode,ic,ifm,id,jd,avg)

c	if(night.eq.40)ifmin = 13	!remove 10 minute early start for night 40
c	if(night.eq.30)ifmin = 322	!sky reversed
c	if(night.eq.30)ifmax = 580	!sky reversed
c	ra0 = 279.425D0

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

c	if(night.eq.41) ra0 = 279.425D0 + 0.005D0       		!Andy's value with small ad hoc adjustment
c	if(night.eq.42) ra0 = 279.425D0 + 0.97D0        		!Aaron's value with ad hoc adjustment (of -0.015626)
c	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  !!TMO Lattitude
	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.D0

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

	open(12,file=listfile,status='old')
	do k=1,225

	  read(12,12)framename
12	  format(a28)

C	  write(infile,'(A19,A28)') 'DATA_SMEI/20030227/',framename
C	  write(infile,'(A19,A28)') 'DATA_SMEI/20030228/',framename
C	  write(infile,'(A19,A28)') 'DATA_SMEI/20030301/',framename
	  write(infile,'(A19,A28)') 'DATA_SMEI/20039999/',framename
	  print*,k,': ', infile

		if(ic.eq.0)then
		  write(infile,'(A8,I2.2,I4.4,A4)') 'DATA/TMO',night,iframe,'.nic'
		  if(iframe.eq.ifmin)print*, infile
		endif

		call read_framecrx(infile,mode,ic,id,jd,frame,qdum)
		!! input quaternion is of form ai+bj+ck+d, not d+ai+bj+cK -> swap with new quaternion
		qsc(1)=qdum(4)
		qsc(2)=qdum(1)
		qsc(3)=qdum(2)
		qsc(4)=qdum(3)

		!! multiply camera quaternion into spacecraft(frame) quaternion and take the inverse
		!! of camera quaternion because it goes from sky to ccd
		qc(1)=qcam(ic,1)
		qc(2)=-qcam(ic,2)
		qc(3)=-qcam(ic,3)
		qc(4)=-qcam(ic,4)

		call qmultiply(qsc,qc,qdum)
		call qmultiply(qdum,qandy,qfinal)
		!!take inverse of final quaternion: going from ccd to sky, not the other way around

C		qfinal(1)=qfinal(1)
		qfinal(2)=-qfinal(2)
		qfinal(3)=-qfinal(3)
		qfinal(4)=-qfinal(4)

		if(ic.eq.0) then
		  gamma=(-dfloat(iframe)*degperminute)-ra0
		  euler(1)=alfa
		  euler(2)=beta
		  euler(3)=gamma
		endif

C		find ra/dec of center of FOV for a frame
C		  if(iframe.eq.ifmin)then
C		    !!find ra of center of the first frame; use to calculate ra bounds in final image
C		    !!car/cdec calculated from thetax=thetay=0
C		    cra=0.D0
C		    cdec=90.D0
C	            call rotateq(qfinal,cra,cdec) !!rotate to standard system
C		    print*, '			RA/DEC:  ',cra,cdec
c		  endif

c		call pedestalr(mode,ic,ifm,id,jd,iframemask,frame,ped)
c		call darkr(mode,ic,ifm,id,jd,iframemask,frame,dark,adu2e,idarkcr,ihit,ihist)
c		call cosmicra(mode,ic,ifm,id,jd,iframemask,frame,avg,adu2e,icrcnt,ihist,ihit,neighbors)
c		print 1000,iframe,ped,dark,idarkcr,icrcnt,neighbors
1000	format(' Frame ',i6,', Ped = ',f8.2,', Dark = ',f8.2,', cosmic rays ',3i5)
c		write(11,11)iframe,ped,dark,idarkcr,icrcnt,neighbors
11		format(i6,2f10.2,3i7)
c		call lsff(mode,ic,ifm,id,jd,iframemask,framelsff,frame)
c		call ssff(mode,ic,ifm,id,jd,iframemask,framessff,frame)

c		print*, 'Indexing frame...'

		call indexframe(mode,ic,id,jd,level,TOP_KMAX,node,nodeid,hits,night,iframe,frame,affine,qfinal,k)

	enddo  !!end frame do loop

	close(12)

	print*, 'Calling average night'
	call averagenight_allsky(node,nodeid,hits,TOP_KMAX,level,image1,image2,image3, npolarimage, spolarimage)
	print*, 'Night Averaged, Writing Files...'

	!!Write image1(sky atlas) to file
	write (outfile1,'(A21)') 'EqPlot_000_to_120.grd'
	open (10,file = outfile1, form = 'formatted')
	write(10,10)
10	format('DSAA'/'1200 1200'/'0 120'/'-60 60'/'0 2000')
	do j=1,1200
	  write(10,'(20i6)')(image1(i,j),i=1,1200)
	enddo
	close(10)
	
	write (outfile2,'(A21)') 'EqPlot_120_to_240.grd'
	open (13,file = outfile2, form = 'formatted')
	write(13,13)
13	format('DSAA'/'1200 1200'/'120 240'/'-60 60'/'0 2000')
	do j=1,1200
	  write(13,'(20i6)')(image2(i,j),i=1,1200)
	enddo
	close(13)
	
	write (outfile3,'(A21)') 'EqPlot_240_to_360.grd'
	open (14,file = outfile3, form = 'formatted')
	write(14,14)
14	format('DSAA'/'1200 1200'/'240 360'/'-60 60'/'0 2000')
	do j=1,1200
	  write(14,'(20i6)')(image3(i,j),i=1,1200)
	enddo
	close(14)

	write (outfile4,'(I2.2,A2,I2.2,A11)') night,'_L',level,'_npolar.grd'
	open (15,file = outfile4, form = 'formatted')
	write(15,15)
15	format('DSAA'/'1600 1600'/'-40 40'/'-40 40'/'0 2000')
	do j=1,1600
	  write(15,'(20i6)')(npolarimage(i,j),i=1,1600)
	enddo
	close(15)

	write (outfile5,'(I2.2,A2,I2.2,A11)') night,'_L',level,'_spolar.grd'
	open (16,file = outfile5, form = 'formatted')
	write(16,16)
16	format('DSAA'/'1600 1600'/'-40 40'/'-40 40'/'0 2000')
	do j=1,1600
	  write(16,'(20i6)')(spolarimage(i,j),i=1,1600)
	enddo
	close(16)

99	print*, 'Fortran out!'


      END