C+ C NAME: C A_mod_mainf C PURPOSE: C This is the main Fortran program to index TMO/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 pedestalr !removes pedestal from empty pixels C darkr !removes dark current from covered pixels C cosmicra !cosmic ray removal 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 ! Dimensions for node, nodeid, hits above must be at least 2*(4**(level+1)) ! For LEVEL = 12, KMAX = 134217728 ! For LEVEL = 8, KMAX = 524288 integer*2 TOP_LEVEL parameter ( TOP_LEVEL = 12 ) integer*8 TOP_KMAX parameter ( TOP_KMAX = 2*(4**(TOP_LEVEL+1)) ) 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 / integer*2 level /TOP_LEVEL/ integer*8 kmax /TOP_KMAX/ 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 degperminute 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 iframemask (1272,256) /325632*1/ integer*4 ihit(1272,256) integer*4 ihist(250) 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*8 k c c AARON'S DECLARATIONS c character infile*47,listfile*23,framename*28 character outfile1*20,outfile2*20,outfile3*20, outfile4*18, outfile5*18 integer*4 idarkcr, icra, meascnt, date, hour, maxk, mink integer*4 i, j, ss real*8 phi, rlat, f C real*8 x(4), y(4), ra(4), dec(4) real*8 qcam(3,4), qc(4), qsc(4), qandy(4), qdum(4), qfinal(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 data ihist /250*0/ data ihit /81408*0/ meascnt=0 print*, 'Fortran In ...' c print*, 'Level: ', level, ' (array size: ', kmax,')' adu2e=4.5*65536./hroom c call set_flatfields(ic,framessff,framelsff,framelsffu) !note dimensions should be checked before ff is done... c call set_cr_backgnd(mode,ic,ifm,id,jd,avg) !! loop to index many orbits separately do ss=8,8 if(ss.eq.1) then date=0227 hour=99 maxk=462 mink=310 write(listfile,'(A15,I4.4,A4)') 'DATA_SMEI/list_',date,'.txt' elseif(ss.eq.2)then date=0227 hour=18 maxk=1390 mink=1387 write(listfile,'(A15,I4.4,A4)') 'DATA_SMEI/list_',date,'.txt' elseif(ss.eq.3)then date=0301 hour=03 maxk=2790 mink=1862 write(listfile,'(A15,I4.4,A4)') 'DATA_SMEI/list_',date,'.txt' elseif(ss.eq.4)then date=0301 hour=08 maxk=3692 mink=2794 write(listfile,'(A15,I4.4,A4)') 'DATA_SMEI/list_',date,'.txt' elseif(ss.eq.5)then date=0301 hour=10 maxk=4610 mink=3692 write(listfile,'(A15,I4.4,A4)') 'DATA_SMEI/list_',date,'.txt' elseif(ss.eq.6)then date=0301 hour=12 maxk=5537 mink=4610 write(listfile,'(A15,I4.4,A4)') 'DATA_SMEI/list_',date,'.txt' elseif(ss.eq.7)then date=0301 hour=13 maxk=6464 mink=5537 write(listfile,'(A15,I4.4,A4)') 'DATA_SMEI/list_',date,'.txt' elseif(ss.eq.8)then date=0227 hour=99 maxk=3641 mink=1 write(listfile,'(A15,I4.4,A4)') 'DATA_SMEI/list_',date,'.txt' elseif(ss.eq.9)then date=0301 hour=99 maxk=7307 mink=1 write(listfile,'(A15,I4.4,A4)') 'DATA_SMEI/list_',date,'.txt' endif do k=1,kmax node(k)=0.0 nodeid(k)=0 hits(k)=0 enddo open(12,file=listfile,status='old') do k=1,maxk read(12,12)framename 12 format(a28) if(k.ge.mink.and.mod(k,15).eq.0)then write(infile,'(A14,I4.4,A1,A28)') 'DATA_SMEI/2003',date,'/',framename C print*,k,': ', infile call read_framecrx(infile,mode,ic,id,jd,frame,qdum,meascnt) !! 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 car/cdec values calculated from thetax=thetay=0 cra=0.D0 cdec=90.D0 call rotateq(qfinal,cra,cdec) !!rotate to standard system if(cdec.lt.-60.0) print*, framename,cra,cdec c call pedestalr(mode,ic,ifm,id,jd,iframemask,frame,ped) c call darkr(mode,ic,ifm,id,jd,iframemask,frame,dark,adu2e,idarkcr,ihit,ihist) c call cosmicra(mode,ic,ifm,id,jd,iframemask,frame,avg,adu2e,icrcnt,ihist,ihit,neighbors) c print 1000,ped,dark,idarkcr,icrcnt,neighbors c1000 format(' Frame ',i6,', Ped = ',f8.2,', Dark = ',f8.2,', cosmic rays ',3i5) c write(11,11)ped,dark,idarkcr,icrcnt,neighbors c11 format(i6,2f10.2,3i7) c call lsff(mode,ic,ifm,id,jd,iframemask,framelsff,frame) c call ssff(mode,ic,ifm,id,jd,iframemask,framessff,frame) C call indexframe(mode,ic,id,jd,level,TOP_KMAX,node,nodeid,hits,frame,affine,qfinal) endif !!end if block to include only one orbit enddo !!end k frame loop close(12) print*, 'Calling average night' C call averagenight_allsky(node,nodeid,hits,TOP_KMAX,level,image1,image2,image3, npolarimage, spolarimage) print*, 'Night Averaged, Writing Files...' if(ss.lt.0)then !!Write image1(sky atlas) to file write (outfile1,'(I4.4,A1,I2.2,A13)') date,'_',hour,'_EqPlot_1.grd' open (10,file = outfile1, 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,'(I4.4,A1,I2.2,A13)') date,'_',hour,'_EqPlot_2.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,'(I4.4,A1,I2.2,A13)') date,'_',hour,'_EqPlot_3.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,'(I4.4,A1,I2.2,A11)') date,'_',hour,'_npolar.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,'(I4.4,A1,I2.2,A11)') date,'_',hour,'_spolar.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) endif !!don't waste time writing files when testing enddo !!end ss loop to separately index many orbits print*, ' ' 99 print*, 'All Orbits Finished!!!' END