C+ C NAME: C A_mod_mainf C PURPOSE: C This is the main Fortran program to index SMEI data using utilizing C Spatial Index functions from the c++ file A_mod_L12_all.cpp C CATEGORY: C SMEI data analysis C CALLING SEQUENCE: C INPUTS: C Various parameters are presently specified by data statements at the top C This program presently combines data from a single SMEI camera C OUTPUTS: C A data file listing pedestal and dark current values versus frame number C A grid file containing the resulting sky map C CALLS: C set_flatfields !sets the appropriate flatfields for the selected camera C read_frame !reads a frame of data C lsff !large scale flatfield C ssff !small scale flatfield C rotate8 !euler angle rotation C indexframe !frame indexing routine C SEE ALSO: C INCLUDE: C COMMON BLOCKS: C PROCEDURE: C This program uses the hierarchical triangle coordinate system from Johns Hopkins, C as a base upon which the photometric measurements from a single SMEI camera are C combined and then averaged to produce a surface-brightness sky map. The "level" C parameter governs the granularityof the triangular coordinate frame. The C coordinates of the four corners of a pixel (or bin) are transformed down into C this frame, and the photometric surface brightness as measured by that pixel is C added to all triangles whose centers lie within the rectangle. C This surface brightness is corrected by 1/cosine of the pixel's angle to the camera C axis, and by r/r0 the pixel's1 fractional distance to the FOV's rotational center C When the data sequence is finished, all triangle response sums are C averaged and the resulting sky map binned into an output surface brightness map. C MODIFICATION HISTORY: C 2002,2003 Aaron Smith (UCSD/CASS) C 2002->Original version C 2003->removed most C++ functions and re-wrote in fortran as subroutines C ->add quaternion capabilities C DEC-2002, Paul Hick (UCSD/CASS; pphick@ucsd.edu) C Defined LEVEL and KMAX as parameters. KMAX is determined from LEVEL. C KMAX is then used to dimension arrays node, nodeid and hits C JAN-2003, Andrew Buffington (UCSD/CASS) C Integrated standard pedestal and dark current routines, FORTRAN C Rendered names compliant with those of B.V. Jackson C Incorporated 1/cosine surface-brightness correction for aperture perspective (currently in indexframe routine) C Incorporated large and small scale flatfield corrections C Incorporated FORTRAN read subroutine C- PROGRAM A_mod_mainf implicit none c ! Dimensions for node, nodeid, hits above must be at least 2*(4**(level+1)) ! For LEVEL = 12, KMAX = 134217728 ! For LEVEL = 8, KMAX = 524288 c integer*2 TOP_LEVEL parameter ( TOP_LEVEL = 11 ) c integer*8 TOP_KMAX parameter ( TOP_KMAX = 2*(4**(TOP_LEVEL+1)) ) c 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 / c integer*2 level /TOP_LEVEL/ integer*8 kmax /TOP_KMAX/ c real*4 frame (318,64) real*4 framessff (318,64) real*4 framelsff (318,64) real*4 framelsffu (318,64) real*4 pattern (318,64) real*4 cutout (720,360) /259200*0./ real*4 cutcount (720,360) /259200*0./ real*4 cutcounter (720,360) real*4 skytime (720,360) real*4 ped /0./ real*4 dark /0./ real*4 avg (318,64) /0./ real*8 t0 /86201.9D0/ !time of day in day 365 of 2002 (23h 56m 41.9s) real*8 orbitperiod !orbital period in seconds real*8 orbitbegintime real*8 orbitendtime real*8 affine(6) /6*0.0D0/ ! Affine parameters zero for now real*8 cra real*8 cdec real*8 arg c integer*2 coarseskyhits(720,360) /259200*0/ integer*4 coarsesky(720,360) /259200*0/ !coarse all sky map integer*4 BadPix(318,64) /20352*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 image1(1200,1200) /1440000*0/ integer*4 image2(1200,1200) /1440000*0/ integer*4 image3(1200,1200) /1440000*0/ integer*4 npolarimage(1600,1600) /2560000*0/ integer*4 spolarimage(1600,1600) /2560000*0/ integer*4 s,k integer*4 offset integer*8 tcount integer*8 ttotal c c AARON'S DECLARATIONS character infile*59 character outfile1*44 character outfile2*44 character outfile3*44 character outfile4*44 character outfile5*44 character outfile6*44 character outfile7*44 character outfile8*44 character outfile9*44 character outfile10*44 character listfile*32 character hhead(5)*10 character framename*25 character fname1*6 character fname2*1 character fname3*1 character fname4*4 character dum c integer*4 i,j,ii,jj,iii,jjj integer*4 ss integer*4 iarg integer*4 jarg integer*4 listtop integer*4 meascnt integer*4 icorrcnt integer*4 iffix integer*4 doy integer*4 fyear integer*4 fdoy integer*4 fh integer*4 fm integer*4 fs integer*4 leaps integer*4 globloctest integer*4 orbyear integer*4 orbdoy integer*4 orbhh integer*4 orbmm integer*4 orbss integer*4 firstorbit integer*4 firstdataorbit integer*4 firstlocalorbit integer*4 lastorbit integer*4 lastlocalorbit integer*4 lss integer*4 globalorbit 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 orbitfraction real*4 toppercent real*4 bottompercent c real*8 phi real*8 rlat real*8 f real*8 pixsum real*8 pixdiff real*8 ftime C real*8 x(4), y(4), ra(4), dec(4) real*8 qcam(3,4) real*8 qc(4) real*8 qsc(4) real*8 qandy(4) real*8 qdum(4) real*8 qignore(4) real*8 qfinal(4) 963 continue print*,'Please enter camera number(1 or 2)' read 88,ic globloctest=0 88 format(i6) if(globloctest.eq.1)then print*, 'please enter first orbit number(>64)' read 88,firstorbit print*,'please enter last orbit number(<92)' read 88,lastorbit elseif(globloctest.eq.0)then c print*, 'please enter first global orbit number of first orbit in data section' c print*, '2035 for may2003' c print*, '4185 for Oct2003' c print*, '4622 for Nov2003' c print*, '7218 for may2004' c read 88,firstdataorbit firstdataorbit=2035 print*,'please enter first local orbit number(>64)' read 88,firstlocalorbit print*,'please enter last local orbit number(<92)' read 88,lastlocalorbit lastorbit=firstdataorbit+(lastlocalorbit-1) firstorbit=firstdataorbit+(firstlocalorbit-1) else print*,'invalid entry' go to 963 endif !! 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 print*, 'Fortran in...' if(ic.eq.1)then offset=17 open(10,file='DATA_SMEI/newfltgrd1.GRD',readonly) open(11,file='DATA_SMEI/pattern1h_143.grd',readonly) endif if(ic.eq.2)then offset=11 open(10,file='DATA_SMEI/newfltgrd2.GRD',readonly) open(11,file='DATA_SMEI/pattern2h_144.grd',readonly) endif read(11,9)dum read(11,9)dum read(11,9)dum read(11,9)dum read(11,9)dum do j=1,jd read(11,11)(pattern(i,j),i=1,id) enddo if(ic.eq.1)patdark=2.78865 if(ic.eq.2)patdark=8.18373 9 format(a10) 11 format(20f9.3) close(11) 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 firstorbit = 7218 !Orbit number in all smei- year: 2004 / day: 145 / time: 05h 06m 22.4s !read off begining of list open(100,file='DATA_SMEI/start_times2.txt',readonly) do s=1,firstorbit-1 read(100,1000)globalorbit,orbyear,orbdoy,orbhh,orbmm,orbss,orbitperiod,orbitbegintime 1000 format(i7,2x,i6,2x,i5,4x,i2.2,1x,i2.2,1x,i2.2,2x,f10.3,2x,f17.3) enddo C**************************************************************** !! loop to index many orbits separately do ss=1,lastorbit-(firstorbit-1) ! last orbit that CAN be done is 184 (when all data frames are unzipped from day 152 through day 157 )(???) globalorbit = firstorbit + ss - 1 lss=firstlocalorbit+ss-1 c if(mod(lss,5).eq.1)then if(lss.ge.0)then !!!!use this when no orbits are to be skipped read(100,1000)globalorbit,orbyear,orbdoy,orbhh,orbmm,orbss,orbitperiod,orbitbegintime !should be able to take the values being read in above orbitperiod = 6096.49D0-0.84D-04*float(globalorbit) orbitbegintime=t0+dfloat(globalorbit-1)*6096.49D0-0.84D-04*dfloat((globalorbit-1)*globalorbit)/2.D0 orbitendtime=orbitbegintime+orbitperiod do i=1,720 do j=1,360 coarsesky(i,j)=0 coarseskyhits(i,j)=0 skytime(i,j)=0. cutout(i,j)=0. cutcount(i,j)=0. cutcounter(i,j)=0. if(i.le.id.and.j.le.jd)BadPix(i,j)=0 enddo enddo print 1004,globalorbit !ss 1004 format(' Beginning orbit # ',i6) if(ic.eq.1)then write(listfile,'(A32)') 'DATA_SMEI/cam1_144to157_CRXa.txt' listtop=1000000 !This just has to be at least as big as the number of frames in the list; bigger is okay. !263905 !total lines(frames) in list cam1_296to310_CRXa.txt endif if(ic.eq.2)then write(listfile,'(A32)') 'DATA_SMEI/cam2_144to157_CRXa.txt' listtop=1000000 !see above note !265951 !total lines(frames) in list cam2_296to310_CRXa.txt endif do k=1,kmax node(k)=0.0 nodeid(k)=0 hits(k)=0 enddo open(12,file=listfile,status='old') C open(12,file=listfile,readonly) do k=1,listtop read(12,12)fname1,fyear,fname2,fdoy,fname3,fh,fm,fs,fname4,doy,ped,dark,qdum,meascnt,pixsum,pixdiff 12 format(a6,I4,a1,I3,a1,3I2.2,a4,I3,2f10.3,4f13.9,I6,2f10.2) leaps=int(float((fyear-1)-2000)/4.)-int((float(fyear-1)-2000)/100.)+int((float(fyear-1)-2000)/400.) ftime = (365.D0)*24.D0*3600.D0*(dfloat(fyear)-2002.D0)+24.D0*3600.D0*(dfloat(fdoy)-365.D0+dfloat(leaps)) & +3600.D0*dfloat(fh)+60.D0*dfloat(fm)+dfloat(fs) c**** c**** Andy's note here: what do we think/do about leap year??? **************************************** c**** if(ftime.gt.orbitendtime)then go to 99 !!get out of k frame loop endif orbitfraction = (ftime-orbitbegintime)/orbitperiod c !!use these to pick portion of orbit to do, 1% ~ 15 frames bottompercent=0.00 !! 0.15 to start 15% through orbit toppercent=1.00 !! 0.75 to end 75% through orbit if(bottompercent.le.orbitfraction.and.orbitfraction.le.toppercent)then write(framename,'(a6,I4.4,a1,I3.3,a1,3I2.2,a4)')fname1,fyear,fname2,fdoy,fname3,fh,fm,fs,fname4 if(dark.le.-9999.)print*, 'BAD DARK CURRENT VALUE IN FRAME ',framename if(ped.le.-9999.)print*, 'BAD PEDESTAL VALUE IN FRAME ',framename cutoff = pixdiff + float(meascnt)/25. if(dark.gt.-9999..and.ped.gt.-9999..and.cutoff.le.120.)then !!NOTE: cutoff is still being used as of 5-1-04 if(ic.eq.1) write(infile,'(A34,A25)') '/zappa/zone/2003_149_150/data1/c1/',framename if(ic.eq.2) write(infile,'(A34,A25)') '/zappa/zone/2003_149_150/data1/c2/',framename if(mod(k,100).eq.0) print 1001,lss,k,infile,meascnt,pixdiff,cutoff,orbitfraction 1001 format(i3,i7,': ',a64,i6,2f10.2,5x,f5.3) call read_frame_q(infile,mode,ic,id,jd,frame,qignore) !qignore quaternions never used. Replaced with Andy's quaternions patcorr=dark/patdark 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 if(i.gt.4.and.i.lt.318)arg=arg-pattern(i,j)*patcorr !Remove dark current and pattern ! was still subtracting dark separately as of 4-22-04 (i.e. up to orbit 120 for ic = 1 & 2) 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 !! 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 !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*, day144name,' [',cra,',',cdec,']' call indexframe(mode,ic,id,jd,level,TOP_KMAX,node,nodeid,hits,frame,affine,qfinal,tcount,ttotal,BadPix, & cutout,cutcount,cutcounter,pattern,patdark,coarsesky,coarseskyhits,skytime,orbitfraction) endif !! end if block to specify fractions within an orbit endif !! ftime must be greater than orbitbegintime enddo !!end k frame loop 99 continue close(12) print 1002, lss 1002 format('Orbit #',I3,' indexed. Calling average night') call averagenight_allsky(node,nodeid,hits,TOP_KMAX,level,image1,image2,image3, npolarimage, spolarimage) print*, 'Night Averaged, Writing Files...' if(lss.lt.50000)then do i=1,720 do j=1,360 if(coarseskyhits(i,j).ne.0)then coarsesky(i,j)=coarsesky(i,j)/coarseskyhits(i,j) skytime(i,j)=skytime(i,j)/float(coarseskyhits(i,j)) endif if(cutcount(i,j).ne.0.)cutout(i,j)=cutout(i,j)/cutcount(i,j) if(cutcounter(i,j).ne.0.)cutcount(i,j)=cutcount(i,j)/cutcounter(i,j) enddo enddo 1234 format(A20,I1.1,A4,I4,A1,I3.3,A1,3i2.2,A4) !!Write image1(sky atlas) to file write (outfile1,1234) '/zappa/zone/Orbits/c',ic,'eq1_',orbyear,'_',orbdoy,'_',orbhh,orbmm,orbss,'.grd' open (10,file = outfile1, status = 'new', form = 'formatted') write(10,10) 10 format('DSAA'/'1200 1200'/'0 120'/'-60 60'/'0 1400') do j=1,1200 write(10,'(20i6)')(image1(i,j),i=1,1200) enddo close(10) write (outfile2,1234) '/zappa/zone/Orbits/c',ic,'eq2_',orbyear,'_',orbdoy,'_',orbhh,orbmm,orbss,'.grd' open (13,file = outfile2, form = 'formatted') write(13,13) 13 format('DSAA'/'1200 1200'/'120 240'/'-60 60'/'0 1400') do j=1,1200 write(13,'(20i6)')(image2(i,j),i=1,1200) enddo close(13) write (outfile3,1234) '/zappa/zone/Orbits/c',ic,'eq3_',orbyear,'_',orbdoy,'_',orbhh,orbmm,orbss,'.grd' open (14,file = outfile3, form = 'formatted') write(14,14) 14 format('DSAA'/'1200 1200'/'240 360'/'-60 60'/'0 1400') do j=1,1200 write(14,'(20i6)')(image3(i,j),i=1,1200) enddo close(14) write (outfile4,1234) '/zappa/zone/Orbits/c',ic,'npl_',orbyear,'_',orbdoy,'_',orbhh,orbmm,orbss,'.grd' open (15,file = outfile4, form = 'formatted') write(15,15) 15 format('DSAA'/'1600 1600'/'-40 40'/'-40 40'/'0 1400') do j=1,1600 write(15,'(20i6)')(npolarimage(i,j),i=1,1600) enddo close(15) write (outfile5,1234) '/zappa/zone/Orbits/c',ic,'spl_',orbyear,'_',orbdoy,'_',orbhh,orbmm,orbss,'.grd' open (16,file = outfile5, form = 'formatted') write(16,16) 16 format('DSAA'/'1600 1600'/'-40 40'/'-40 40'/'0 1400') do j=1,1600 write(16,'(20i6)')(spolarimage(i,j),i=1,1600) enddo close(16) write (outfile6,1234) '/zappa/zone/Orbits/c',ic,'bpx_',orbyear,'_',orbdoy,'_',orbhh,orbmm,orbss,'.grd' open (17,file = outfile6, form = 'formatted') write(17,17) 17 format('DSAA'/'318 64'/'0 318'/'0 64'/'0 20000') do j=1,64 write(17,'(20i8)')(BadPix(i,j),i=1,318) enddo close(17) write (outfile7,1234) '/zappa/zone/Orbits/c',ic,'cut_',orbyear,'_',orbdoy,'_',orbhh,orbmm,orbss,'.grd' open (18,file = outfile7, form = 'formatted') write(18,18) 18 format('DSAA'/'720 360'/'0 360'/'-90 90'/'0 1400') do j=1,360 write(18,'(20f10.2)')(cutout(i,j),i=1,720) enddo close(18) write (outfile8,1234) '/zappa/zone/Orbits/c',ic,'cnt_',orbyear,'_',orbdoy,'_',orbhh,orbmm,orbss,'.grd' open (19,file = outfile8, form = 'formatted') write(19,19) 19 format('DSAA'/'720 360'/'0 360'/'-90 90'/'0 10') do j=1,360 write(19,'(20f8.2)')(cutcount(i,j),i=1,720) enddo close(19) write (outfile9,1234) '/zappa/zone/Orbits/c',ic,'asm_',orbyear,'_',orbdoy,'_',orbhh,orbmm,orbss,'.grd' open (20,file = outfile9, form = 'formatted') write(20,20) 20 format('DSAA'/'720 360'/'0 360'/'-90 90'/'0 1400') do j=1,360 write(20,'(20i8)')(coarsesky(i,j),i=1,720) enddo close(20) write (outfile10,1234) '/zappa/zone/Orbits/c',ic,'tmp_',orbyear,'_',orbdoy,'_',orbhh,orbmm,orbss,'.grd' open (21,file = outfile10, form = 'formatted') write(21,21) 21 format('DSAA'/'720 360'/'0 360'/'-90 90'/'0 1') do j=1,360 write(21,'(20f8.3)')(skytime(i,j),i=1,720) enddo close(21) endif !!don't waste time writing files when testing endif enddo !!end ss loop to separately index many orbits close(100) print*, ' ' 98 continue print*, 'All Orbits Finished!!!' END