C+
C NAME:
C	A_mainf_centroid
C PURPOSE:
C	average over transits of a star to get an average of that star
C CATEGORY:
C	SMEI Data processing
C INPUTS:
C
C OUTPUTS:
C
C CALLS:
C       none
C PROCEDURE:
C
C MODIFICATION HISTORY:
C	April, 2003 A. Smith (UCSD/CASS)
C-
	program A_mainf_centroid
	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		ra0
	real*8		arg
	real*8		arg1
	real*8		dusk

	integer*4	mode			/1/			!onboard binning mode
	integer*4	ic			/2/			!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	image(200,200)		/40000*0/

	character	outfile*19,stroutfile*23,listfile*23,framename*28,infile*47
c
c	AARON'S DECLARATIONS
c
	integer*4 idarkcr, ii, jj, icra, meascnt, fcnt, badangle, orbitnum, date, maxk
	integer*8 k, kk, frametest, framemax, starloop, dateloop, orbitloop
	real*8 rq, dr, thetax, thetaxsq, thetay, cx, cy, cz, crsq
	real*8 rccd, x0, y0, r0, tempx, avgthetarot, thetarot
	real*8 response, arg, arg1, phi, rlat, xra(2), ydec(2)
	real*8 x(4), y(4), ra(4), dec(4), avgthetax, thetaxsum
	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)

	!!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
c	qcam(2,1) = -0.27454405D0
c	qcam(2,2) = 0.017651734D0
c	qcam(2,3) = 0.96032546D0
c	qcam(2,4) = 0.045705368D0
	qcam(3,1) = -0.21298450D0
	qcam(3,2) = -0.43870382D0
	qcam(3,3) = 0.86674816D0
	qcam(3,4) = 0.10451884D0
	!! updated test quaternions (e1=-86.10, e2=59.05, e3=-6.75)
	qcam(2,1) = -0.268131D0
	qcam(2,2) = 0.0170381D0
	qcam(2,3) = 0.962053D0
	qcam(2,4) = 0.0476403D0

	!!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 ...'
	adu2e=4.5*65536./hroom
c	call set_flatfields(ic,framessff,framelsff,framelsffu)	!note dimensions should be checked before ff is done...
	call set_cr_backgnd(mode,ic,ifm,id,jd,avg)

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

	do dateloop=3,3
	  if(dateloop.eq.1)then
	    date=0227
	    maxk=3641
	    framemax=14
	  elseif(dateloop.eq.2)then
	    date=0228
	    framemax=14
	  elseif(dateloop.eq.3)then
	    date=0301
	    framemax=14
	    maxk=7307
	  elseif(dateloop.eq.4)then
	    date=0529
	    framemax=30
	    ic=1
	    mode=4
	  endif
	  write(listfile,'(A15,I4.4,A4)') 'DATA_SMEI/list_',date,'.txt'
	  if(date.eq.0529) write(listfile,'(A15,I4.4,A4)') 'DATA_SMEI/cam1_',date,'.txt'
	do orbitloop=1,1
	  orbitnum=orbitloop
	do starloop=1,1

	 open(13,file=listfile,status='old')
	 frametest=14
	 avgthetax=0.D0
	 avgthetarot=0.D0
	 thetarot=0.D0
	 thetaxsum=0.D0
	 fcnt=0

	 do k=1,maxk

 	  read(13,13)framename
13	  format(a28)

	  frametest=frametest+1

	if(starloop.eq.1)then
CCC	  Frame bounds for Aldebaran
	  cra=68.9802D0
	  cdec=16.5093D0
	  if(date.eq.0301)then
	    if(k.eq.6096.and.orbitnum.eq.1) frametest=1
	  elseif(date.eq.0228)then
	    if(k.eq.254.and.orbitnum.eq.1) frametest=1
	    if(k.eq.1164.and.orbitnum.eq.2) frametest=1
	    if(k.eq.2090.and.orbitnum.eq.3) frametest=1
	    if(k.eq.11623.and.orbitnum.eq.4) frametest=1
	  elseif(date.eq.0227)then
	    if(k.eq.355.and.orbitnum.eq.1) frametest=1
	    if(k.eq.1251.and.orbitnum.eq.2) frametest=1
	    if(k.eq.2122.and.orbitnum.eq.3) frametest=1
	    if(k.eq.2997.and.orbitnum.eq.4) frametest=1
	  endif
	  write (outfile,'(A8,I4.4,A1,I2.2,A4)') 'Aldebar_',date,'_',orbitnum,'.grd'

	elseif(starloop.eq.2)then
CCC	Frame bounds for Betelgeuse
	  cra=88.75D0
	  cdec=7.41D0
	  if(date.eq.0227)then
	    if(k.eq.329.and.orbitnum.eq.1) frametest=1
	    if(k.eq.1225.and.orbitnum.eq.2) frametest=1
	    if(k.eq.2096.and.orbitnum.eq.3) frametest=1
	    if(k.eq.2972.and.orbitnum.eq.4) frametest=1
	  endif
	  write (outfile,'(A8,I4.4,A1,I2.2,A4)') 'Betel___',date,'_',orbitnum,'.grd'

	elseif(starloop.eq.3)then
CCC	  Frame bounds for Sirius (HD48915)
	  cra=101.28717D0
	  cdec=-16.71611D0
	  if(date.eq.0301)then
	    if(k.eq.6540.and.orbitnum.eq.1) frametest=1
	  elseif(date.eq.0228)then
	    if(k.eq.155.and.orbitnum.eq.1) frametest=1
	    if(k.eq.1068.and.orbitnum.eq.2) frametest=1
	    if(k.eq.1990.and.orbitnum.eq.3) frametest=1
	    if(k.eq.11524.and.orbitnum.eq.4) frametest=1
	  elseif(date.eq.0227)then
	    if(k.eq.254.and.orbitnum.eq.1) frametest=1
	    if(k.eq.1150.and.orbitnum.eq.2) frametest=1
	    if(k.eq.2020.and.orbitnum.eq.3) frametest=1
	    if(k.eq.2895.and.orbitnum.eq.4) frametest=1
	  endif
	  write (outfile,'(A8,I4.4,A1,I2.2,A4)') 'Sirius__',date,'_',orbitnum,'.grd'

	elseif(starloop.eq.4)then
CCC	  Frame bounds for Vega
	  cra=279.23475D0
	  cdec=38.78361D0
	  if(date.eq.0301)then
	    if(k.eq.6928.and.orbitnum.eq.1) frametest=1
	  elseif(date.eq.0228)then
	    if(k.eq.544.and.orbitnum.eq.1) frametest=1
	    if(k.eq.1463.and.orbitnum.eq.2) frametest=1
	    if(k.eq.11922.and.orbitnum.eq.3) frametest=1
	  elseif(date.eq.0227)then
	    if(k.eq.648.and.orbitnum.eq.1) frametest=1
	    if(k.eq.1538.and.orbitnum.eq.2) frametest=1
	    if(k.eq.2417.and.orbitnum.eq.3) frametest=1
	    if(k.eq.3289.and.orbitnum.eq.4) frametest=1
	  endif
	  write (outfile,'(A8,I4.4,A1,I2.2,A4)') 'Vega____',date,'_',orbitnum,'.grd'

	elseif(starloop.eq.5)then
CCC	  Frame bounds for Alcaid
	  cra=206.8852D0
	  cdec=49.3133D0
	  if(date.eq.0301)then
	    if(k.eq.6919.and.orbitnum.eq.1) frametest=1
	  elseif(date.eq.0228)then
	    if(k.eq.534.and.orbitnum.eq.1) frametest=1
	    if(k.eq.11911.and.orbitnum.eq.2) frametest=1
	  elseif(date.eq.0227)then
	    if(k.eq.639.and.orbitnum.eq.1) frametest=1
	    if(k.eq.1529.and.orbitnum.eq.2) frametest=1
	    if(k.eq.2408.and.orbitnum.eq.3) frametest=1
	    if(k.eq.3280.and.orbitnum.eq.4) frametest=1
	  endif
	  write (outfile,'(A8,I4.4,A1,I2.2,A4)') 'Alcaid__',date,'_',orbitnum,'.grd'

	elseif(starloop.eq.6)then
CCC	Frame bounds for Eta Draconis
	  cra=245.99786D0
	  cdec=61.51421D0
	  if(date.eq.0301)then
	    if(k.eq.6897.and.orbitnum.eq.1) frametest=1
	  elseif(date.eq.0228)then
	    if(k.eq.509.and.orbitnum.eq.1) frametest=1
	    if(k.eq.11887.and.orbitnum.eq.2) frametest=1
	  elseif(date.eq.0227)then
	    if(k.eq.615.and.orbitnum.eq.1) frametest=1
	    if(k.eq.1504.and.orbitnum.eq.2) frametest=1
	    if(k.eq.2383.and.orbitnum.eq.3) frametest=1
	    if(k.eq.3255.and.orbitnum.eq.4) frametest=1
	  endif
	  write (outfile,'(A8,I4.4,A1,I2.2,A4)') 'EtaDrac_',date,'_',orbitnum,'.grd'
	  
	  elseif(starloop.eq.7)then
CCC	Frame bounds for Almach
	  cra=30.975D0
	  cdec=42.330D0
	  if(date.eq.0227)then
	    if(k.eq.439.and.orbitnum.eq.1) frametest=1
	    if(k.eq.1330.and.orbitnum.eq.2) frametest=1
	    if(k.eq.2205.and.orbitnum.eq.3) frametest=1
	    if(k.eq.3078.and.orbitnum.eq.4) frametest=1
	  endif
	  write (outfile,'(A8,I4.4,A1,I2.2,A4)') 'Almach__',date,'_',orbitnum,'.grd'

	elseif(starloop.eq.8)then
CCC	  Frame bounds for HD 20902
	  cra=51.0807D0
	  cdec=49.8611D0
	  if(date.eq.0301)then
	    if(k.eq.6721.and.orbitnum.eq.1) frametest=1
	  elseif(date.eq.0228)then
	    if(k.eq.336.and.orbitnum.eq.1) frametest=1
	    if(k.eq.1250.and.orbitnum.eq.2) frametest=1
	    if(k.eq.11710.and.orbitnum.eq.3) frametest=1
	  elseif(date.eq.0227)then
	    if(k.eq.445.and.orbitnum.eq.1) frametest=1
	    if(k.eq.1334.and.orbitnum.eq.2) frametest=1
	    if(k.eq.2210.and.orbitnum.eq.3) frametest=1
	    if(k.eq.3083.and.orbitnum.eq.4) frametest=1
	  endif
	  write (outfile,'(A8,I4.4,A1,I2.2,A4)') 'HD20902_',date,'_',orbitnum,'.grd'

	elseif(starloop.eq.9)then
C	  print*, 'NSV 6846...'
	  cra=222.682D0
	  cdec=61.75D0
	  if(k.eq.1273) frametest=1
	  write (outfile,'(A8,I4.4,A1,I2.2,A4)') 'NSV6846_',date,'_',orbitnum,'.grd'
	elseif(starloop.eq.10)then
C	  print*, 'HD 124897...'
	  cra=213.752D0
	  cdec=19.183D0
	  if(k.eq.1014) frametest=1
	  write (outfile,'(A8,I4.4,A1,I2.2,A4)') 'HD12489_',date,'_',orbitnum,'.grd'
	elseif(starloop.eq.11)then
C	  print*, 'HD 164058...'
	  cra=269.153D0
	  cdec=51.489D0
	  if(k.eq.1403) frametest=1
	  write (outfile,'(A8,I4.4,A1,I2.2,A4)') 'HD16405_',date,'_',orbitnum,'.grd'
	elseif(starloop.eq.12)then
C	  print*, 'HD 123139...'
	  cra=211.671D0
	  cdec=-22.336D0
	  if(k.eq.756) frametest=1
	  write (outfile,'(A8,I4.4,A1,I2.2,A4)') 'HD12313_',date,'_',orbitnum,'.grd'
	endif

	write (stroutfile,'(A15,A8)') outfile(1:15),'_str.grd'

	if(frametest.le.framemax)then

	      write(infile,'(A14,I4.4,A1,A28)') 'DATA_SMEI/2003',date,'/',framename
	      if(date.eq.0529) write(infile,'(A14,I4.4,A1,A28)') 'DATA_SMEI/cam1',date,'/',framename

	          if(frametest.eq.1) then
		    print*, ' '
		    print*, stroutfile
		  endif
c		  print*, k, infile

		call read_framecrx(infile,mode,ic,id,jd,frame,qdum,meascnt)

cc		if(meascnt.lt.55000)then

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

C		  find ra/dec of center of FOV for a frame
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*, k,' RA/DEC: ',cra,cdec

		  CALL indexframe(mode,ic,id,jd,level,TOP_KMAX,node,nodeid,hits,frame,affine,qfinal,cra,cdec,avgthetarot,thetax,badangle,xra,ydec)
		  if(badangle.eq.0) then
		    fcnt=fcnt+1
		    thetaxsum=thetaxsum+thetax
		  endif

cc	   	endif
	 endif
	 enddo
	 
c	 if(fcnt.le.2) goto 98  !!too many frames were saturated with cosmic rays

	 print*, 'Fcnt: ',fcnt
	 thetarot=avgthetarot/float(fcnt)
	 avgthetax=thetaxsum/float(fcnt)
	 print*, 'Thetax: ',avgthetax
	 print*, 'Thetarot: ',thetarot

	 print*, 'Calling average night...'
c	 call averagenight(node,nodeid,hits,TOP_KMAX,level,cra,cdec,image)
	 call averagenight_rotate(node,nodeid,hits,TOP_KMAX,level,cra,cdec,thetarot,avgthetax,image)

	!!Write imagecnt(sky atlas) to file
	open (10,file = stroutfile, form = 'formatted')
	write(10,10)
10	format('DSAA'/'200 200'/'0 200'/'0 200'/'0 65000')
	do j=1,200
	  write(10,'(20i6)')(image(i,j),i=1,200)
	enddo
	close(10)

98	close(13)
	enddo  !!end starloop
	enddo  !!end orbitloop
	enddo  !!end dateloop

99	print*, 'Fortran out!'

	end
