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 glaremap (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 cra real*8 cdec 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 fisheye(360,360) /129600*0/ !coarse North fisheye map integer*4 ifisheye(360,360) /129600*0/ 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,kmax(4),mloop,ifish,jfish integer*4 offset integer*4 subtsky(720,360) c c AARON'S DECLARATIONS character infile*50 !53 character hhead(5)*10 character framename*25 c integer*4 i,j,ii,jj,iii,jjj,icorrcnt integer*4 iarg integer*4 jarg integer*4 iffix,ibin,jbin,jjbin integer*4 mapuse(90,180) 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 rx real*4 ry c real*4 pixsum c real*4 pixdiff real*4 glare,glareuse,glarepix,r,theta,x,y,glarecam2(90,180) c real*8 phi real*8 rlat 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) data mapuse /16200*0/ ic=2 c c print*,'Please enter camera number' c read 88,ic c88 format(f7.4) !! 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 c open (13,file='c2pat_2003_326_4x4.grd',readonly) !DOY327 c print *,'reading in c2pat_2003_326_4x4.grd' !DOY327 c open (13,file='c2pat_2003_144_4x4.grd',readonly) !DOY151 c print *,'reading in c2pat_2003_144_4x4.grd' !DOY151 open (13,file='c2pat_2004_092_4x4.grd',readonly) !DOY110-2004 print *,'reading in c2pat_2004_092_4x4.grd' !DOY110-2004 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 if(ic.eq.2) then open(37,file='glare151_13.grd',readonly) print *,'reading in glare151_13.grd' read(37,13)hhead read(37,37)((glaremap(i,j),i=1,id),j=1,jd) 37 format(20f9.2) close(37) endif 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) c print *,'reading in G:\Sean\cam2allz.grd' open(20,file='cam2all.grd') read(20,1)hhead do j=1,360 read(20,201)(subtsky(i,j),i=1,720) 201 format(20i8) enddo close(20) c print *,'reading in G:\Sean\glaremapcam2.grd' open(22,file='glaremapcam2.grd') 22 format(20f8.1) !empirical map c print *,'reading in G:\Sean\aglaremodel.grd' c open(22,file='aglaremodel.grd') c22 format(30f5.2) !glare model read(22,1)hhead do j=1,180 read(22,22)(glarecam2(i,j),i=1,90) enddo close(22) c c Day 168 1 to 713 in G c Day 174 714 to 1407 in G c Day 180 4898 to 5610 in G kmax(1)= 44542 kmax(2)= 9160 kmax(3)= 23675 kmax(4)= 5610 !7038 kmin = 4897 !00 do i=1,720 do j=1,360 coarsesky(i,j)=0 coarseskyhits(i,j)=0 if(i.le.360)then fisheye(i,j)=0 ifisheye(i,j)=0 endif enddo enddo c c***************************Begin main loops c do mloop = 4,4 if(mloop.eq.1)open(17,file='glarecrx2b.txt',readonly) if(mloop.eq.2)open(17,file='glarecrx2c.txt',readonly) if(mloop.eq.3)open(17,file='glarecrx2f.txt',readonly) if(mloop.eq.4)open(17,file='glarecrx2g.txt',readonly) c c c wind off unwanted frames at the beginning c if(kmin.gt.0)then do k=1,kmin read(17,170)infile enddo endif c do k=kmin+1,kmax(mloop) read(17,170)infile,ped,dark,glare,rx,ry,glarepix 170 format(a50,6f10.3) c c if(glare.gt.50.)glare=50. !questionable c if(glare.gt.50.)go to 99 !bail if bad glare value c if(glarepix.gt.40.)glarepix=40. !equally questionable c if(glarepix.gt.40.)go to 99 !bail if bad glarepix value c c 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 c cutoff = pixdiff + float(meascnt)/25. c if(dark.gt.-9999..and.ped.gt.-9999.)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 ibin=1+int(abs(rx)) jbin=91+int(abs(ry)) !top half of table only, reflected c jbin=int(91-ry) !table NOT reflected, bad idea! glareuse=glarecam2(ibin,jbin) mapuse(ibin,jbin)=1 if(glareuse.eq.0.)then print *,k,', zero glareuse: rx,ry = ',rx,ry jjbin=90-int(abs(ry)) if(glarecam2(ibin,jjbin).ne.0.)then glareuse=glarecam2(ibin,jjbin) print *,'OK, other glareuse = ',glareuse endif endif c print *,k,glare/15.5,glarepix/15.5,glareuse,rx,ry c do i=1,id do j=1,jd if(mloop.eq.1.and.k.ge.31413.and.k.le.34230)frame(188,10)=0. !kill a bad pixel... 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(flatcorr(i,j).ne.0.)arg=arg/flatcorr(i,j) !Flat Field Correction c arg=arg-glarepix*glaremap(i,j)/15.5 !Glare correction, 2-pix map c arg=arg-glare*glaremap(i,j)/15.5 !Glare correction, 180-pix band arg=arg-glareuse*glaremap(i,j) 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) c endif 99 continue enddo !end k frame loop close(17) enddo !end mloop on 4 data sections c do i=1,720 do j=1,360 if(coarseskyhits(i,j).ne.0.and.subtsky(i,j).ne.0)then coarsesky(i,j)=coarsesky(i,j)/coarseskyhits(i,j) coarsesky(i,j)=coarsesky(i,j)-subtsky(i,j) r=float(j)/2.-180. theta=float(i)/2. x=r*cosd(theta) y=r*sind(theta) ifish=nint(180.+2.*x) jfish=nint(180.+2.*y) if(ifish.gt.0.and.ifish.le.360.and.jfish.gt.0.and.jfish.le.360)then ifisheye(ifish,jfish)=ifisheye(ifish,jfish)+1 fisheye(ifish,jfish)=fisheye(ifish,jfish)+coarsesky(i,j) endif else coarsesky(i,j)=0 endif enddo enddo open (21,file = 'mapuse180.grd', form = 'formatted') write(21,21) 21 format('DSAA'/'90 180'/'0 90'/'-90 90'/'0 1') do j=1,180 write(21,'(90i2)')(mapuse(i,j),i=1,90) enddo close(21) open (20,file = 'cam2sub180.grd', form = 'formatted') write(20,20) 20 format('DSAA'/'720 360'/'0 360'/'-90 90'/'0 50') do j=1,360 write(20,'(20i8)')(coarsesky(i,j),i=1,720) enddo close(20) open (30,file = 'fish2sub180.grd', form = 'formatted') write(30,30) 30 format('DSAA'/'360 360'/'-90 90'/'-90 90'/'0 50') do j=1,360 do i=1,360 if(ifisheye(i,j).gt.0)fisheye(i,j)=fisheye(i,j)/ifisheye(i,j) enddo write(30,'(20i8)')(fisheye(i,j),i=1,360) enddo close(30) 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*50 !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) 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*8 x(4), y(4), ra(4), dec(4), euler(4), q(4), affine(6), aq, bq, cq 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) 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 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 c coarsesky(ira,jdec)=coarsesky(ira,jdec)+int(response) c coarseskyhits(ira,jdec)=coarseskyhits(ira,jdec)+1 c*********************************This section for forming a least skymap if(coarseskyhits(ira,jdec).eq.0)then coarsesky(ira,jdec)=int(response) coarseskyhits(ira,jdec)=coarseskyhits(ira,jdec)+1 elseif(int(response).lt.coarsesky(ira,jdec))then coarsesky(ira,jdec)=int(response) endif 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