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 (636,128) real*4 pattern (636,128) 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 /2/ !onboard binning mode integer*4 ic !Camera # integer*4 ifm /0/ !ROI mask (0=no,1=yes) integer*4 id /636/ !1280 for TMO integer*4 jd /128/ !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 character hhead(5)*10 character framename*25 character*7 ydh,ydhsav,ydhlist(4),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(4) integer*4 ismallmaps(4,11,11,3) c real*4 acorr(1280,301) real*4 acorrsum real*4 flatcorr(636,128) real*4 patdark real*4 patcorr real*4 cutoff !,glare,glaremap(318,64),glarescale real*4 racut(4),decut(4) real*4 smallmaps(4,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/Buffer3/Nic288/c3m1_2004_DOY_hhmmss.raw'/ data infile /'G:/Grayburst/Orbit#/c3m1_200Y_DOY_hhmmss.nic'/ data kmin /4*0/ data smallmaps /1452*0./ data ismallmaps /1452*0/ data nines /999.9/ data ydhlist /'3_085_0','4_363_0','5_029_1','5_065_0'/ data phicknum /64,2,81,88/ data one,two/'1','2'/ data decut / -11.72,3.05,-3.06,-9.16/ data racut / 292.97,336.64,252.8,282.32/ ic=3 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 if(ic.eq.3)offset=13 c open (13,file='patt3_365_2x2.grd',readonly) print *,'reading in patt3_365_2x2.grd' read(13,13)hhead 13 format(a10) do j=1,128 read(13,15)(pattern(i,j),i=1,636) 15 format(20f9.3) enddo read(13,1500)patdark 1500 format(f10.4) print *,' patdark = ',patdark close(13) c open (10,file='D:\Ccd\birmingham\flatfield\newfltgrd3.GRD',readonly) print *,'reading in D:\Ccd\birmingham\flatfield\newfltgrd3.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 (12,file='G:\grayburst\framelist3.txt') c c glare151_13.grd c c open(10,file='glare151_13.grd') c read (10,1)hhead c read(10,110)((glaremap(i,j),i=1,318),j=1,64) 110 format(20f9.2) c close(10) c glarescale=1.5 c kmax(1)=214 kmax(2)=216 kmax(3)=205 c kmax(1)= 1500 c kmax(2)= 3000 !3150 c kmax(3)= 4500 !4648 c kmax(4)= 4648 c kmin(1) = 0 c kmin(2) = 1501 !1550 c kmin(3) = 3001 !3150 c kmin(4) = 0 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.or.mloop.eq.3)then carg=ydhlist(1) write(carg(7:7),'(a1)')one ydhlist(1)=carg carg=ydhlist(2) write(carg(7:7),'(a1)')one ydhlist(2)=carg if(mloop.eq.3)then carg=ydhlist(3) write(carg(7:7),'(a1)')two ydhlist(3)=carg endif endif ksect=1 ydhsav='0000000' pnum=1 write(11,1001) 1001 format(/) do i=1,720 do j=1,360 coarsesky(i,j)=0 coarseskyhits(i,j)=0 enddo enddo c open(17,file='cam3_287to288_CRX.txt',readonly) if(mloop.eq.1)open(17,file='G:\grayburst\cam3_grayburst1_CRX.txt',readonly) if(mloop.eq.2)open(17,file='G:\grayburst\cam3_grayburst2_CRX.txt',readonly) if(mloop.eq.3)open(17,file='G:\grayburst\cam3_grayburst3_CRX.txt',readonly) c if(mloop.eq.4)open(17,file='G:\grayburst\cam3_grayburst4_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)framename,ped,dark,glare read(17,170)infile,ped,dark !grayburst c170 format(1x,a24,3x,3f10.3) 170 format(a44,2x,2f10.3) c if(k.ne.3119.and.k.ne.4421)go to 99 !does individual frames c write(infile(24:44),'(a21)')framename !new naming convention c write(ydh,'(a7)')infile(29:35) if(ydh.ne.ydhsav.or.k.eq.kmax(mloop))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 if(ydh.ne.ydhlist(pnum))pnum=pnum+1 c print *,'Beginning ',ydh endif endif racutd=dfloat(racut(pnum)) decutd=dfloat(decut(pnum)) cutoff=0. 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,10).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) c the camera 3 gain doubled at 17:30 UT on 2-24-05, 2005_055... if(pnum.ge.4)arg=0.5*arg c arg = arg/0.75 !be careful here, omit headroom adjustment arg=arg-ped !Remove pedestal c if(i.ge.9.and.i.le.634)arg=arg-pattern(i,j)*patcorr !Remove dark current and pattern c if(i.ge.9.and.i.le.634)arg=arg-0.66*dark*(1.+float(j)/128.) !linear slope removal, naive... if(i.ge.9.and.i.le.634)arg=arg-0.5*dark*(1.+2.*float(j)/128.) !linear slope removal, some steeper if(flatcorr(i,j).ne.0.)arg=arg/flatcorr(i,j) !Flat Field Correction c ii=(i+1)/2 c jj=(j+1)/2 c if(frame(i,j).ne.0..and.glaremap(ii,jj).gt.0.)arg=arg-glarescale*glaremap(ii,jj) !remove glare if(frame(i,j).ne.0.)frame(i,j)=arg enddo enddo c c 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 !! 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\cam3gray_orbit1a.grd', form = 'formatted') c if(mloop.eq.2)open (20,file = 'g:\grayburst\cam3gray_orbit2a.grd', form = 'formatted') c if(mloop.eq.3)open (20,file = 'g:\grayburst\cam3gray_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\cam3_orbit1g1p5.grd', form = 'formatted') c if(mloop.eq.2)open (20,file = 'g:\Sean\Eclipse\cam3_orbit2g1p5.grd', form = 'formatted') c if(mloop.eq.3)open (20,file = 'g:\Sean\Eclipse\cam3_orbit3g1p5.grd', form = 'formatted') c if(mloop.eq.4)open (20,file = 'g:\Sean\Eclipse\cam3_orbit3e.grd', form = 'formatted') c write(20,20) 20 format('DSAA'/'720 360'/'0 360'/'-90 90'/'0 500') 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,4 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\cam3smallmaps.grd') c write(20,21) 21 format('DSAA'/'48 35'/'0.5 4.5'/'0.5 3.5'/'0 500') do jj=1,3 do j=1,11 c write(20,22)((smallmaps(ii,i,j,jj),i=1,11),nines,ii=1,4) 22 format(20f6.1) enddo c if(jj.lt.3)write(20,22)(nines,i=1,48) 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(636,128),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(4,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(4,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 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