C+
C NAME:	AaronLite
C	Dec-2004, Andrew Buffington (UCSD/CASS)
C	Quick & Dirty indexing to sky without HTM stuff
C-
	PROGRAM Aaronlite
	implicit none
c
	real*4		frame      (318,64)
	real*4		pattern    (318,64)
	real*4  	ped			/0./
	real*4		dark			/0./
	real*8		t0			/86201.9D0/ 		!time of day in day 365 of 2002 (23h 56m 41.9s)
	real*8		arg
	real*8		affine(6)		/6*0.0D0/		! Affine parameters zero for now
c
	integer*2	coarseskyhits(720,360)	/259200*0/
	integer*4	coarsesky(720,360)	/259200*0/		!coarse all sky map
	integer*4	mode			/4/			!onboard binning mode
	integer*4	ic						!Camera #
	integer*4	ifm			/0/			!ROI mask (0=no,1=yes)
	integer*4	id			/318/			!1280 for TMO
	integer*4	jd			/64/			!600 for TMO
	integer*4	k,kmin(4),kmax(4),mloop,ksect
	integer*4	offset
c
c	AARON'S DECLARATIONS
	character iinfile*44
	character infileprev*44
	character infile*44		!grayburst
c	character infile*47	!44	!Don't forget, infile size also changes in readframe
	character hhead(5)*10
	character framename*25
	character*7 ydh,ydhsav,ydhlist(26),carg
	character*1 one,two
c
	integer*4 i,j,ii,jj,iii,jjj,icorrcnt
	integer*4 iarg
	integer*4 jarg
	integer*4 iffix,pnum,phicknum(26)
	integer*4 ismallmaps(26,11,11,3)
c
	real*4 acorr(1280,301)
	real*4 acorrsum
	real*4 flatcorr(318,64)
	real*4 patdark
	real*4 patcorr
	real*4 cutoff
	real*4 racut(26),decut(26)
	real*4 smallmaps(26,11,11,3),nines
c
	real*8 qcam(3,4)
	real*8 qc(4)
	real*8 qsc(4)
	real*8 qandy(4)
	real*8 qdum(4)
	real*8 qfinal(4)
	real*8 racutd,decutd
c	data infile   /'G:/Sean/Buffer2/Nic288/c2m2_2004_DOY_hhmmss.raw'/
	data infile   /'G:/Grayburst/Orbit#/c2m2_200Y_DOY_hhmmss.nic'/
	data kmin /4*0/
	data smallmaps /9438*0./
	data ismallmaps /9438*0/
	data nines /999.9/
	data ydhlist /'3_079_0','3_084_1','3_085_0','3_187_0','3_206_0','3_337_2','4_006_1','4_083_1','4_094_0',
     &	'4_176_0','4_183_1','4_212_0','4_223_1','4_225_0','4_247_1','4_352_0','4_353_1',
     &	'4_354_0','4_354_1','4_354_1','4_361_1','5_017_1','5_017_1','5_024_1','5_028_0','5_050_1'/
        data phicknum /68,65,64,72,50,41,39,35,34,29,28,26,24,23,19,11,10,9,8,7,3,1,0,78,80,86/
        data one,two/'1','2'/
        data decut / -25.3,-19.1,-11.72,-48.5,-50.77,-39.85,-46.77,-52.35,68.22,-3.58,-40.23,-56.45,-35.06,
     &	-44.72,-25.25,-17.95,71.33,62.83,-33.45,-76.8,73.35,65.93,-70.35,13.03,-34.77,-57.77/
        data racut / 267.88,70.779,292.97,342.9,308.45,120.625,178.06,208.47,115.225,195.04,311.94,238.32,
     &  358.53,246.52,270.84,164.79,24.775,6.11,167.675,343.96,79.54,358.425,243.38,192.88,219.59,81.275/
	ic=2
c	
	!! 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
	!! refined camera 1 quaternions [-28.5,70.4,-21.5]
	qcam(1,1) = -0.237689D0
	qcam(1,2) = 0.467706D0
	qcam(1,3) = 0.848236D0
	qcam(1,4) = 0.0724612D0
	!! test quaternion for camera 3 [-143.3, 63.43, -0.42]
	qcam(3,1) = 0.203781D0
	qcam(3,2) = 0.437286D0
	qcam(3,3) = -0.869462D0
	qcam(3,4) = -0.106259D0
	!!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
	if(ic.eq.1)offset=17
	if(ic.eq.2)offset=11
c
	open (13,file='c2pat_2003_144_4x4.grd',readonly)	!DOY144
	print *,'reading in c2pat_2003_144_4x4.grd'		!DOY144
c
	read(13,13)hhead
13	format(a10)
	do j=1,64
	read(13,15)(pattern(i,j),i=1,318)
15	format(20f9.3)
	enddo
	read (13,16)patdark
16	format(f10.4)
	print *,'patdark = ',patdark	
	close(13)
c	
	open (10,file='D:\Ccd\birmingham\flatfield\newfltgrd2.GRD',readonly)
	print *,'reading in D:\Ccd\birmingham\flatfield\newfltgrd2.GRD'
	read(10,1)hhead
1	format(a10)
	read(10,3)((acorr(i,j),i=1,1280),j=1,301)
3	format(20f10.5)
	iffix=0
	do i=5,id-5
	do j=5,jd-5
	    iarg=2+mode*(i-1)
	    jarg=offset+mode*(j-1)
	    if(acorr(iarg,jarg).eq.0.)then
	        acorrsum=0.
	        icorrcnt=0
		do ii=-1,1
		do jj=-1,1
		    if(ii.ne.0.or.jj.ne.0)then
		        if(acorr(iarg+mode*ii,jarg+mode*jj).ne.0.)then
			    acorrsum=acorrsum + acorr(iarg+mode*ii,jarg+mode*jj)
			    icorrcnt=icorrcnt+1
		        endif
		    endif
		enddo
		enddo
		if(icorrcnt.ge.6)then
		    acorrsum=acorrsum/float(icorrcnt)
		    iffix=iffix+1
		    do iii=0,mode-1
		    do jjj=0,mode-1
		        acorr(iarg+iii,jarg+jjj)=acorrsum
		    enddo
		    enddo
		endif
	    endif
	enddo
	enddo
	print *,iffix,'large scale flatfield entries filled in with neighbors average.'
	do i=1,id
	do j=1,jd
	    iarg=2+mode*(i-1)
	    jarg=offset+mode*(j-1)
	    flatcorr(i,j)=acorr(iarg,jarg)
	enddo
	enddo
	close(10)
	open (11,file='G:\grayburst\temp.txt')
	open (12,file='G:\grayburst\framelist2.txt')
c
	kmax(1)= 1309	!grayburst
	kmax(2)= 1299	!grayburst
	kmax(3)= 1139	!grayburst
	kmax(4)= 0	!grayburst
c	kmax(1)= 1500	!21
c	kmax(2)= 3000	!67	!20
c	kmax(3)= 4500
c	kmax(4)= 6000
c	kmin(1)=  0
c	kmin(2)= 1500
c	kmin(3)= 3000
c	kmin(4)= 4500
c	kmin(1)= 360	!PPH#35
c	kmin(2)= 360	!"
c	kmax(1)= 442	!"
c	kmax(2)= 442	!"
c	kmin(2)=416	!just 130447
c	kmax(2)=417	!"
c
c***************************Begin main loops
c
	do mloop = 1,3	
c
c	Fix ydhlist entries where later orbits advance 10's of hours column: a big pain!
c	
	if(mloop.eq.2)then
		carg=ydhlist(1)
		write(carg(7:7),'(a1)')one
		ydhlist(1)=carg
		carg=ydhlist(3)
		write(carg(7:7),'(a1)')one
		ydhlist(3)=carg
		carg=ydhlist(5)
		write(carg(7:7),'(a1)')one
		ydhlist(5)=carg
		carg=ydhlist(21)
		write(carg(7:7),'(a1)')two
		ydhlist(21)=carg
		carg=ydhlist(26)
		write(carg(7:7),'(a1)')two
		ydhlist(26)=carg
	endif
	if(mloop.eq.3)then
		carg=ydhlist(1)
		write(carg(7:7),'(a1)')one
		ydhlist(1)=carg
		carg=ydhlist(3)
		write(carg(7:7),'(a1)')one
		ydhlist(3)=carg
		carg=ydhlist(5)
		write(carg(7:7),'(a1)')one
		ydhlist(5)=carg
		carg=ydhlist(7)
		write(carg(7:7),'(a1)')two
		ydhlist(7)=carg
		carg=ydhlist(20)
		write(carg(7:7),'(a1)')two
		ydhlist(20)=carg
		carg=ydhlist(21)
		write(carg(7:7),'(a1)')two
		ydhlist(21)=carg
		carg=ydhlist(23)
		write(carg(7:7),'(a1)')two
		ydhlist(23)=carg
		carg=ydhlist(26)
		write(carg(7:7),'(a1)')two
		ydhlist(26)=carg
	endif
c	
	  ksect=1	
	  ydhsav='0000000'
	  pnum=1
c	  pnum=7	!PPH# 35
	  write(11,1001)
1001	format(/)	  
	  do i=1,720
	  do j=1,360
	      coarsesky(i,j)=0
	      coarseskyhits(i,j)=0
	  enddo
	  enddo
	if(mloop.eq.1)open(17,file='G:\grayburst\cam2_grayburst1_CRX.txt',readonly)
	if(mloop.eq.2)open(17,file='G:\grayburst\cam2_grayburst2_CRX.txt',readonly)
	if(mloop.eq.3)open(17,file='G:\grayburst\cam2_grayburst3_CRX.txt',readonly)
c	if(mloop.eq.4)open(17,file='G:\grayburst\cam2_grayburst4_CRX.txt',readonly)
c	open(17,file='G:\Sean\eclipse288c2_CRX.txt',readonly)
c	
c	wind off unwanted frames at the beginning
c	  
	  if(kmin(mloop).gt.0)then
	  do k=1,kmin(mloop)
	  read(17,170)iinfile
	  enddo
	  endif
c
	  do k=kmin(mloop)+1,kmax(mloop)
c		read(17,170)iinfile,ped,dark
		read(17,170)infile,ped,dark	!grayburst
170		format(a44,2x,2f10.3)
		if(mloop.eq.1.and.k.ge. 555.and.k.le. 565)go to 99	!kill bad CRs frames
		if(mloop.eq.2.and.k.ge.1116.and.k.le.1122)go to 99	!kill bad CRs frames
c	    write(infile(1:44),'(a44)')iinfile	!other than grayburst
c	    print *,infile
c		
c	Load up list of time frames of GRBs
c
	        write(ydh,'(a7)')infile(29:35)
		if(mloop.eq.3.and.k.eq.649)then
			carg=ydhlist(16)
			write(carg(7:7),'(a1)')one
			ydhlist(16)=carg
			ydhsav=ydh
		endif
	if(ydh.ne.ydhsav.or.k.eq.kmax(mloop).or.(pnum.eq.19.and.k-ksect.eq.50).or.(pnum.eq.22.and.k-ksect.eq.50))then
	        	if(k.ne.1)then
	        	write(11,1000)k,ydhsav,k-ksect,pnum,phicknum(pnum)
c	        	print 1000,k,ydhsav,k-ksect,pnum,phicknum(pnum)
	        	write(12,1002)phicknum(pnum),iinfile,infileprev
	        	print 1002,phicknum(pnum),iinfile,infileprev
1000	format(' frame ',i4,', 200',a7,'HMMSS, ',i3,' frames in window for #',i3,', P.Hick # = ',i3)
1002	format(' PPH #',i3,', frames ',a44,' to ',a44)
			endif
	        	ydhsav=ydh
	        	ksect=k
	        	iinfile=infile
	        	if(k.ne.kmax(mloop))then
	        		if(k.ne.1)pnum=pnum+1
	        		if(ydh.ne.ydhlist(pnum))print *,'advancing past #',pnum,', at k = ',k
	        		if(ydh.ne.ydhlist(pnum))pnum=pnum+1
c	        		print *,'Beginning ',ydh
	        	endif
	        endif
	    racutd=dfloat(racut(pnum))
	    decutd=dfloat(decut(pnum))
	    cutoff=0.
c	    
c	    if(cutoff.eq.0.)go to 99	!quick bail
c
	    if(dark.le.-9999.)print*, 'BAD DARK CURRENT VALUE IN FRAME ',framename
	    if(ped.le.-9999.)print*, 'BAD PEDESTAL VALUE IN FRAME ',framename
c	    
	    if(dark.gt.-9999..and.ped.gt.-9999..and.cutoff.le.120.)then	!bail on funny ped, dark
c
		if(mod(k,100).eq.0)print *,k,' ',infile
	      call read_frame_q(infile,mode,ic,id,jd,frame,qdum)
	      patcorr=dark/patdark
c
	      do i=1,id
	      do j=1,jd
	        arg=frame(i,j)
	        arg = arg/0.75		                           !Headroom adjustment (may not comply with engineering data and cam3)
		arg=arg-ped                                        !Remove pedestal
c		if(i.gt.4.and.i.lt.318)arg=arg-pattern(i,j)*patcorr   !Remove dark current and pattern
		if(i.ge.5.and.i.le.317)arg=arg-0.66*dark*(1.+float(j)/64.)	!linear slope removal, naive...
		if(flatcorr(i,j).ne.0.)arg=arg/flatcorr(i,j)	   !Flat Field Correction
		if(frame(i,j).ne.0.)frame(i,j)=arg
	      enddo
	      enddo
c
c		Make SURFER frame
c
c		open (12, file='G:\grayburst\c2m2_2004_083_130447.grd')
c		write (12,12)
12	format('DSAA'/'318 64'/'0 1272'/'0 256'/'0 1000')
c		do jj=1,64
c		j=65-jj
c		write(12,120)(frame(i,j),i=1,318)
120	format(20f9.1)		
c		enddo
c		close(12)		
c
c		input quaternion is of form ai+bj+ck+d, not d+ai+bj+cK -> swap with new quaternion
c
		  qsc(1)=qdum(4)
		  qsc(2)=qdum(1)
		  qsc(3)=qdum(2)
		  qsc(4)=qdum(3)
c
		  !! 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)
c
		  call qmultiply(qsc,qc,qdum)
		  call qmultiply(qdum,qandy,qfinal)
c
c		  !!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
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*, infile,'	[',cra,',',cdec,']'
c
	  call indexframe(mode,ic,id,jd,frame,affine,qfinal,patdark,coarsesky,coarseskyhits,racutd,decutd,pnum,
     &    smallmaps,ismallmaps,mloop)
c
	  endif
	  infileprev=infile
99	  continue	  
	  enddo  	!end k frame loop
	  close(17)
c
	  do i=1,720
	    do j=1,360
	      if(coarseskyhits(i,j).ne.0)then
	        coarsesky(i,j)=coarsesky(i,j)/coarseskyhits(i,j)
	      else
	        coarsesky(i,j)=0
	      endif
	    enddo
	  enddo
c	  if(mloop.eq.1)open (20,file = 'g:\grayburst\cam2gray_orbit1aa.grd', form = 'formatted')
c	  if(mloop.eq.2)open (20,file = 'g:\grayburst\cam2gray_orbit2aaz.grd', form = 'formatted')
c	  if(mloop.eq.3)open (20,file = 'g:\grayburst\cam2gray_orbit3a.grd', form = 'formatted')
c	  if(mloop.eq.4)open (20,file = 'g:\grayburst\cam2gray_orbit4.grd', form = 'formatted')
c	  if(mloop.eq.1)open (20,file = 'g:\Sean\eclipse\cam2_orbit1.grd', form = 'formatted')
c	  if(mloop.eq.2)open (20,file = 'g:\Sean\eclipse\cam2_orbit2.grd', form = 'formatted')
c	  if(mloop.eq.3)open (20,file = 'g:\Sean\eclipse\cam2_orbit3.grd', form = 'formatted')
c	  if(mloop.eq.4)open (20,file = 'g:\Sean\eclipse\cam2_orbit4.grd', form = 'formatted')
c	  write(20,20)
20	  format('DSAA'/'720 360'/'0 360'/'-90 90'/'0 1400')
	  do j=1,360
c	    write(20,'(20i8)')(coarsesky(i,j),i=1,720)
	  enddo
c	  close(20)
	  enddo		!end mloop on 4 data sections
	  do ii=1,26
	  do jj=1,3
	  do i=1,11
	  do j=1,11
	  if(ismallmaps(ii,i,j,jj).ne.0)smallmaps(ii,i,j,jj)=smallmaps(ii,i,j,jj)/float(ismallmaps(ii,i,j,jj))
	  if(smallmaps(ii,i,j,jj).gt.999.)smallmaps(ii,i,j,jj)=999.
	  if(smallmaps(ii,i,j,jj).lt.-99.)smallmaps(ii,i,j,jj)=-99.
	  enddo
	  enddo
	  enddo
	  enddo
c	  open (20,file='g:\grayburst\cam2smallmaps.grd')
c	  write(20,21)
21	  format('DSAA'/'312 35'/'0.5 26.5'/'0.5 3.5'/'0 700')
	  do jj=1,3
	  do j=1,11
c	  write(20,22)((smallmaps(ii,i,j,jj),i=1,11),nines,ii=1,26)
22	  format(20f6.1)
	  enddo
c	  if(jj.lt.3)write(20,22)(nines,i=1,312)
	  enddo
98	continue
	print*, 'All Orbits Finished!!!'
	END
C+
C NAME:
C	read_frame(infile,mode,ic,id,jd,frame)
C PURPOSE:
C	To read a frame of SMEI data and convert the unsigned short integers to floating 
C CATEGORY:
C	SMEI Data processing
C CALLING SEQUENCE:
C	call read_frame(infile,mode,ic,id,jd,frame)
C INPUTS:
C	infile		   character	Input file name (currently 18 characters)
C***					This will have to change...
C	mode		   integer	1 - "Engineering" frame records, no onboard binning
C					2 - "Hi-res" frame records, binned 2 x 2 pixels
C					4 - "Normal Data" frame records, binned 4 x 4 pixels
C	ic		   integer	Camera 0,1,2,3 or 4
C	id		   integer	data frame max x-value (along rows)
C	jd		   integer	data frame max y-value (down columns)
C OUTPUTS:
C	frame(id,jd)   	   real		data frame 
C FUNCTIONS/SUBROUTINES:
C       none
C PROCEDURE:
C	Wind off header, read in one data frame, convert to short unsigned integer, float to real
C	Note: ***This procedure probably needs significant alteration for the SMEI data
C	coming to us on the weekly CD's.***
C	In particular, the FITS header probably contains desirable information which will
C	want to be passed out through this subroutine's argument list. Also, the unsigned
C	short integer format of our *.nic files is probably not Hanscom's convention.
C	yoffset is the number of rows skipped by EGSE before starting its 256-row readout rectangle
C MODIFICATION HISTORY:
C	Jan, 2003 A. Buffington (UCSD/CASS)
C	Feb, 2003 A. Smith (UCSD/CASS)
C	Dec. 2004 cut down version for AaronLite program
C-
	subroutine read_frame_q(infile,mode,ic,id,jd,frame,q)
	implicit none
	integer*4 i,j,id,jd,itemp,ic,mode
	integer*1 trail1(72),trail2(152)
	real*8  q(4)
	character*1 a3(3)
	integer*2 idata(318,64),head(3)		!can't see why id,jd don't work here...
	real*4 frame(id,jd),yoffset(5)
	character infile*44			!not compatable with TMO data, or much of anything else
c	character infile*47			!not compatable with TMO data, or much of anything else
	data yoffset /0.,65.,59.,62.,55./
c
c	Note: this routine does not yet read flight-data files, only TMO (RAL format) data.
c	Besides winding off the proper header, this will also need to include the 
c	yoffset provision above -- requires care since 2 x 2 and 4 x 4 binning not lined up
c	well with the above offsets. Have a nice day!
c
	if(ic.gt.0)then
	   if(mode.ne.1.and.mode.ne.2.and.mode.ne.4)then
	     print *,' Bummer!!, bad mode = ',mode
	     return
	   endif
	   open(10,file=infile,form='binary',readonly)
           read(10)head,((idata(i,j),i=1,id),j=1,jd),a3,trail1,q,trail2
           close(10)
c
c	   convert from short unsigned integer to real*4
c
  	   do j=1,jd
	      do i=1,id
	        itemp=int4(idata(i,j))
	        if(itemp.lt.0)itemp=itemp+65536
	        frame(i,j)=float(itemp)
c
c	Now fix covered pixels
c
		if(mode.eq.1.and.(i.eq.21.or.i.eq.1264))frame(i,j)=0.0
		if(mode.eq.2.and.(i.eq.11.or.i.eq.632))frame(i,j)=2.*frame(i,j)
		if(mode.eq.4.and.(i.eq.6.or.i.eq.316))frame(i,j)=1.333*frame(i,j)
	      enddo
	   enddo

	else
           open(10,file=infile,form='binary',readonly)
           read(10)head,((idata(i,j),i=1,1280),j=1,600)
           close(10)
c
c	convert from short unsigned integer to real*4
c
  	   do j=1,600
	      do i=1,1280
	        itemp=int4(idata(i,j))
	        if(itemp.lt.0)itemp=itemp+65536
	        frame(i,j)=float(itemp)
	      enddo
	   enddo
	endif
c
	return
	end
C+
C NAME:
C	qmultiply
C PURPOSE:
C	multiply two quaternions together
C CATEGORY:
C	Math
C CALLING SEQUENCE:
C	call qmultiply(Q1, Q2, Q)
C INPUTS:
C	Q1(4)	real	quaternion on the left of the multiplication
C	Q2(4)	real	quaternion on the right of the multiplication
C OUTPUTS:
C	Q(4)	real	final quaternion
C MODIFICATION HISTORY:
C	2003 March, Aaron (UCSD/CASS)
C-
	subroutine qmultiply(q1,q2,q)
	real*8 q1(4), q2(4), q(4)
	q(1) = q1(1)*q2(1) - q1(2)*q2(2) - q1(3)*q2(3) - q1(4)*q2(4)
	q(2) = q1(2)*q2(1) + q1(1)*q2(2) - q1(4)*q2(3) + q1(3)*q2(4)
	q(3) = q1(3)*q2(1) + q1(4)*q2(2) + q1(1)*q2(3) - q1(2)*q2(4)
	q(4) = q1(4)*q2(1) - q1(3)*q2(2) + q1(2)*q2(3) + q1(1)*q2(4)
	return
	end
C+
C NAME:
C	indexframe_q_test2
C PURPOSE:
C	To index an entire frame of SMEI data using a generalized NxN binning scheme
C CATEGORY:
C	Data processing
C CALLING SEQUENCE:
C	call indexframe(mode,ic,id,jd,level,TOP_KMAX,node,nodeid,hits,frame,affine,q,badcnt,totcnt,BadPix,
C     &	cutout,cutcnt,cutcounter,pattotal,patdark,coarsesky,coarseskyhits,skytime,orbitfraction)
C INPUTS:
C	mode			camera binning style; mode=N -> NxN binning
C	ic			camera number(0-4); 0 -> TMO camera
C	id			data frame max x-value (along rows)
C	jd			data frame max y-value (down columns)
C	level			spatial index level to be used
C	TOP_KMAX		dimension of the three triangle arrays(node,nodeid, hits)
C	node(TOP_KMAX)		array storing summed response for a given spatial triangle
C	nodeid(TOP_KMAX)	array storing triangle name in spatial index format(uint64 in C++)
C	hits(TOP_KMAX)		array storing amount of frames contibuting to a triangle's response
C	frame(id,jd)		2-D array storing the pixelated ccd picture
C	affine(6)		array storing parameters for affine transformation
C	q(4)			quaternion for coordinate rotation
C OUTPUTS:
C	node, nodeid, and hits are the only function arguments that this subroutine changes
C	and returns top the main program
C FUNCTIONS/SUBROUTINES:
C	rotate8(alfa, beta, gamma, phi,rlat)
C	rotateq(q, phi, rlat)
C	INDEXPIXEL(level,node,nodeid,hits,response,ra,dec,tcount,ttotal)
C PROCEDURE:
C	find the 4 vertices of a pixel, transform those ccd coordinates to sky coordinates,
C	and send the 4 ra/dec vertices to the C++ INDEXPIXEL routine.
C MODIFICATION HISTORY:
C	Jan-May, 2003 Aaron Smith (UCSD)
C	  note: not compliant with TMO data anymore
C-
c	subroutine indexframe(mode,ic,id,jd,level,TOP_KMAX,node,nodeid,hits,frame,affine,q,badcnt,totcnt,BadPix,
c     &	cutout,cutcnt,cutcounter,pattotal,patdark,coarsesky,coarseskyhits,skytime,orbitfraction)
	subroutine indexframe(mode,ic,id,jd,frame,affine,q,patdark,coarsesky,coarseskyhits,racutd,decutd,pnum,
     &  smallmaps,ismallmaps,mloop)
	implicit none
c	integer*8 TOP_KMAX, kk
c	integer*8 nodeid(TOP_KMAX),tcount, ttotal, totcnt, badcnt
c	integer*2 hits(TOP_KMAX), coarseskyhits(720,360),iflag,mflag
c	real*4 node(TOP_KMAX),patdark,patarg,z(3,3),z0,zcut,sumzsqr,skytime(720,360),orbitfraction
	integer*4 xmin, xmax, ymin, ymax, mode, ic, xoff, id, jd, ira, jdec
	real*4 patdark,patarg,z(3,3),z0,zcut,sumzsqr	!,skytime(720,360),orbitfraction
c	real*4 frame(id,jd), yoffset(5), xoffset(5), minr,cutcnt(720,360),cutcounter(720,360),cutout(720,360),pattotal(636,128)
	real*4 frame(id,jd), yoffset(5), xoffset(5), pattotal(636,128)
	real*4 smallmaps(26,11,11,3)
	real*8 x(4), y(4), ra(4), dec(4), euler(4), q(4), affine(6), aq, bq, cq, racutd, decutd, cutrad, delra
	real*8 rq, dr, thetax, thetaxsq, thetay, cx, cy, cz, crsq, czsum, drmin, drmax, thetsum
	real*8 x0(5), y0(5),x00(5), y00(5), yref(5), r0, tempx, phi, rlat, rccd, arg, arg1, response, i, j
c	integer*4 ii, iii,jjj,jj, BadPix(id,jd),coarsesky(720,360)
	integer*4 ii, iii,jjj,jj, coarsesky(720,360)
	integer*4 ismallmaps(26,11,11,3),pnum,mloop
c	integer*4 itemp, jtemp, itempmin, itempmax, jtempmin, jtempmax, rcount
	integer*2 k, level, iflag, mflag, coarseskyhits(720,360)
	integer*8 badcnt, tcount, ttotal, totcnt

C	numbers from Andy's memo
	data yoffset /0.,65.,59.,61.,55./
	data xoffset /0.,1.,1.,1.,1./
	data x00 /630.D0,634.D0,628.7D0,637.7D0,633.D0/
	data y00 /1291.38D0,1306.1D0,1297.1D0,1304.9D0,1294.3D0/
	if(ic.eq.1) r0=1193.5D0
	if(ic.eq.2) r0=1188.5D0
	if(ic.eq.3) r0=1198.5D0
	totcnt=0
	badcnt=0
	iflag=0
	mflag=0

	x0(ic+1)=x00(ic+1)-xoffset(ic+1)
	y0(ic+1)=y00(ic+1)-yoffset(ic+1)

	if(mode.eq.1)then
	  ymin=1
	  ymax=jd
	  xmin=22
	  xmax=1263 !!(1264-1) to remove the last half covered pixel
	  xoff=0
	  zcut=111.	!500 electron cut
	  if(ic.eq.1.or.ic.eq.2)then
	    drmin=-27.D0
	    drmax=22.D0		!was 32.D0: should be 22.D0 according to 11 July 2001 memo, page 5
	  elseif(ic.eq.3)then
	    drmin=-26.D0
	    drmax=20.D0
	  endif
	else if(mode.eq.2)then
	  ymin=1
	  ymax=jd
	  xmin=11
	  xmax=632
	  xoff=0
	  zcut=27.76	!500 electron cut
	  if(ic.eq.1.or.ic.eq.2)then
	    drmin=-27.D0
	    drmax=22.D0
	  elseif(ic.eq.3)then
	    drmin=-26.D0
	    drmax=20.D0
	  endif
	else if(mode.eq.4)then
	  ymin=1
	  ymax=jd
	  xmin=6	!**************do entire frame
c	  xmin=161	!**************do only half-frame away from Sun
	  xmax=316
	  xoff=0
	  zcut=6.94	!500 electron cut
	  if(ic.eq.1.or.ic.eq.2)then
	    drmin=-26.D0
	    drmax=21.D0		!was 31.D0: see note above
	  elseif(ic.eq.3)then
	    drmin=-25.D0
	    drmax=19.D0
	  endif
	endif

	do jj=ymin,ymax
	do ii=xmin,xmax
	  iflag=0
	  if(ic.eq.3)then
	     patarg=pattotal(ii,jj)-patdark*0.4*(1.+0.01*float(jj))	!remove dark current slope from hot/flipper pixel flag
	     if(patarg.gt.50.) iflag=1
	  endif
	  mflag=0
	  response = frame(ii,jj)
	  if(response.gt.0.)then
	     if(ic.eq.3.and.ii.gt.1.and.ii.lt.xmax.and.jj.gt.1.and.jj.lt.ymax)then	!? apply to cam 1 & 2?
	     do iii=1,3
	     do jjj=1,3
	     	z(iii,jjj)=frame(ii+iii-2,jj+jjj-2)
	     enddo
	     enddo
	     call plane(z,z0)
	     if(abs(z(2,2)).gt.zcut)then
	     	sumzsqr=0.
	     	do iii=1,3
	     	do jjj=1,3
	     	if(iii.ne.2.or.jjj.ne.2)sumzsqr=sumzsqr+z(iii,jjj)**2
	     	enddo
	     	enddo
	     	sumzsqr=sumzsqr/5.
c		if((z(2,2)*z(2,2))/sumzsqr.gt.5.)then
c		print 1000,response,sumzsqr,z(2,2),z0,z(2,2)**2/sumzsqr,ii,jj
c		endif
1000	format('response = ',f8.2,'  sumzsqr = ',f10.2,'  z(2,2) = ',f7.2,'  z0 = ',f7.2,'  z(2,2)**2/sumzsqr = ',f10.2,
     &		'  (i,j) = (',i3,',',i3,')')
		if((z(2,2)*z(2,2))/sumzsqr.gt.25.)mflag=1
		if(iflag.eq.1.and.z(2,2).lt.0.)mflag=1
	     endif
	     endif
CC	     find unbinned ccd coordinates
CC	     -1.0 because ccd pixels should start 1 less to make compliant with old C version
	     i = dfloat( ii*mode ) - dfloat(mode-1)/2.D0 - 1.D0 + xoff
	     j = dfloat( jj*mode ) - dfloat(mode-1)/2.D0 - 1.D0

CC	     find 4 vertex coordinates for pixel
	     x(1) = i - dfloat(mode)/2.D0
	     if(mode.gt.1.and.ii.eq.xmin)x(1) = x(1) + 1.D0	!fix geometry of stamped pixel
	     y(1) = j - dfloat(mode)/2.D0
	     x(2) = x(1)
	     y(2) = j + dfloat(mode)/2.D0
	     y(3) = y(2)
	     x(3) = i + dfloat(mode)/2.D0
	     if(mode.gt.1.and.ii.eq.xmax)x(3) = x(3) - 1.D0	!fix geometry of stamped pixel
	     x(4) = x(3)
	     y(4) = y(1)

CC	     this rccd/dr are only preliminary quantities for determining
CC	     if the pixel is in the FOV
	     rccd = dsqrt( (i-x0(ic+1))*(i-x0(ic+1)) + (j-y0(ic+1))*(j-y0(ic+1)) )
	     dr = rccd - r0

	     !!this dr test is for cameras 1, 2, and 3
	     if( dr.lt.drmax.and.dr.gt.drmin)then
	      if( (dabs( (i-x0(ic+1))/(j-y0(ic+1)) ) ).lt.0.57735D0 )then !!less than tan 30

	       response=response*(rccd/r0)	!first order radial binsize correction
c
c	       loop for each vertex of a pixel
c
	       czsum=0.D0
	       thetsum=0.D0
c	       do k=1,4
	       do k=1,1		!*******************Note doing only one vertex, see consequences below
c
	         !!affine transformation to std. sky coordinates
		 tempx = affine(1) + (1.D0+affine(2))*x(k) + affine(3)*y(k);
		 y(k) = affine(4) + affine(5)*x(k) + (1.D0+affine(6))*y(k);
		 x(k) = tempx;
	         rccd = dsqrt((x(k) - x0(ic+1))*(x(k) - x0(ic+1)) + (y(k) - y0(ic+1))*(y(k) - y0(ic+1)));
	         rq = (rccd - r0)
		 thetax = dasind( (x(k) - x0(ic+1)) / rccd );

		 !! quadratic coefficients from Andy's July 2001 memo
		 aq = 0.1664D0 - (0.00001*thetax*thetax)
		 bq = (0.000036D0*thetax*thetax) - 18.0226D0
		 cq = (0.000936D0*thetax*thetax) - rq
		 thetay = ( -bq-sqrt((bq*bq)-(4.D0*aq*cq)) ) / ( 2.D0*aq )
		 thetsum=thetsum+thetay
C	         thetay = rq / -18.0226D0

	         cx = dsind(thetax)
	         cy = dsind(thetay)
	         crsq = (cx*cx)+(cy*cy)
	         if(crsq.gt.0.D0) then
cc		   cx=-cx
cc		   cy=-cy
	           phi = datan2d(cy,cx)
	         else
	           phi = 0.D0
	         endif
	         cz = dsqrt(1.D0-(crsq))
		 czsum = czsum + cz
	         rlat = dasind(cz)

		 call rotateq(q, phi, rlat) !!rotate to standard system

		 ra(k)=phi
		 dec(k)=rlat
		 cutrad=dabs((ra(k)-racutd)*dcosd(dec(k)))**2 + dabs(dec(k)-decutd)**2
c		 if(cutrad.gt.4.D0)go to 99
c		 if(dabs(ra(k)-racutd).gt.2.0D0/dcosd(dec(k)).or.dabs(dec(k)-decutd).gt.2.0D0)go to 99
		 ira=720-int(2.*ra(1))
		 jdec=1+int(180.+2.*dec(1))
		 if(ira.gt.720) then
		   ira=720
		 elseif(ira.lt.1)then
		   ira=1
		 endif
		 if(jdec.gt.360)then
		   jdec=360
		 elseif(jdec.lt.1)then
		   jdec=1
		 endif

	       enddo !!end vertex loop
c	       if(ii.eq.159.and.abs(dr).lt.0.5)print *,phi,rlat
c
c	       response=response*(1.D0-0.0092D0*thetsum/4.D0)	!radial binsize, nonlinear term
c	       if(czsum.gt.0.0) response=4.D0*response/czsum  !!Here correct surface brightness for aperture effective area
c******************************Note consequences here of not having loop k=1,4
	       response=response*(1.D0-0.0092D0*thetsum)	!radial binsize, nonlinear term
	       if(czsum.gt.0.0) response=response/czsum  !!Here correct surface brightness for aperture effective area
c******************************	       

c	       tcount=0
c	       ttotal=0
c	       if(ii.le.250)then
c		 CALL INDEXPIXEL(level,node,nodeid,hits,response,ra,dec,tcount,ttotal,iflag,mflag)
c	       endif
c	       BadPix(ii,jj)=BadPix(ii,jj)+tcount
c	       totcnt=totcnt+ttotal
c	       badcnt=badcnt+tcount
c	       if(ttotal.ne.0)then
c	         cutout(ira,jdec)=cutout(ira,jdec)+response*float(tcount)/float(ttotal)
c	         cutcnt(ira,jdec)=cutcnt(ira,jdec)+float(tcount)/float(ttotal)
c	         cutcounter(ira,jdec)=cutcounter(ira,jdec)+1.
c	       endif
c*********************************This section for forming an average skymap
	       coarsesky(ira,jdec)=coarsesky(ira,jdec)+int(response)
	       coarseskyhits(ira,jdec)=coarseskyhits(ira,jdec)+1
c*********************************End section for forming an average skymap
c*********************************This section for forming a least skymap
c	       if(coarseskyhits(ira,jdec).eq.0)then
c	       		coarsesky(ira,jdec)=int(response)
c	       		coarseskyhits(ira,jdec)=coarseskyhits(ira,jdec)+1
c	       elseif(int(response).lt.coarsesky(ira,jdec))then
c	       		coarsesky(ira,jdec)=int(response)
c	       endif
c*********************************End section for forming a least skymap
	       delra=racutd-ra(1)
	       if(delra.gt.350.D0)delra=delra-360.D0
	       if(delra.lt.-350.D0)delra=delra+360.D0
	       ira=idint(delra*dcosd(dec(1))+5.5D0)
	       if(ira.lt.1.or.ira.gt.11)go to 99
	       jdec=idint(dec(1)-decutd+5.5D0)
	       if(jdec.lt.1.or.jdec.gt.11)go to 99
	       ismallmaps(pnum,ira,jdec,mloop)=ismallmaps(pnum,ira,jdec,mloop)+1
	       smallmaps(pnum,ira,jdec,mloop)=smallmaps(pnum,ira,jdec,mloop)+response
99	       continue	       
c*********************************	       
c	       arg=skytime(ira,jdec)/float(coarseskyhits(ira,jdec))
c	       if(orbitfraction/arg.gt.30.)skytime(ira,jdec)=-skytime(ira,jdec)
c	       if(skytime(ira,jdec).lt.0.)skytime(ira,jdec)=skytime(ira,jdec)-orbitfraction
c	       if(skytime(ira,jdec).ge.0.)skytime(ira,jdec)=skytime(ira,jdec)+orbitfraction
	     endif
	     endif !!end FOV if-block
	   endif  !!end non-zero response if-block
	enddo
	enddo !!end pixel by pixel double loop
	
	return
	end
C+
C NAME:
C	rotateq
C PURPOSE:
C	Calculate polar angles in rotated coordinate system
C CATEGORY:
C	Math: coordinate transformation
C CALLING SEQUENCE:
C	call rotateq(Q,PHI,RLAT)
C CALLS:
C	DACOSD,DATAN2D,DCOSD,DSIND,DSQRT,DMOD
C INPUTS:
C	Q(4)	real	quaternion by which rotation is to be performed
C	PHI	real	phase angle in old coordinate system
C	RLAT	real	latitude in old coordinate system (-90<=RLAT<=90)
C OUTPUTS:
C	PHI	real	phase angle in new coordinate system (0<=PHI<360)
C	RLAT	real	latitude in new coordinate system
C PROCEDURE:
C	given a unit quaternion, it specifies a rotation about a unit vector.  This rotation 
C	can be put into matrix form and multiplied by cartesian coordinates to get rotated
C	cartesian coordinates.
C	All angles are in degrees.
C MODIFICATION HISTORY:
C	2003 March, Aaron(UCSD/CASS)
C
C-

	subroutine rotateq(q, phi, rlat)
	implicit none

	real*8 q(4), phi, rlat, xold, yold, zold, xnew, ynew, znew
	real*8 r(3,3), q1sq, q2sq, q3sq, q4sq
	
	!!find cartesian coordinates from ra/dec
	xold = dcosd(rlat)*dcosd(phi)
	yold = dcosd(rlat)*dsind(phi)
	zold = dsind(rlat)

	!!set up rotation matrix r(3,3)
	q1sq=q(1)*q(1)
	q2sq=q(2)*q(2)
	q3sq=q(3)*q(3)
	q4sq=q(4)*q(4)
	r(1,1) = q1sq + q2sq - q3sq - q4sq
	r(1,2) = 2.D0*( q(2)*q(3) + q(1)*q(4) )
	r(1,3) = 2.D0*( q(2)*q(4) - q(1)*q(3) )
	r(2,1) = 2.D0*( q(2)*q(3) - q(1)*q(4) )
	r(2,2) = q1sq - q2sq + q3sq - q4sq
	r(2,3) = 2.D0*( q(3)*q(4) + q(1)*q(2) )
	r(3,1) = 2.D0*( q(2)*q(4) + q(1)*q(3) )
	r(3,2) = 2.D0*( q(3)*q(4) - q(1)*q(2) )
	r(3,3) = q1sq - q2sq - q3sq + q4sq
	
C	print*, 'DCM Matrix:'
C	print*, r(1,1), R(1,2), r(1,3)
C	print*, r(2,1), R(2,2), r(2,3)
C	print*, r(3,1), R(3,2), r(3,3)

	!!rotate X,Y,Z old by matrix multiplication
	xnew = r(1,1)*xold + r(1,2)*yold + r(1,3)*zold
	ynew = r(2,1)*xold + r(2,2)*yold + r(2,3)*zold
	znew = r(3,1)*xold + r(3,2)*yold + r(3,3)*zold

	!!at poles |znew| becomes greater than 1 (by less than 10**-6)
	if(znew.gt.1.D0) znew=1.D0
	if(znew.lt.-1.D0) znew=-1.D0

	phi = datan2d(ynew,xnew)
	phi = dmod(dmod(phi,360.D0)+360.D0,360.D0) !!0<phi<360
	rlat = dasind(znew)

	return
	end
C+
C NAME:
C	plane
C PURPOSE:
C	Removes tilt and offset from a 3 x 3 array of values
C CATEGORY:
C	Data processing: called by subroutine "cosmicra.f"
C CALLING SEQUENCE:
C	call plane(z,z0)
C INPUTS:
C	z(3,3)   	   real		a 3x3 square array of values
C OUTPUTS:
C	z(3,3)		   real		the above array with LSQ fit removed
C	z0		   real		the average value removed from all 9
C FUNCTIONS/SUBROUTINES:   none
C PROCEDURE:
C	Finds a least-squares planar fit through the 8 outside z values,
C	then reduces all 9 z's to values relative to the plane
C MODIFICATION HISTORY:
C	Jan, 2003 A. Buffington (UCSD)
C-
	subroutine plane(z,z0)
	implicit none
c
c	This subroutine solves for a least-squares planar fit through the 8 z-values
c	surrounding [2,2], then reduces all 9 z's relative to the plane,
c	and returns with the reduced values.
c
	real*4 z(3,3),sumz,sumx,sumy,a,b,z0
	integer*4 i,j
c
	sumz=z(1,2)+z(3,2)
	sumx=0.
	sumy=0.
	do i=1,3
		sumz=sumz+z(i,1)+z(i,3)
		sumx=sumx+z(1,i)-z(3,i)
		sumy=sumy+z(i,1)-z(i,3)
	enddo
	a=sumx/6.
	b=sumy/6.
	z0=sumz/8.
	do i=1,3
	    do j=1,3
		z(i,j)=z(i,j)-z0+a*float(i-2)+b*float(j-2)
	    enddo
	enddo
	return
	end
