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 = 11 ) 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/ C real*4 frame (1272,256) C real*4 framessff (1272,256) C real*4 framelsff (1272,256) C real*4 framelsffu (1272,256) real*4 frame (318,64) real*4 frame2 (318,64) real*4 framessff (318,64) real*4 framelsff (318,64) real*4 framelsffu (318,64) real*4 pattern (318,64) real*4 ped /0./ real*4 dark /0./ C real*4 avg (1272,256) /0./ real*4 avg (318,64) /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 /4/ !onboard binning mode integer*4 ic /2/ !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 icrcnt /0/ integer*4 neighbors /0/ C integer*4 iframemask (1272,256) /325632*1/ C integer*4 ihit(1272,256) integer*4 iframemask (318,64) /20352*1/ integer*4 ihit(318,64) 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*4 ikbad /1/ integer*4 badk,imaxaurora,offset,k integer*8 tcount !changed k from integer*8 to Integer*4 integer*8 ttemp /0/ character kbad(10000)*6 real*4 response(10000),tempresp c c AARON'S DECLARATIONS c character infile*44,infile2*47, listfile*31,framename*25, mytime*6, badtime*6, dum1*2, dum2*2, dum3*2, dum character outfile1*20,outfile2*20,outfile3*20, outfile4*23, outfile5*23 character hhead(5)*10 integer*4 idarkcr, icra, meascnt, date, hour, maxk, mink(72) integer*4 i, j, ss, fcnt,iarg,jarg,ii,jj,iii,jjj,icorrcnt,iffix,doy real*4 acorr(1280,301),acorrsum,flatcorr(318,64) real*8 phi, rlat, f, pixsum, pixdiff C real*8 x(4), y(4), ra(4), dec(4) real*8 qcam(3,4), qc(4), qsc(4), qandy(4), qdum(4), qignore(4), qfinal(4) c**********************This data for camera 2 data mink / 1, 950, 2437, 3926, 5415, 6901, 8351, 9839, 11326, 12814, & 14300, 15787, 17270, 18701, 20150, 21635, 23121, 24608, 26097, 27582, & 29049, 30535, 32024, 33510, 34996, 36483, 37949, 39413, 40850, 42309, & 43777, 45263, 46750, 48239, 49708, 51178, 52669, 54159, 55648, 57136, & 58621, 60058, 61486, 62954, 64438, 65926, 67412, 68900, 70369, 71835, & 73325, 74815, 76303, 77790, 79275, 80725, 82157, 83626, 85114, 86601, & 88090, 89576, 91046, 92513, 94003, 95493, 96980, 98468, 99954,101410,102839,104299/ c**********************This data for camera 1 c data mink / 1, 976, 2498, 4019, 5544, 7066, 8552, 10076, 11600, 13124, c & 14648, 16171, 17695, 19166, 20645, 22160, 23681, 25205, 26728, 28250, c & 29757, 31263, 32786, 34309, 35833, 37357, 38862, 40368, 41846, 43329, c & 44819, 46342, 47864, 49389, 50893, 52397, 53922, 55445, 56968, 58492, c & 60017, 61489, 62955, 64442, 65967, 67492, 69015, 70539, 72046, 73550, c & 75075, 76599, 78123, 79646, 81169, 82660, 84127, 85616, 87139, 88662, c & 90186, 91694, 93219, 94723, 96248, 97770, 99294,100819,102343,103844,105293,106773/ c*********************** data ihist /250*0/ C data ihit /81408*0/ data ihit /20352*0/ !! 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 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) c imaxaurora=0 c if(ic.eq.1)then c imaxaurora=2747 c open(19,file='aurora1.txt') c endif c if(ic.eq.2)then c imaxaurora=5018 !2686 c open(19,file='aurora.txt') c endif c if(imaxaurora.gt.0)then c do i=1,imaxaurora c read(19,19)dum1, dum2, dum3,tempresp ! Not reading in date, only time of day. If no aurora for >24 hours, could have problem. c19 format(10x,a2,1x,a2,1x,a2,12x,f7.2) c write(badtime,'(a2,a2,a2)')dum1,dum2,dum3 c kbad(i)=badtime c response(i)=tempresp c tempresp=0 cc print*, response(i) c enddo c close(19) c endif if(ic.eq.1)then offset=17 open(10,file='DATA_SMEI/newfltgrd1.GRD',readonly) open(11,file='DATA_SMEI/pattern1.grd',readonly) endif if(ic.eq.2)then offset=11 open(10,file='DATA_SMEI/newfltgrd2.GRD',readonly) open(11,file='DATA_SMEI/pattern2.grd',readonly) endif C if(ic.eq.3)then C offset=?? C open(10,file='DATA_SMEI/newfltgrd3.GRD',readonly) C open(11,file='DATA_SMEI/pattern3.grd',readonly) C endif read(11,9)dum read(11,9)dum read(11,9)dum read(11,9)dum read(11,9)dum do j=1,64 read(11,11)(pattern(i,j),i=1,318) enddo 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) !! loop to index many orbits separately do ss= 26,26 !!22,71 ! last orbit that CAN be done is 71 print 1004,ss 1004 format(' Beginning orbit # ',i2) maxk = mink(ss+1)-1 if(ic.eq.1)write(listfile,'(A31)') 'DATA_SMEI/cam1_148to152_CRX.txt' if(ic.eq.2)write(listfile,'(A31)') 'DATA_SMEI/cam2_148to152_CRX.txt' if(ic.eq.3)write(listfile,'(A31)') 'DATA_SMEI/cam3_148to152_CRX.txt' 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,maxk read(12,12)framename,doy,ped,dark,qdum,meascnt,pixsum,pixdiff 12 format(a25,I3,2f10.3,4f13.9,I6,2f10.2) if(k.ge.mink(ss).and.dark.le.-9999.)print*, 'BAD DARK CURRENT VALUE IN FRAME ',framename if(k.ge.mink(ss).and.ped.le.-9999.)print*, 'BAD PEDESTAL VALUE IN FRAME ',framename badk=0 write(mytime,'(a6)')framename(16:21) C if(mytime.eq.kbad(ikbad))then C if(response(ikbad).gt.5.) badk=1 C if(ikbad.lt.imaxaurora) ikbad=ikbad+1 C endif C if(k.ge.mink(ss).and.badk.eq.1) print*, mytime, response(ikbad-1) c if(k.ge.mink(ss).and.badk.eq.0)then !Choose badk.eq.1 to do aurora frames C if(k.ge.mink(ss).and.dark.gt.-9999..and.ped.gt.-9999..and.mytime.eq.'000003'.and.doy.eq.148)then if(k.ge.mink(ss).and.dark.gt.-9999..and.ped.gt.-9999.)then write(infile,'(A13,I1.1,A1,I3.3,A1,A25)') 'DATA_SMEI/cam',ic,'_',doy,'/',framename C write(infile2,'(A47)') 'DATA_SMEI/cam2_148/c2frm_2003_148_000003crx.bin' if(mod(k,100).eq.0) print 1001,ss,k,infile 1001 format(i2,i7,': ',a47) call read_frame_q(infile,mode,ic,id,jd,frame,qignore) !qignore quaternions never used. Replaced with Andy's quaternions C call read_framecrx_mode(infile2,mode,ic,id,jd,frame2,qdum,meascnt) C call read_framecrx_old_i2(infile,mode,ic,id,jd,frame,qdum) do i=1,id do j=1,jd arg=frame(i,j) arg = arg/0.75 !Headroom adjustment (may noy comply with engineering data and cam3) arg=arg-ped !Remove pedestal if(i.gt.4.and.i.lt.318)arg=arg-dark-pattern(i,j) !Remove dark current and pattern if(flatcorr(i,j).ne.0.)arg=arg/flatcorr(i,j) !Flat Field Correction if(frame(i,j).ne.0.)frame(i,j)=arg C if(i.eq.159)then C print 1003, frame(i,j), frame2(i,j), i, j,frame(i,j)-frame2(i,j),pattern(i,j) C1003 format('NIC = ',f9.2, ', BIN = ',f9.2 ', (i,j) = (',I3,',',I3,') ,NIC-BIN = ',f9.2,', Pattern =',f9.2) C endif enddo enddo C go to 98 c tcount=0 c do i=1,318 c do j=1,64 c if(frame(i,j).gt.0) then c tcount=tcount+1 c if(mod(tcount,10).eq.0) print*, i,j,frame(i,j) c endif c enddo c enddo c print*, 'Count: ',tcount !! 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*, k,' RA/DEC: ',cra,cdec tcount=1 call indexframe(mode,ic,id,jd,level,TOP_KMAX,node,nodeid,hits,frame,affine,qfinal,tcount) c ttemp=ttemp+tcount c if(mod(k,100).eq.0)then c print*, ' 100 Frame Bad Triangle Average: ',float(ttemp)/100. c ttemp=0 c endif C print*, ' Bad Triangles: ',tcount endif !!end if block to include only one orbit enddo !!end k frame loop close(12) print 1002, ss 1002 format('Orbit #',I2,' indexed. Calling average night') call averagenight_allsky(node,nodeid,hits,TOP_KMAX,level,image1,image2,image3, npolarimage, spolarimage) print*, 'Night Averaged, Writing Files...' if(ss.lt.1000)then !!Write image1(sky atlas) to file C write (outfile1,'(I4.4,A1,I2.2,A13)') date,'_',hour,'_EqPlot_1.grd' write (outfile1,'(I3.3,A7,I2.2,A8)') doy,'_Orbit_',ss,'_Eq1.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) C write (outfile2,'(I4.4,A1,I2.2,A13)') date,'_',hour,'_EqPlot_2.grd' write (outfile2,'(I3.3,A7,I2.2,A8)') doy,'_Orbit_',ss,'_Eq2.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) C write (outfile3,'(I4.4,A1,I2.2,A13)') date,'_',hour,'_EqPlot_3.grd' write (outfile3,'(I3.3,A7,I2.2,A8)') doy,'_Orbit_',ss,'_Eq3.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) C write (outfile4,'(I4.4,A1,I2.2,A11)') date,'_',hour,'_npolar.grd' write (outfile4,'(I3.3,A7,I2.2,A11)') doy,'_Orbit_',ss,'_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) C write (outfile4,'(I4.4,A1,I2.2,A11)') date,'_',hour,'_spolar.grd' write (outfile5,'(I3.3,A7,I2.2,A11)') doy,'_Orbit_',ss,'_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*, ' ' 98 continue print*, 'All Orbits Finished!!!' END