program ccd_to_sky
	implicit none
	
	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
	character	listfile*21,infile*47,outfile*17,framename*47
	
	real*8 qcam(3,4), qc(4), qsc(4), qandy(4), qdum(4), qfinal(4)
	real*4 yoffset(5), xoffset(5)
	real*8 x(4), y(4), ra(4), dec(4), q(4), qq(4), xccd, yccd, aq, bq, cq
	real*8 rq, dr, thetax, thetaxsq, thetay, cx, cy, cz, crsq, czsum, cra, cdec
	real*8 x0(5), y0(5),x00(5), y00(5), yref(5), r0, tempx, phi, rlat, rccd, arg, arg1, response, dusk
	integer*4 ii, jj
	integer*4 itemp, jtemp, itempmin, itempmax, jtempmin, jtempmax, rcount
	integer*2 k, kk, level, zz

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/
	
	r0=1194.5D0
	x0(ic+1)=x00(ic+1)-xoffset(ic+1)
	y0(ic+1)=y00(ic+1)-yoffset(ic+1)

	!!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
	!! (9) updated test quaternions (e1=-86.10, e2=59.36, e3=-6.25)
	qcam(2,1) = -0.265453D0
	qcam(2,2) = 0.0183744D0
	qcam(2,3) = 0.962964D0
	qcam(2,4) = 0.0435699D0

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

C	write(listfile,'(A21)') 'DATA_SMEI/00/list.txt'
C	write(listfile,'(A21)') 'DATA_SMEI/00/list.txt'
C	open(12,file=listfile,status='old')

	do kk=1,2

	if(kk.eq.1)then
	  write(framename,'(A47)') 'DATA_SMEI/20030227/c2frm_2003_058_180251crx.bin'
	  x(1) = 732.D0
	  y(1) = 50.D0
	else
	  write(framename,'(A47)') 'DATA_SMEI/20030228/c2frm_2003_059_174515crx.bin'
	  x(1) = 730.D0
	  y(1) = 56.D0
	endif

C	read(12,12)framename
C12	format(a25)

	        write(infile,'(A47)') framename
	        print*, kk, infile

		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)
		
C		find ra/dec of center of FOV for a frame
C		    cra=0.D0
C		    cdec=90.D0
C	            call rotateq(qfinal,cra,cdec) !!rotate to standard system
C		    print*, 'Center RA/DEC:  ',cra,cdec

		do k=1,1
C		 !!affine transformation to std. sky coordinates
C		 tempx = affine(1) + (1.D0+affine(2))*x(k) + affine(3)*y(k);
C		 y(k) = affine(4) + affine(5)*x(k) + (1.D0+affine(6))*y(k);
C		 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) !!/1.013D0 !!y-plate scale adjustment

		 thetax = dasind( (x(k) - x0(ic+1)) / rccd );
		 
		 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 )
c		 print*, 'ThetaY 1: ',thetay

c	         thetay = rq / -18.0226D0
c		 print*, 'ThetaY 2: ',thetay

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

		 call rotateq(qfinal, phi, rlat) !!rotate to standard system using original qq quaternion
		 ra(k)=phi
		 dec(k)=rlat

		 print*, 'X/Y:    ',x(k),y(k)
		 print*, 'RA/DEC: ',ra(k),dec(k)
		 
		 !!inverse quaternion for inverse rotation
		 q(1)=qfinal(1)
	  	 q(2)=-qfinal(2)
	  	 q(3)=-qfinal(3)
	  	 q(4)=-qfinal(4)

		 !!find ccd coordinates of star
	  	 call rotateq_sky_to_ccd(q,ra(k),dec(k), cx, cy, cz)

	  	 thetax=dasind(cx)
	  	 thetay=dasind(cy)

	  	 rccd=-18.2569D0*thetay + r0

	  	 xccd=rccd*dsind(thetax) + x0(ic+1)
	  	 yccd=y0(ic+1)-rccd*dcosd(thetax)

	 	 print*, 'X,Y:    ',xccd,yccd

		enddo

	enddo

C	close(12)
	end