C+
C NAME:
C	centroid()
C PURPOSE:
C	to find the centroid of a star averaged over a transit
C CATEGORY:
C	SMEI Data processing
C CALLING SEQUENCE:
C	call centroid()
C INPUTS:
C
C OUTPUTS:
C
C CALLS:
C       none
C PROCEDURE:
C
C MODIFICATION HISTORY:
C	April, 2003 A. Smith (UCSD/CASS)
C-
	subroutine 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		degperminute
	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(1200,1200)	/1440000*0/
	integer*4	npolarimage(1600,1600)	/2560000*0/
	integer*4	spolarimage(1600,1600)	/2560000*0/
	integer*8	k

	character	infile*38	!not compatable with TMO data
c
c	AARON'S DECLARATIONS
c
	integer*4 iairplane, airplane, idarkcr, ii, jj, dumcnt4, icra
	real*8 rq, dr, thetax, thetaxsq, thetay, cx, cy, cz, crsq
	real*8 rccd, x0, y0, r0, 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)
	

	!!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
	qcam(3,1) = -0.21298450D0
	qcam(3,2) = -0.43870382D0
	qcam(3,3) = 0.86674816D0
	qcam(3,4) = 0.10451884D0
	!!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)

	degperminute = 0.2506845D0

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

	do k=1,19
		if(k.eq.1) iframe=103
		if(k.eq.2) iframe=107
		if(k.eq.3) iframe=115
		if(k.eq.4) iframe=123
		if(k.eq.5) iframe=131
		if(k.eq.6) iframe=135
		if(k.eq.7) iframe=143
		if(k.eq.8) iframe=151
		if(k.eq.9) iframe=155
		if(k.eq.10) iframe=203
		if(k.eq.11) iframe=211
		if(k.eq.12) iframe=219
		if(k.eq.13) iframe=227
		if(k.eq.14) iframe=235
		if(k.eq.15) iframe=243
		if(k.eq.16) iframe=247
		if(k.eq.17) iframe=255
		if(k.eq.18) iframe=303
		if(k.eq.19) iframe=311
		f=dfloat(iframe)

		  ifmin=615
		  if(k.lt.5)then
		    write(infile,'(A31,I3.3,A4)') 'DATA_SMEI/00/c2frm_2003_059_034',iframe,'.nic'
		  else
		    write(infile,'(A31,I3.3,A4)') 'DATA_SMEI/00/c2frm_2003_059_035',iframe,'.nic'
		  endif
		  print*, infile

		call read_frame(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)

C		print*, 'Spacecraft quaternion'
C		print*, qsc(1), qsc(2), qsc(3), qsc(4)

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

		print*, 'Indexing frame...'

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

	enddo

	close(11)
	cra=15.D0*dfloat(nint(cra/15.D0))
	cdec=0.D0
	icra=int(cra)

	print*, "Calling average night"
c	call averagenight(node,nodeid,hits,TOP_KMAX,level,cra,cdec,image, npolarimage, spolarimage)


99	print*, 'Fortran out!'

	return
	end