C+ C NAME: C smei_skyd_sky C PURPOSE: C Indexes a frame of SMEI data. C CATEGORY: C camera/for/sky C CALLING SEQUENCE: subroutine smei_skyd_sky(hdr,iframe,frame) C INPUTS: C hdr(*) double precision frame header C frame(*) real 2-D SMEI frame C OUTPUTS: C IMPLICIT: implicit none C INCLUDE: include 'filparts.h' include 'smei_frm_layout.h' include 'smei_frm_hdr.h' include 'smei_frm_hdrok.h' include 'smei_skyd_dim.h' C CALLS: C smei_hdr_quaternion, quaternion_rotate_xyz C smei_ccd2cam, smei_hdr_time, smei_orbit2, smei_FitPlane, Say C Time2Str, Str2Str, Int2Str, Flt2Str, Dbl2Str C iSetFileSpec, iGetFileSpec, iOSDeleteFile, cInt2Str, ForeignArgSet C smei_skyd_go, iFilePath, itrim, iGetLun, iFreeLun C BadI4, BadR4, BadR8, quaternion_inverse, smei_frm_c3mask_active C SEE ALSO: C smei_skyd_init, smei_skyd_make C PROCEDURE: C The photometric measurements from a single SMEI camera are combined C and then averaged to produce a surface-brightness sky map. C The "level" parameter governs the C granularity of the triangular coordinate frame. The coordinates of C the four corners of a pixel (or bin) are transformed down into this C frame, and the photometric surface brightness as measured by that C pixel is added to all nodes whose centers lie within the rectangle. C This surface brightness is corrected by 1/cosine of the pixels angle C to the camera axis, and by r/rad0 the pixels fractional distance to C the FOVs rotational center. When the data sequence is finished, all C triangle response sums are averaged and the resulting sky map binned C into an output surface brightness map. C MODIFICATION HISTORY: C 2002,2003 Aaron Smith (UCSD/CASS) C 2002->Original C 2003->removed most C++ functions and re-wrote in fortran as subroutines C ->add quaternion capabilities C DEC-2002, Paul Hick (UCSD/CASS) C Defined LEVEL and NMAX as parameters. NMAX is determined from LEVEL. C NMAX 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. C Rendered names compliant with those of B.V. Jackson C Incorporated 1/cosine surface-brightness correction for aperture C perspective (currently in smei_skyd_sky routine) C Incorporated large and small scale flatfield corrections C Incorporated FORTRAN read subroutine C Jan-May, 2003 Aaron Smith (UCSD) C note: not compliant with TMO data anymore C FEB-2005, Paul Hick (UCSD/CASS) C Reduced resolution of polar maps from 1600x1600 to 800x800. C Added theta-x and theta-y maps for glare removal as C low-resolution maps to Fits skymap. C NOV-2005, Paul Hick (UCSD/CASS) C Added another low resolution map to output file (fraction C of orbit in seconds since start of orbit). C The lowres maps for the 'glare angles' are now always calculated. C They used to be set zero if glare subtraction was suppressed. C All lowres maps for camera 2 are run through GridSphere2D to C fill holes near the equatorial poles. C DEC-2005, Paul Hick (UCSD/CASS) C Modified to deal with presence of -overwrite on cmd line C JUN-2006, Paul Hick (UCSD/CASS), V5.1 C Added smei_skyd_size subroutine to control size of C votes arrays. C Added direction cosine angle in long dimension to C skymap file. Note that this changes the layout of the C Fits skymap file. C JUL-2006, Paul Hick (UCSD/CASS), V5.3 C Fixed error in some of the comments for the map description C keyword for the polar maps. C Added Fits keyword STIME and ETIME to extension containing the time C STIME in seconds since the first frame used (these were only specified C in the main header). STIME is the time origin for the array. C MAR-2007, Paul Hick (UCSD/CASS; pphick@ucsd.edu), V6.1 C Added check for "bad pixel" mask presence for camera 3. C The "closed shutter" calibration pattern for camera 3 is different C depending on whether a "bad pixel" mask is in effect. C MAY-2007, Paul Hick (UCSD/CASS; pphick@ucsd.edu), V6.2 C Added GridSphereWeight calls to smooth low resolution maps C where only a few CCD pixels (< 3) contribut a skybin. C- double precision hdr (*) integer iframe real frame(*) character cSaySky * 8 /'skyd_sky' / character cSayInit * 9 /'skyd_init' / character cSayUpdate*11 /'skyd_update'/ character cSayFts * 8 /'skyd_fts' / character cSayMake * 9 /'skyd_make' / character cSeverity /'E'/ ! Arg list of smei_skyd_fts integer ibreak !================================= ! Frame arrays. ! Initialized in smei_skyd_init. ! Filled in smei_skyd_sky. ! Written to disk in smei_skyd_fts. integer badpix(SMEI__FRM_NPIX) save badpix real lsff(SMEI__FRM_NPIX) save lsff real cal_pattern_dark save cal_pattern_dark real cal_pattern(SMEI__FRM_NPIX) save cal_pattern real orb_pattern(SMEI__FRM_NPIX) save orb_pattern real pattern(SMEI__FRM_NPIX) save pattern double precision cal_version /0.0d0/ save cal_version double precision orb_version /0.0d0/ save orb_version !================================= ! High-resolution skymap arrays ! Filled in smei_skyd_make from indexing arrays and written to disk ! in smei_skyd_fts real img(SKYD__N_ALL) save img real xmg(SKYD__N_ALL) save xmg integer mask(SKYD__N_ALL) ! Star mask to identify locations close to stars save mask integer*1 jmg(SKYD__N_ALL) ! Used to store node counts in Fits file !================================= ! Low-resolution skymap arrays ! String ID for Fits file character cFitsID(SKYD__LO_NMAP)*35 & /'Dirty sky brightness', & 'Time (fraction of orbit)', & 'PSF position angle from north', & 'PSF position angle from equinox', & 'Discarded response map', & 'Discarded HTM nodes', & 'Theta-x for glare removal', & 'Theta-y for glare removal', & 'Time (sec since TORIGIN)', & 'Direction cosine angle-x', & 'Contributing pixel count', & ''/ ! Not written to Fits file character cPole(2)*5 /'north','south'/ character cSign(2) /'+' ,'-' / double precision dSign(2) / +1.0d0, -1.0d0/ real lowres(SKYD__NXLO,SKYD__NYLO,SKYD__LO_NMAP) save lowres real lowval(SKYD__LO_NMAP) !================================= ! smei_skyd_init: double precision hdrok(*) double precision bodyexcl(SMEI__HDROK_NBOX,SMEI__HDROK_NBODY) logical bGo integer icurrent_orbit integer TLow (2) integer THigh(2) logical bOverwrite save bOverwrite logical bExists save bExists logical bCheckVersion save bCheckVersion integer iSilent save iSilent integer iType save iType logical bKeepGlitch save bKeepGlitch ! sets : icam, mode, level, cOut, iOut (saved internally) ! initializes: 6 low resolution sky arrays character cVar(*)*(*) character cArg*(*) character cSkyType*7 /'.fts.gz'/ save cSkyType integer thisorbit /-1/ ! Integer orbit nr for internal use save thisorbit double precision thisperiod save thisperiod character cSplit*(FIL__LENGTH)! Build from cVar(4) and thisorbit save cSplit ! Template for name of output files integer iSplit ! Loc in cSplit where file name id is inserted save iSplit integer jSplit ! Effective length of template save jSplit integer level ! Extracted from cmd line cArg save level integer icam ! Copy of icam_ save icam integer mode ! Copy of mode_ save mode integer nbin ! Derived from mode save nbin ! .. binning (1,2,4) integer nX ! Derived from mode save nX ! .. hor frame size (1272,636,318) integer nY ! Derived from mode save nY ! .. vert frame size (256,128,64) integer nfov ! Derived from first SMEI frame save nfov ! .. # pixels inside of SMEI fov integer pfov(SKYD__NFOV) ! Filled from first SMEI frame save pfov ! .. linear pixel indices of nfov pixels inside SMEI fov real geometric(SKYD__NFOV) save geometric double precision runit(3,0:4,SKYD__NFOV) save runit double precision thetax(SKYD__NFOV) save thetax real badvalue save badvalue ! The biggest glare map is a mode 1 frame, but don't change this!!! real glare_map(SMEI__FRM_NPIX) /SMEI__FRM_NPIX*0.0/ save glare_map real glare_factor(SKYD__NX_GLARE,SKYD__NY_GLARE) save glare_factor real cx_glare save cx_glare real cy_glare save cy_glare real dx_glare save dx_glare real dy_glare save dy_glare double precision ra_min save ra_min double precision ra_max save ra_max double precision dec_min save dec_min double precision dec_max save dec_max integer tlowltd(2) save tlowltd integer thighltd(2) save thighltd integer torigin(2) ! Time origin for lowres time map save torigin integer tfirstfrm(2) ! Time of first frame used save tfirstfrm integer tlastfrm(2) ! Time of last frame used save tlastfrm integer tcurrentfrm(2) ! Time of current frame save tcurrentfrm character currentname*(FIL__LENGTH) save currentname ! Name of current frame double precision qq_equ2cam_edge(4) save qq_equ2cam_edge integer nxeq_ integer nyeq_ integer n_eq_ double precision d_eq_ double precision cxeq_ double precision cyeq_ integer nxpl_(2) integer nypl_(2) integer n_pl_(2) double precision d_pl_(2) double precision cxpl_(2) double precision cypl_(2) ! These six are accessed only in smei_skyd_init real zcut(0:2) /111.00, 27.76, 6.94/ !500 electron cut integer near_star integer nhot_pix real zz(-1:1,-1:1) real z0 real rcx real rcy double precision dfov(4) / 0.0d0, 0.0d0, 0.0d0, 0.0d0/ double precision dx(4) /-0.5d0,-0.5d0, 0.5d0, 0.5d0/ !Lower-left,Upper-left,upper-right,lower-right double precision dy(4) /-0.5d0, 0.5d0, 0.5d0,-0.5d0/ character cStr*(FIL__LENGTH) character cTmp*(FIL__LENGTH) double precision x(0:4) double precision y(0:4) double precision dr double precision phi double precision rr double precision tx double precision ty double precision rx double precision ry double precision rz double precision px double precision py double precision pz double precision dval real rval real darkratio integer i integer j integer k integer ii integer jj integer ibeg integer iend integer n integer norb integer nneg integer nbad integer ip integer iU integer istat character cInt2Str*14 logical ForeignArgSet logical smei_skyd_go logical smei_inside_fov logical smei_skyd_checkpix integer smei_frm_c3mask_active integer ibset integer Str2Str integer Int2Str integer Flt2Str integer Dbl2Str integer Time2Str integer iFilePath integer itrim integer iGetLun integer iFreeLun integer iSetFileSpec integer iGetFileSpec integer iOSDeleteFile integer iSearch integer Time2Differ integer BadI4 real BadR4 double precision BadR8 double precision dsind ! Needed for GNU compiler double precision dcosd double precision dasind double precision datan2d logical bSimple /.TRUE./ ! Conforms to FITS standard integer iBitPix integer nAxis / 2/ ! Two dim array integer nAxes(2) logical bExtend /.TRUE./ ! Allows FITS extensions real rbad double precision qq_cam2equ(4) double precision qq_equ2cam(4) double precision qq_cam2cam_edge(4) real ped real dark real fudge logical bOneSky save bOneSky logical bKeepGlare save bKeepGlare integer c3mask save c3mask logical bInRADecOff save bInRaDecOff logical bInRADec logical bglare logical bNewCal logical bNewOrb real glare_multiplier real response1 real response2 integer zero_gain_time(2) /1095,0/ ! 2002-DEC-31 0 UT double precision inverse_zero_gain(3) /0.97,1.00,0.92/ double precision sky_version double precision forbit_ext(2) double precision forbit_ltd(2) integer time (2) integer time_(2) integer hms (4) integer nframe integer ip_pix include 'smei_skyd_fnc.h' ip_pix(i,j) = (j-1)*nX+i ! Used to index SMEI frame arrays call smei_hdr_quaternion(hdr,1,qq_equ2cam) ! Equatorial -> camera frame call quaternion_inverse(qq_equ2cam,qq_cam2equ) ! Camera -> equatorial frame call quaternion_multiply(qq_cam2equ,qq_equ2cam_edge,qq_cam2cam_edge) ped = sngl( hdr(SMEI__HDR_PEDESTAL ) ) dark = sngl( hdr(SMEI__HDR_DARK_MEDIAN) ) ! Set up timing information for current frame: ! lowval(SKYD__LO_TIME) is a fraction of an orbit ! lowval(SKYD__LO_SECS) is the number of seconds since torigin call smei_orbit2(tcurrentfrm,i,dval) ! Current frame time (set in smei_skyd_init) dval = dble(i-thisorbit)+dval ! i could be thisorbit plus/minus 1 lowval(SKYD__LO_TIME) = sngl(dval) ! Fraction in current orbit call Time2Delta(tcurrentfrm,torigin,time) ! Time diff with torigin call Time2HMS(0,time,hms) ! Seconds since torigin lowval(SKYD__LO_SECS) = float(time(1)*86400+(hms(1)*60+hms(2))*60+hms(3)) ! time(1) might be -1; hms(4) should always be zero ! Set up glare subtraction for this frame rx = hdr( SMEI__HDR_SUN+0 ) ry = hdr( SMEI__HDR_SUN+1 ) rz = hdr( SMEI__HDR_SUN+2 ) lowval(SKYD__LO_ROTX) = sngl(datan2d(rx,rz)) lowval(SKYD__LO_ROTY) = sngl(datan2d(ry,rz)) glare_multiplier = 0.0 if (.not. bKeepGlare .and. rz .gt. dsind(1d0)) then! More than 1 deg above horizon ii = nint(cx_glare+dx_glare* lowval(SKYD__LO_ROTX) ) jj = nint(cy_glare+dy_glare*abs(lowval(SKYD__LO_ROTY))) if (1 .le. ii .and. ii .le. SKYD__NX_GLARE .and. 1 .le. jj .and. jj .le. SKYD__NY_GLARE) then glare_multiplier = glare_factor(ii,jj) !write (*,'(A,2F6.1,2I4,F5.1)') currentname(:itrim(currentname)), lowval(SKYD__LO_ROTX), lowval(SKYD__LO_ROTY), ii, jj, glare_multiplier end if end if bglare = glare_multiplier .ne. 0.0 darkratio = dark/cal_pattern_dark norb = 0 nneg = 0 nbad = 0 if (iSilent .le. 2) call Say(cSaySky,'I','','.#frame '//currentname) ! Remove nodes that have moved outside the FOV call smei_skyd_flush(iType,bKeepGlitch,iSilent,level,icam,mode,qq_equ2cam,SKYD__DRFOV_MAX,img,xmg,badpix,lowres) call smei_skyd_sort(iSilent) ! Sort pool call smei_skyd_save() do n=1,nfov ! Loop over all pixels inside fov ! If tcurrentfrm < tlowltd or tcurrentfrm > thighltd ! we need to check whether or not the pixel is inside the ! orbit. Rotate (rx,ry,rz) from the current camera frame ! to the camera 2 frame at the orbit start. ! The sign of the y-component decides whether the pixel ! is inside the orbit or not. ii = 0 if (bOneSky) then ii = 1 else if (Time2Differ(tcurrentfrm,tlowltd) .lt. 0) then ! Start of orbit if ( smei_skyd_checkpix(pfov(n),runit(1,0,n),qq_cam2cam_edge)) ii = 1 else if (Time2Differ(tcurrentfrm,thighltd) .gt. 0) then! End of orbit if (.not. smei_skyd_checkpix(pfov(n),runit(1,0,n),qq_cam2cam_edge)) ii = 1 else ii = 1 end if ip = pfov(n) ! Linear index ! Note that response pretty much always should be positive ! for pixels inside the fov, so this test could be redundant. if (ii .eq. 1 .and. frame(ip) .gt. 0.0) then norb = norb+1 rx = runit(1,0,n) ! Unit vector for pixel in ry = runit(2,0,n) ! .. camera coordinate frame rz = runit(3,0,n) dval = dsqrt(1.0d0-rx*rx) px = 0.0d0 ! Unit vector for PSF orientation py = rz/dval ! .. in camera coordinate frame pz = -ry/dval ! Camera to equatorial frame call quaternion_rotate_xyz(qq_cam2equ,rx,ry,rz) ! Force (rx,ry,rz) to be a unit vector. This is primarily done to avoid ! problems with asind below for vectors close to one of the poles. dval = dsqrt(rx*rx+ry*ry+rz*rz) rx = rx/dval ry = ry/dval rz = rz/dval call quaternion_rotate_xyz(qq_cam2equ,px,py,pz) dval = dsqrt(px*px+py*py+pz*pz) px = px/dval py = py/dval pz = pz/dval ! psfn and psfe are 'position angles' counterclockwise from ! the reference direction (as viewed from outside the unit ! sphere) with equatorial north and the vernal equinox respectively. lowval(SKYD__LO_PSFN) = sngl( datan2d( ry*px-rx*py,pz) ) lowval(SKYD__LO_PSFE) = sngl( datan2d(-rz*px+rx*pz,py) ) k = 0 ! Center of pixel x(k) = datan2d(ry,rx) ! Right ascension x(k) = dmod(x(k)+360.0d0,360.0d0) ! Force into [0,360) rz = max(-1.0d0,min(rz,1.0d0)) ! Safety belt for dasind fnc y(k) = dasind(rz) ! Declination ! Check whether pixel is in selected piece of sky bInRADec = bInRADecOff if (.not. bInRADecOff) then bInRADec = dec_min .le. y(k) .and. y(k) .le. dec_max if (bInRADec) bInRADec = & (ra_min .le. ra_max .and. (ra_min .le. x(k) .and. x(k) .le. ra_max)) .or. & (ra_min .gt. ra_max .and. (ra_min .le. x(k) .or. x(k) .le. ra_max)) end if if (bInRADec) then do k=1,4 ! Corners (lower-right, upper-left, upper-right, lower-right) rx = runit(1,k,n) ry = runit(2,k,n) rz = runit(3,k,n) ! Camera --> equatorial frame call quaternion_rotate_xyz(qq_cam2equ,rx,ry,rz) dval = dsqrt(rx*rx+ry*ry+rz*rz) rx = rx/dval ry = ry/dval rz = rz/dval x(k) = datan2d(ry,rx) ! Right ascension x(k) = dmod(x(k)+360.0d0,360.0d0)! Force into [0,360) rz = max(-1.0d0,min(rz,1.0d0))! Safety belt for dasind fnc y(k) = dasind(rz) ! Declination end do ! end vertex loop i = 1+mod(ip-1,nX) ! Column index j = 1+(ip-1)/nX ! Row index lowval(SKYD__LO_SKY) = frame(ip) ! Take away pedestal and dark-current-scaled pattern. ! Apply flat field correction. response1 = lowval(SKYD__LO_SKY) lowval(SKYD__LO_SKY) = (lowval(SKYD__LO_SKY)-(ped+pattern(ip)*darkratio))/lsff(ip) response2 = lowval(SKYD__LO_SKY) if (bglare) then ! Take out the glare !print *, lowval(SKYD__LO_SKY), i,j lowval(SKYD__LO_SKY) = lowval(SKYD__LO_SKY)-glare_multiplier*glare_map(ip) !print *, lowval(SKYD__LO_SKY) end if ! Print a warning if the response is negative at this point. Note that ! we haven't made the geometric corrections yet, but that involves just a ! multiplication that can't make the response negative. if (iSilent .le. 0) then if (lowval(SKYD__LO_SKY) .lt. 0) then nneg = nneg+1 if (nneg .lt. 20) then ii = 0 ii = ii+Str2Str('I<0 [' ,cStr(ii+1:)) ii = ii+Int2Str(i ,cStr(ii+1:)) ii = ii+Str2Str(',' ,cStr(ii+1:)) ii = ii+Int2Str(j ,cStr(ii+1:)) ii = ii+Str2Str(']:' ,cStr(ii+1:))+1 ii = ii+Flt2Str(response1,-2 ,cStr(ii+1:))+1 ii = ii+Flt2Str(response2,-2 ,cStr(ii+1:))+1 ii = ii+Flt2Str(lowval(SKYD__LO_SKY) ,-2,cStr(ii+1:)) else if (nneg .eq. 20) then cStr = '.... and more' else cStr = ' ' end if if (cStr .ne. ' ') call Say(cSaySky,'W',currentname,cStr) end if end if ! Apply geometric corrections lowval(SKYD__LO_SKY ) = lowval(SKYD__LO_SKY)*geometric(n) lowval(SKYD__LO_HIT ) = 1.0 lowval(SKYD__LO_FOVX) = sngl(thetax(n)) ! Before sending the pixel to the indexing routine determine whether ! or not it is close to a bright star. jj = nint(origin(1,SKYD__NYEQ)+y(0)*SKYD__D_EQ) if (1 .le. jj .and. jj .le. SKYD__NYEQ) then! Latitude inside equatorial map ii = nint(origin(1,0)+x(0)*SKYD__D_EQ)! (ii,jj) is nearest bin in equat map if (ii .lt. 1 .or. ii .gt. SKYD__NXEQ) then! Better not happen! write (*,*) ' RA=',x(0),' Dec=',y(0) stop 'smei_skyd_sky: bad RA, darn' end if ii = ip_eq(1,ii,jj) near_star = mask(ii) ! 0=away from star; 1=near star else ! Use polar map ! The equatorial map and polar maps overlap, so if a point is not in the ! equatorial map, it must be in one of the polar maps. dval = SKYD__D_PL*(90.0d0-dabs(y(0))) ii = nint(origin(1,SKYD__NXPL)+dval*dcosd(x(0))) jj = nint(origin(1,SKYD__NYPL)+dval*dsind(x(0))) ! (ii,jj) is nearest bin in polar map if (ii .lt. 1 .or. ii .gt. SKYD__NXPL .or. jj .lt. 1 .or. jj .gt. SKYD__NYPL) then write (*,*) ' RA=',x(0),' Dec=',y(0)! Better not happen! stop 'smei_skyd_sky: bad dec, darn' end if ii = ip_pl(1,ii,jj,(3-sign(1,int(y(0))))/2)! 1 for North, 2 for South near_star = mask(ii) ! 0=away from star; 1=near star end if nhot_pix = 0 if (icam .eq. 3 .or. near_star .eq. 1) then do jj=-1,1 do ii=-1,1 k = ip_pix(i+ii,j+jj) zz(ii,jj) = (frame(k)-(ped+pattern(k)*darkratio))/lsff(k) end do end do call smei_FitPlane(zz,zz,z0,rcx,rcy,rval) ! If close to bright star, then check slope. If it is less than ! the threshold slope SKYD__MAX_SLOPE, switch off the near_star flag if (near_star .eq. 1) then if (sqrt(rcx*rcx+rcy*rcy) .le. SKYD__MAX_SLOPE) near_star = 0 end if if (icam .eq. 3) then if (abs(zz(0,0)) .gt. zcut(mode)) then !if((z(0,0)*z(0,0))/rval.gt.5.)then !write (*,'(A,F8.2,A,F10.2,2(A,F7.2),A,F10.2,2(A,I3),A) !& 'response = ',lowval(SKYD__LO_SKY),' sumzq = ',rval,' zz(0,0) = ',zz(0,0),' z0 = ',z0, !& ' zz(0,0)**2/sumzsqr = ',zz(0,0)*zz(0,0)/rval,' (i,j) = (',i,',',j,')' !endif ! Remove dark current slope from hot/flipper pixel flag if (zz(0,0)*zz(0,0) .gt. 25.0*rval .or. (pattern(ip)-cal_pattern_dark*0.4*(1.0+0.01*j) & .gt. 50.0 .and. zz(0,0) .lt. 0.0)) nhot_pix = 1 end if end if end if ii = 0 if (near_star .eq. 1) ii = ibset(ii,0) if (nhot_pix .eq. 1) ii = ibset(ii,1) call smei_skyd_fill(iSilent,level,icam,mode,x,y,ip,lowval(SKYD__LO_SKY),ii) ! Determine in which low-res bin to drop the current pixel. ! 0 <= x < 360 -> 1 <= ii <= NXLO , i.e. ii is always valid ! -90 <= y <= +90 -> 1 <= jj <= NYLO+1, i.e. jj is invalid for y=+90 ii = nint(origin(1, 0)+x(0)*SKYD__D_LO)! Pixel center jj = nint(origin(1,SKYD__NYLO)+y(0)*SKYD__D_LO) jj = min(jj,SKYD__NYLO) ! Only used if y(0) exactly +90.0d0 if (ii .lt. 1 .or. ii .gt. SKYD__NXLO .or. jj .lt. 1 .or. jj .gt. SKYD__NYLO) then write (*,*) ' RA=',x(0),' Dec=',y(0) ! Better not happen! stop 'skyd_sky: bad low res bin, darn' end if ! Sum low resolution maps call smei_skyd_addvalue(lowval(SKYD__LO_FOVX),360.0 ,lowres(ii,jj,SKYD__LO_HIT),lowres(ii,jj,SKYD__LO_FOVX))! Direction cosine angle call smei_skyd_addvalue(lowval(SKYD__LO_PSFN),360.0 ,lowres(ii,jj,SKYD__LO_HIT),lowres(ii,jj,SKYD__LO_PSFN))! PSF rotation call smei_skyd_addvalue(lowval(SKYD__LO_PSFE),360.0 ,lowres(ii,jj,SKYD__LO_HIT),lowres(ii,jj,SKYD__LO_PSFE)) call smei_skyd_addvalue(lowval(SKYD__LO_ROTX),360.0 ,lowres(ii,jj,SKYD__LO_HIT),lowres(ii,jj,SKYD__LO_ROTX))! Glare angles call smei_skyd_addvalue(lowval(SKYD__LO_ROTY),360.0 ,lowres(ii,jj,SKYD__LO_HIT),lowres(ii,jj,SKYD__LO_ROTY)) call smei_skyd_addvalue(lowval(SKYD__LO_TIME), 1.0 ,lowres(ii,jj,SKYD__LO_HIT),lowres(ii,jj,SKYD__LO_TIME))! Fractions call smei_skyd_addvalue(lowval(SKYD__LO_SECS),sngl(thisperiod),lowres(ii,jj,SKYD__LO_HIT),lowres(ii,jj,SKYD__LO_SECS))! Seconds lowres(ii,jj,SKYD__LO_SKY) = lowres(ii,jj,SKYD__LO_SKY)+lowval(SKYD__LO_SKY) ! Sidereal skymap at low resolution lowres(ii,jj,SKYD__LO_HIT) = lowres(ii,jj,SKYD__LO_HIT)+lowval(SKYD__LO_HIT) end if else if (ii .eq. 1) then norb = norb+1 nbad = nbad+1 end if end do ! Last frame time for frame with pixels inside orbit ! (saved for Fits header) if (norb .gt. 0) then if (iframe .eq. 0) then call ArrI4Copy(2,tcurrentfrm,tfirstfrm) call Say(cSaySky,'I',currentname,'starts orbit '//cInt2Str(thisorbit)) end if iframe = iframe+1 call ArrI4Copy(2,tcurrentfrm,tlastfrm) end if if (iSilent .le. 0) call smei_skyd_print() if (iSilent .le. 1) then if (nneg .ne. 0) then i = 0 i = i+Str2Str('intensity in', cStr(i+1:))+1 i = i+Int2Str(nneg , cStr(i+1:)) i = i+Str2Str('/' , cStr(i+1:)) i = i+Int2Str(norb , cStr(i+1:)) i = i+Str2Str(' pixels' , cStr(i+1:)) call Say(cSaySky,'W','negative',cStr) end if if (nbad .ne. 0) then i = 0 i = i+Int2Str(nbad , cStr(i+1:)) i = i+Str2Str('/' , cStr(i+1:)) i = i+Int2Str(norb , cStr(i+1:)) i = i+Str2Str(' bad pixels', cStr(i+1:)) call Say(cSaySky,'W','ignored',cStr) end if end if return C+ C NAME: C smei_skyd_init C PURPOSE: C Initialize variables for smei_skyd_sky. C CALLING SEQUENCE: entry smei_skyd_init(sky_version,icurrent_orbit,TLow,THigh,hdr,cVar,cArg,bGo) C INPUTS: C sky_version double precision smei_skyd version number C icurrent_orbit integer orbit number of the current orbit C being processed (the 'hdr' for the C the current frame might still belong C to the previous of next orbit). C hdr(*) double precision frame header array C cVar(*) character*(*) cVar(2): pattern file or SMEIDB? C cVar(4): output directory C cArg character*(*) command line keywords C OUTPUTS: C bGo logical bGo is returned .FALSE. only C if the first frame for a new skymap C is processed, and the Fits file for C the skymap already exists. C CALLS: C ForeignI4Arg, smei_orbit2, smei_orbit_time2, Str2Str, Int2Str, Time2Str C iFilePath, ArrR4Zero, ArrI8Zero, ArrI2Zero, ArrI4Zero, Say, cInt2Str C smei_cal_get, smei_get_glare, smei_get_lsff, smei_get_starmask C ForeignArgSet, smei_hdr_time C PROCEDURE: C Entry point in href=smei_skyd_sky=. C MODIFICATION HISTORY: C NOV-2004, Paul Hick (UCSD/CASS) C DEC-2005, Paul Hick (UCSD/CASS) C Modified to deal with presence of -overwrite on cmd line C JUL-2007, Paul Hick (UCSD/CASS; pphick@ucsd.edu) C Added check for cmd line argument -c3mask=0/1 C Added argument ipick to smei_cal_get calls C- bGo = .TRUE. call smei_hdr_time(hdr,tcurrentfrm) ! Frame time, saved for Fits header if (thisorbit .eq. -1) then ! First frame of first orbit ! If bOneSky is set then this section is ! .. only passed once. bOneSky = ForeignArgSet(cArg,'onesky') ! First frame of first orbit: do stuff that needs to be done only once icam = nint(hdr(SMEI__HDR_CAMERA)) ! Save: camera mode = nint(hdr(SMEI__HDR_MODE )) ! Save: mode nbin = 2**mode ! Save: binning call ForeignI4Arg(cArg,'level',-1,level) if (level .eq. -1) then ! No level specified level = 3 ! Science modes @ level 11 if (mode .eq. 0) level = 5 ! Engineering mode @ level 12 end if if (level .gt. SKYD__LEVEL) call Say(cSayInit,'E','level '//cInt2Str(level), & 'is too high; max is SKYD__LEVEL='//cInt2Str(SKYD__LEVEL)) level = max(1,min(level,SKYD__LEVEL))! Save: level call ForeignI4Arg(cArg,'silent',0,iSilent) bKeepGlare = ForeignArgSet(cArg,'keepglare' ) bKeepGlitch = ForeignArgSet(cArg,'keepglitch') if (ForeignArgSet(cArg,'mean')) then iType = 0 else if (ForeignArgSet(cArg,'median')) then iType = 1 else call ForeignI4Arg(cArg,'type',0,iType) end if if (icam .eq. 3 .and. mode .eq. 1) then call ForeignI4Arg(cArg,'c3mask',-1,c3mask) if (c3mask .ne. 0 .and. c3mask .ne. 1) c3mask = -1 else c3mask = 0 end if ! Limit construction of skymap to box of sky in RA and dec call ForeignR8Arg(cArg,'ra_min' , -1.0, ra_min ) call ForeignR8Arg(cArg,'ra_max' , -1.0, ra_max ) call ForeignR8Arg(cArg,'dec_min',-91.0, dec_min) call ForeignR8Arg(cArg,'dec_max',-91.0, dec_max) ! All four must be set. If one is omitted the others are silently ignored bInRADecOff = ra_min .lt. 0.0 .or. ra_max .lt. 0.0 .or. & dec_min .lt. -90.0 .or. dec_max .lt. -90.0 if (.not. bInRADecOff) then ! Adjust input values making sure they are located at pixel ! edges (at half-integer array index values. ra_min = mod(mod(ra_min,360.0)+360.0,360.0) ra_max = mod(mod(ra_max,360.0)+360.0,360.0) dval = origin(1,0) ra_min = (nint(dval+ra_min*SKYD__D_EQ)-0.5-dval)/SKYD__D_EQ ra_max = (nint(dval+ra_max*SKYD__D_EQ)+0.5-dval)/SKYD__D_EQ dval = origin(1,SKYD__NYEQ) dec_min = (nint(dval+dec_min*SKYD__D_EQ)-0.5-dval)/SKYD__D_EQ dec_max = (nint(dval+dec_max*SKYD__D_EQ)+0.5-dval)/SKYD__D_EQ dec_min = max(-90.0,min(dec_min,90.0)) dec_max = max(-90.0,min(dec_max,90.0)) i = 0 i = i+Str2Str('RA range (deg):' , cStr(i+1:))+1 i = i+Dbl2Str(ra_min, 3 , cStr(i+1:)) i = i+Str2Str(' to' , cStr(i+1:))+1 i = i+Dbl2Str(ra_max, 3 , cStr(i+1:)) i = i+Str2Str(', Dec range (deg)', cStr(i+1:))+1 i = i+Dbl2Str(dec_min, 3 , cStr(i+1:)) i = i+Str2Str(' to' , cStr(i+1:))+1 i = i+Dbl2Str(dec_max, 3 , cStr(i+1:)) call Say(cSayInit,'I','Limit',cStr) end if if (bKeepGlare) then call Say(cSayInit,'I','glare','subtraction skipped') else bKeepGlare = icam .eq. 1 .or. (icam .eq. 2 .and. mode .eq. 1) if (bKeepGlare) then call Say(cSayInit,'I','glare','removed only for c2m0, c2m2 and c3') else ! glare_map is a 318x 64 array for c2m2 and 636x128 for c3m1 ! glare_factor is a 90x180 array call smei_get_glare(icam,mode,glare_map,glare_factor,cx_glare,cy_glare,dx_glare,dy_glare) end if end if badvalue = 0.0 ! Save: bad value flag (put in Fits header) nX = SMEI__FRM_NX/nbin ! Save nY = SMEI__FRM_NY/nbin ! Save ! Get the star masks (0 = away from star; 1 is near star) call smei_get_starmask(SKYD__N_EQ,SKYD__N_PL,mask) ! Get large-scale flatfield correction. Zeroes are set to one. call smei_get_lsff(mode,icam,lsff) ! Run through the first frame to make a list of all pixels inside the fov ! that need to be processed. ibeg = (SMEI__FRM_LEFT-1)/nbin+1 ! Partially covered column iend = SMEI__FRM_RIGHT /nbin if (mode .eq. 0) then ! For mode 0 drop the half-covered column on either side ibeg = ibeg+1 iend = iend-1 end if nfov = 0 do j=1,nY do i=ibeg,iend rx = dble(i) ry = dble(j) ! Convert from ccd coordinates to sky coordinates in the camera ! coordinate frame (with z along the optical axis and x along the long ! dimension of the fov). call smei_axis_ccd(2,icam,mode,rx,ry,dr,phi,rr) ! dr : radial distance from optical axis (in mode 0 pixels) ! phi: angle from optical axis ! rr : (radial distance from origin)/(radial dist of opt axis) if (smei_inside_fov(icam,mode,dfov,dr,phi)) then ! Pixel inside FOV ? k = ip_pix(i,j) ! Linear array index ! Exclude pixels along the edge. For mode 1 and 2 these could ! be bins including the partially covered pixels on left or right. ! For camera 3 this also avoids problems when filling the ! zz(-1:1,-1:1) array. None of these should be inside the fov, ! so the check is probably redundant. if (i .ne. ibeg .and. i .ne. iend .and. j .ne. 1 .and. j .ne. nY) then nfov = nfov+1 ! Save: number of pixels inside SMEI fov if (nfov .gt. SKYD__NFOV) call Say(cSayInit,'E','SKYD__NFOV='//cInt2Str(SKYD__NFOV),'is too small') pfov(nfov) = k ! Save: 1-based linear array index into SMEI frame call smei_axis_cam(0,icam,mode,dr,phi,runit(1,0,nfov),runit(2,0,nfov),runit(3,0,nfov),thetax(nfov),ty) ! Geometric corrections to response using pixel center position. ! - First order radial binsize correction ! - Radial binsize, nonlinear term ! - Correct surface brightness for aperture effective area geometric(nfov) = sngl(rr*(1d0-0.0092d0*ty)/runit(3,0,nfov)) ! Unit vectors in camera frame for all four corners of the pixel do n=1,4 call smei_axis_ccd(0,icam,mode,rx+dx(n),ry+dy(n),dr,phi,rr) call smei_axis_cam(0,icam,mode,dr,phi,runit(1,n,nfov),runit(2,n,nfov),runit(3,n,nfov),tx,ty) end do else k = 0 k = k+Str2Str('[', cStr(k+1:)) k = k+Int2Str( i , cStr(k+1:)) k = k+Str2Str(',', cStr(k+1:)) k = k+Int2Str( j , cStr(k+1:)) k = k+Str2Str(']', cStr(k+1:)) call Say(cSayInit,'W',cStr,'ditched edge pixel inside fov') end if end if end do end do call Say(cSayInit,'I',cInt2Str(nfov),'pixels inside the FOV') ! Build file name template for output files. ! cVar(4) stores the output directory. if (bOneSky) then ! Store time of first good frame ! (used in file name construction, is stored in Fits header ! and is used as time origin for the orbit time (in secs) call ArrI4Copy(2,tcurrentfrm,torigin) i = 0 i = i+Str2Str ('c' , cStr(i+1:)) i = i+Int2Str (icam, cStr(i+1:)) i = i+Str2Str ('_' , cStr(i+1:)) iSplit = i ! Save i = i+Time2Str(SMEI__UT_FORMAT,torigin, cStr(i+1:)) ! cVar(4) is output directory jSplit = iFilePath(cVar(4),0,'',cStr,cSplit) !Save iSplit = iSplit+jSplit-i ! Save ! bOverwrite: .TRUE. if -overwrite is set on command line ! bCheckversion: .TRUE. if -checkversion is set on command line ! bGo : signals to main program whether to keep processing frames ! bExists: used in smei_frm_fts to delete a file if -overwrite is set bOverwrite = ForeignArgSet(cArg,'overwrite') bCheckVersion = ForeignArgSet(cArg,'checkversion') cStr = cSplit(:iSplit-1)//'sky'//cSplit(iSplit:jSplit)//cSkyType bGo = smei_skyd_go(bOverwrite,bCheckVersion,cStr,sky_version,bExists) call ArrR4Zero(SKYD__N_ALL,img) call ArrR4Zero(SKYD__N_ALL,xmg) call ArrI4Zero(SMEI__FRM_NPIX/(nbin*nbin),badpix) call ArrR4Zero(SKYD__N_LO*SKYD__LO_NMAP,lowres)! Initialize low resolution maps i = 0 i = i+Time2Str(SMEI__UT_FORMAT,torigin,cStr(i+1:)) i = i+Str2Str(', camera', cStr(i+1:))+1 i = i+Int2Str(icam , cStr(i+1:)) i = i+Str2Str(', mode' , cStr(i+1:))+1 i = i+Int2Str(mode , cStr(i+1:))+1 i = i+Str2Str('(' , cStr(i+1:)) i = i+Int2Str(nbin , cStr(i+1:)) i = i+Str2Str('x' , cStr(i+1:)) i = i+Int2Str(nbin , cStr(i+1:))+1 i = i+Str2Str('binned)' , cStr(i+1:))+1 i = i+Str2Str('@ level' , cStr(i+1:))+1 i = i+Int2Str(level , cStr(i+1:)) call Say(cSayInit,'I','start',cStr) call smei_orbit2(torigin,thisorbit,dval) ! Orbit number call smei_orbit_period2(thisorbit,0.0d0,thisperiod) ! Orbit period thisperiod = thisperiod*86400.0d0 call Time2smei_quaternion(torigin,2,1,qq_equ2cam_edge) end if end if if (.not. bOneSky) then ! icurrent_orbit is the loop variable used in the calling program ! to loop over all orbits. ! Note that the time stored in the frame hdr is probably still ! part of the previous orbit when icurrent_orbit changes if (icurrent_orbit .ne. thisorbit) then! Start new orbit call ArrI4Copy(2,TLow ,tlowltd ) call ArrI4Copy(2,THigh,thighltd) ! First frame of new orbit: do stuff that needs to be done ! for each orbit. thisorbit = icurrent_orbit ! New orbit number ! Orbit start time, rounded to the nearest second call smei_orbit_time2(thisorbit,0.0d0,torigin) torigin(2) = nint(dble(torigin(2))/1000.0d0)*1000 ! Orbital period call smei_orbit_period2(thisorbit,0.0d0,thisperiod) thisperiod = thisperiod*86400.0d0 ! Get quaternion that rotates from equatorial coordinates ! to camera 2 coordinates at the start of the orbit. ! This quaternion is used to draw a 'line on skey' defining ! the start/end of the orbit. ! Note that we use camera 2 to define this edge for all ! three cameras. This makes for a better match when the ! skymaps for the three cameras are put together. call Time2smei_quaternion(torigin,2,1,qq_equ2cam_edge) ! Save: "edge-of-the-map" quaternion ! For the first frame all pixels should lie outside orbit ! thisorbit. Check, and print a warning if this is not the case. ! (note that there are valid reasons for this to happen, i.e. ! a data gap near the beginning of the orbit). !call smei_hdr_quaternion(hdr,2,qq_cam2equ) ! Camera -> equatorial frame ! Get quaternion that rotates from the camera frame for the ! current CCD frame, to the camera frame for camera 2 at the ! start of the orbit. !call quaternion_multiply(qq_cam2equ,qq_equ2cam_edge,qq_cam2cam_edge) ! Build file name template for output files. ! cVar(4) stores the output directory. i = 0 i = i+Str2Str ('c' , cStr(i+1:)) i = i+Int2Str (icam, cStr(i+1:)) i = i+Str2Str ('_' , cStr(i+1:)) iSplit = i ! Save i = i+Time2Str(SMEI__UT_FORMAT,torigin,cStr(i+1:)) ! cVar(4) is output directory jSplit = iFilePath(cVar(4),0,'',cStr,cSplit) !Save iSplit = iSplit+jSplit-i ! Save ! bOverwrite: .TRUE. if -overwrite is set on command line ! bCheckversion: .TRUE. if -checkversion is set on command line ! bGo : signals to main program whether to keep processing frames ! bExists: used in smei_frm_fts to delete a file if -overwrite is set bOverwrite = ForeignArgSet(cArg,'overwrite') bCheckVersion = ForeignArgSet(cArg,'checkversion') cStr = cSplit(:iSplit-1)//'sky'//cSplit(iSplit:jSplit)//cSkyType bGo = smei_skyd_go(bOverwrite,bCheckVersion,cStr,sky_version,bExists) call ArrR4Zero(SKYD__N_ALL,img) call ArrR4Zero(SKYD__N_ALL,xmg) call ArrI4Zero(SMEI__FRM_NPIX/(nbin*nbin),badpix) call ArrR4Zero(SKYD__N_LO*SKYD__LO_NMAP,lowres)! Initialize low resolution maps i = 0 i = i+Str2Str('camera' , cStr(i+1:))+1 i = i+Int2Str(icam , cStr(i+1:))+1 i = i+Str2Str('in mode' , cStr(i+1:))+1 i = i+Int2Str(mode , cStr(i+1:))+1 i = i+Str2Str('(' , cStr(i+1:)) i = i+Int2Str(nbin , cStr(i+1:)) i = i+Str2Str('x' , cStr(i+1:)) i = i+Int2Str(nbin , cStr(i+1:))+1 i = i+Str2Str('binned)' , cStr(i+1:))+1 i = i+Str2Str('@ level' , cStr(i+1:))+1 i = i+Int2Str(level , cStr(i+1:)) call Say(cSayInit,'I','orbit '//cInt2Str(thisorbit),cStr) end if end if ! Finally, do stuff that needs to be done for each frame. ! The "closed shutter" calibration pattern is either stored in cVar(2) ! or if cVar(2) is 'SMEIDB?' it is selected from the SMEI database ! based on the time stored in tcurrentfrm. ! For camera 3 in mode 1, cVar(5) contains the orbital "on-the-fly" ! difference pattern to be added to the "closed shutter" pattern. if (icam .eq. 3 .and. mode .eq. 1) then if (c3mask .eq. -1) then i = Time2Str(SMEI__UT_FORMAT,tcurrentfrm,cStr) i = smei_frm_c3mask_active(0,cStr(:i), & nint(hdr(SMEI__HDR_FLAT_ENABLED)),cStr,fudge) else i = c3mask end if call smei_cal_get(cVar(2),0,icam,mode,i,tcurrentfrm,cStr,cal_pattern,cal_pattern_dark,bNewCal,dval) if (bNewCal) cal_version = dval call smei_orb_get(cVar(5),icam,mode,cStr,orb_pattern,bNewOrb,dval) if (bNewOrb) orb_version = dval if (bNewCal .or. bNewOrb) call ArrR4PlusArrR4(nX*nY,cal_pattern,orb_pattern,pattern) else call smei_cal_get(cVar(2),0,icam,mode,0,tcurrentfrm,cStr,pattern,cal_pattern_dark,bNewCal,dval) if (bNewCal) cal_version = dval end if if (nint(hdr(SMEI__HDR_CAMERA)) .ne. icam .or. & nint(hdr(SMEI__HDR_MODE )) .ne. mode) & call Say(cSayInit,'E','wrong','camera or mode') i = 0 i = i+Str2Str('c' ,currentname(i+1:)) i = i+Int2Str(icam,currentname(i+1:)) i = i+Str2Str('m' ,currentname(i+1:)) i = i+Int2Str(mode,currentname(i+1:)) i = i+Str2Str('_ ',currentname(i+1:)) i = i+Time2Str(SMEI__UT_FORMAT,tcurrentfrm,currentname(i+1:)) return C+ C NAME: C smei_skyd_make C PURPOSE: C Takes the three node arrays (node, nodeid, and hits), C averages them and puts them into a two-dimensional skymap C CATEGORY: C camera/for/sky C CALLING SEQUENCE: entry smei_skyd_make() C OUTPUTS: C Output goes to three equatorial .grd files, two polar .grd files C and a FTS file containing both equatorial and polar maps. C CALLS: C Say, ArrI4Zero, ArrR4Zero, TRIANGLECENTER, smei_skyd_node C ArrR4DivideByArrI4, ArrR4DivideByArrR4, ArrR4Mask C PROCEDURE: C > Entry point in href=smei_skyd_sky=. See that procedure for more information. C MODIFICATION HISTORY: C April 2003, Aaron Smith (UCSD/CASS) C Updated to include capability for all sky imaging in 3 equatorial C plots, and two polar plots C- call Say(cSayMake,'I',currentname,'ends orbit '//cInt2Str(thisorbit)) ! Put all remaining nodes on the skymaps rval = BadR4() call smei_skyd_flush(iType,bKeepGlitch,iSilent,level,icam,mode,qq_equ2cam,rval,img,xmg,badpix,lowres) call smei_skyd_sort(iSilent) ! Not really necessary (clears n_nodes entries) call smei_skyd_print_nodes(icam,mode) ! The division of xmg and img will set all bins with img=0 to BadR4(). ! Set the bad values to zero before writing to file. call ArrR4DivideByArrR4(SKYD__N_ALL,xmg,img,xmg) ! All low resolution maps. The division sets lowres elements where ! lowres(*,*,SKYD__LO_HIT)=0.0 to BadR4(). call ArrR4DivideByArrR4(SKYD__N_LO,lowres(1,1,SKYD__LO_SKY ),lowres(1,1,SKYD__LO_HIT),lowres(1,1,SKYD__LO_SKY )) call ArrR4DivideByArrR4(SKYD__N_LO,lowres(1,1,SKYD__LO_TIME),lowres(1,1,SKYD__LO_HIT),lowres(1,1,SKYD__LO_TIME)) call ArrR4DivideByArrR4(SKYD__N_LO,lowres(1,1,SKYD__LO_SECS),lowres(1,1,SKYD__LO_HIT),lowres(1,1,SKYD__LO_SECS)) call ArrR4DivideByArrR4(SKYD__N_LO,lowres(1,1,SKYD__LO_FOVX),lowres(1,1,SKYD__LO_HIT),lowres(1,1,SKYD__LO_FOVX)) call ArrR4DivideByArrR4(SKYD__N_LO,lowres(1,1,SKYD__LO_PSFN),lowres(1,1,SKYD__LO_HIT),lowres(1,1,SKYD__LO_PSFN)) call ArrR4DivideByArrR4(SKYD__N_LO,lowres(1,1,SKYD__LO_PSFE),lowres(1,1,SKYD__LO_HIT),lowres(1,1,SKYD__LO_PSFE)) call ArrR4DivideByArrR4(SKYD__N_LO,lowres(1,1,SKYD__LO_ROTX),lowres(1,1,SKYD__LO_HIT),lowres(1,1,SKYD__LO_ROTX)) call ArrR4DivideByArrR4(SKYD__N_LO,lowres(1,1,SKYD__LO_ROTY),lowres(1,1,SKYD__LO_HIT),lowres(1,1,SKYD__LO_ROTY)) call ArrR4DivideByArrR4(SKYD__N_LO,lowres(1,1,SKYD__LO_CUT ),lowres(1,1,SKYD__LO_CNT),lowres(1,1,SKYD__LO_CUT )) call ArrR4DivideByArrR4(SKYD__N_LO,lowres(1,1,SKYD__LO_FRAC),lowres(1,1,SKYD__LO_CNT),lowres(1,1,SKYD__LO_FRAC)) ! The low res maps have small areas near the poles which do not get filled, ! i.e. lowres(*,*,SKYD__LO_HIT) remains zero. The GridSphere3D calls should ! fill these. This helps where a missing value would cause problems (e.g. in ! the star subtraction. ! The DO loop can be avoided, but that would increase the required size of the ! scratch arrays in GridSphere3D. ! Only do this for camera 2 (the other two don't get close enough to the pole ! to cause problems, I hope). ! The GridSphereWeight calls force averaging of all skybins with a count of ! 1 and 2 (< 3.0), in addition to the bins with a count of 1 (set to BadR4) by ! the division above. This reduces the noise for skybins where only a few ! pixel contribute. if (icam .eq. 2) then i = 2+10 ! Fill holes, use open grid rval = 360.0/SKYD__NXLO ! Gaussian halfwidth in degrees call GridSphere2D(1.0,SKYD__NXLO,SKYD__NYLO,1,lowres(1,1,SKYD__LO_SKY ),rval,i,0.0,0.0) call GridSphereRange(0.0,1.0) call GridSphereWeight(SKYD__NXLO,SKYD__NYLO,1,lowres(1,1,SKYD__LO_HIT ),3.0) call GridSphere2D(1.0,SKYD__NXLO,SKYD__NYLO,1,lowres(1,1,SKYD__LO_TIME),rval,i,0.0,0.0) call GridSphereRange(0.0,sngl(thisperiod)) call GridSphereWeight(SKYD__NXLO,SKYD__NYLO,1,lowres(1,1,SKYD__LO_HIT ),3.0) call GridSphere2D(1.0,SKYD__NXLO,SKYD__NYLO,1,lowres(1,1,SKYD__LO_SECS),rval,i,0.0,0.0) call GridSphereRange(-180.0,180.0) call GridSphereWeight(SKYD__NXLO,SKYD__NYLO,1,lowres(1,1,SKYD__LO_HIT ),3.0) call GridSphere2D(1.0,SKYD__NXLO,SKYD__NYLO,1,lowres(1,1,SKYD__LO_FOVX),rval,i,0.0,0.0) call GridSphereRange(-180.0,180.0) call GridSphereWeight(SKYD__NXLO,SKYD__NYLO,1,lowres(1,1,SKYD__LO_HIT ),3.0) call GridSphere2D(1.0,SKYD__NXLO,SKYD__NYLO,1,lowres(1,1,SKYD__LO_PSFN),rval,i,0.0,0.0) call GridSphereRange(-180.0,180.0) call GridSphereWeight(SKYD__NXLO,SKYD__NYLO,1,lowres(1,1,SKYD__LO_HIT ),3.0) call GridSphere2D(1.0,SKYD__NXLO,SKYD__NYLO,1,lowres(1,1,SKYD__LO_PSFE),rval,i,0.0,0.0) call GridSphereRange(-180.0,180.0) call GridSphereWeight(SKYD__NXLO,SKYD__NYLO,1,lowres(1,1,SKYD__LO_HIT ),3.0) call GridSphere2D(1.0,SKYD__NXLO,SKYD__NYLO,1,lowres(1,1,SKYD__LO_ROTX),rval,i,0.0,0.0) call GridSphereRange(-180.0,180.0) call GridSphereWeight(SKYD__NXLO,SKYD__NYLO,1,lowres(1,1,SKYD__LO_HIT ),3.0) call GridSphere2D(1.0,SKYD__NXLO,SKYD__NYLO,1,lowres(1,1,SKYD__LO_ROTY),rval,i,0.0,0.0) !call GridSphere2D(1.0,SKYD__NXLO,SKYD__NYLO,1,lowres(1,1,SKYD__LO_FRAC),rval,i,0.0,0.0) !call GridSphere2D(1.0,SKYD__NXLO,SKYD__NYLO,1,lowres(1,1,SKYD__LO_CUT ),rval,i,0.0,0.0) !call GridSphere2D(1.0,SKYD__NXLO,SKYD__NYLO,1,lowres(1,1,SKYD__LO_HIT ),rval,i,0.0,0.0) end if ! Ugly, but no way to avoid this, I think: set BadR4() values to badvalue (=zero). rbad = BadR4() call ArrR4Mask(SKYD__N_ALL,xmg,rbad,0.0,0.0,0.0,1.0,xmg) call ArrR4Mask(SKYD__N_LO*SKYD__LO_NMAP,lowres,rbad,0.0,0.0,0.0,1.0,lowres) ! The rotation angles could be slightly outside the range (-180,+180]. ! Make sure they are strictly inside this range. ! Note that this has to be done AFTER BadR4 has been replaced by badvalue (=zero) do j=1,SKYD__NYLO do i=1,SKYD__NXLO call smei_skyd_fixvalue(lowres(i,j,SKYD__LO_FOVX),-180.0,180.0)! Redundant ????? call smei_skyd_fixvalue(lowres(i,j,SKYD__LO_PSFN),-180.0,180.0) call smei_skyd_fixvalue(lowres(i,j,SKYD__LO_PSFE),-180.0,180.0) call smei_skyd_fixvalue(lowres(i,j,SKYD__LO_ROTX),-180.0,180.0) call smei_skyd_fixvalue(lowres(i,j,SKYD__LO_ROTY),-180.0,180.0) end do end do return C+ C NAME: C smei_skyd_fts C PURPOSE: C Writes orbital sky maps to FTS file C CATEGORY: C ucsd/camera/for/sky C CALLING SEQUENCE: entry smei_skyd_fts(sky_version,forbit_ext,forbit_ltd,hdrok,nframe,iframe,ibreak) C CALLS: C iGetLun, iFreeLun, ArrR4Mask, Say, say_fts C iSetFileSpec, iGetFileSpec, Time2Split, itrim C Time2System, Time2Str, smei_orbit_time2, smei_orbit_period2 C FTINIT, FTPHPR, FTPKYS, FTPKYJ, FTPKYE, FTPKYD, FTPKYG C FTPPRE, FTCRHD, FTCLOS, FTPPRB, FTPPRJ, iOSDeleteFile C PROCEDURE: C > Entry point in href=smei_skyd_sky=. C C > The equatorial maps use a linear scale to from array index i,j to ra,dec C i=1,nxeq, ra=[0,360), j=1,nyeq, dec=[-60,60) C i = cxeq+ra *d_eq ra = (i-cxeq)/d_eq C j = cyeq+dec*d_eq dec = (j-cyeq)/d_eq C > The polar maps use a polar transformation C i=1,nxeq, j=1,nyeq C i = cxpl+d_pl*(90-abs(dec))*cos(ra) C j = cypl+d_pl*(90-abs(dec))*sin(ra) C ra = atan( (j-cypl)/(i-cxpl) ) C 90-abs(dec) = sqrt( (i-cxpl)^2+(j-cypl)^2 ) C > The equatorial map is written as primary data, the north and C south polar maps as first and second extension, respectively. C > Added entries to FTS header: C C Ext 0 (equatorial map): C SMEI_SKY character smei_skyd software version number C DATE character File creation time C NAME character File name (no directory, no type) C CAMERA integer Camera C ORBIT double precision Orbit number C MAP character String describing map C CX double precision Array index for RA=0 C CY double precision Array index for DEC=0 C BD double precision Number of bins per degree C C Ext 1 (north pole map) C MAP character String describing map C CX double precision Array index for north pole C CY double precision Array index for north pole C BD double precision Number of bins per degree C C Ext 2 (south pole map) C MAP character String describing map C CX double precision Array index for south pole C CY double precision Array index for south pole C BD double precision Number of bins per degree C C Note that the values for both poles are the same. C C > After the high-resolution skymaps follow a group of (currently) six C low resolution maps: C C Ext 3-8 ('dirty' sky brightness, time, psf orientation w.r.t. equatorial C north and w.r.t. vernal equinox, fraction of discarded response, and C fraction of bad nodes) C MAP character String describing map C CX double precision Array index for RA=0 C CY double precision Array index for DEC=0 C BD double precision Number of bins per degree C MODIFICATION HISTORY: C NOV-2004, Paul Hick (UCSD/CASS) C APR-2005, Paul Hick (UCSD/CASS) C New set of camera quaternions for cam 1 and 2. C Ad hoc corrections to tx and ty in smei_ccd2cam. C Added 'onesky' keyword. C Extra keyword ONEORBIT to identify orbital maps. C NOV-2005, Paul Hick (UCSD/CASS) C Modified to accepts different glare maps for cameras 2 and 3. C DEC-2005, Paul Hick (UCSD/CASS) C Modified to deal with presence of -overwrite on cmd line. C Added Fits keyword MODE. C Added processing of SMEI_HTM keyword for existing sky maps C to force overwrite of maps with lower version number. C APR-2007, Paul Hick (UCSD/CASS) C Added keyword CLNEDGE, set to .FALSE.. Should be set .TRUE. C when straylight from sky just outside the fov is removed. C JUL-2007, Paul Hick (UCSD/CASS) C Renamed key SMEI_HTM to key SMEI_SKY (for version number) C DEC-2007, Paul Hick (UCSD/CASS) C Replaced keywords CX*,CY*,D_* keywords by WCS keywords C CRPIX*, CDELT*. Added other WCS keywords to make the C header WCS compliant: CRVAL*, CTYPE*, CUNIT*, RADESYS, C EQUINOX, MJD-OBS and (polar maps only) LONPOLE. C Added keyword BAD_DATA to headers for lowres maps. C Added keyword GAIN with camera gain based on time C set in MJD-OBS and TIME. C APR-2008, Paul Hick (UCSD/CASS; pphick@ucsd.edu) C Removed comma from comment string for keyword MAP C- if (bInRADecOff) then d_eq_ = SKYD__D_EQ n_eq_ = SKYD__N_EQ nxeq_ = SKYD__NXEQ nyeq_ = SKYD__NYEQ cxeq_ = origin(1,0) cyeq_ = origin(1,SKYD__NYEQ) do i=1,2 d_pl_(i) = SKYD__D_PL n_pl_(i) = SKYD__N_PL nxpl_(i) = SKYD__NXPL nypl_(i) = SKYD__NYPL cxpl_(i) = origin(1,SKYD__NXPL) cypl_(i) = origin(1,SKYD__NYPL) end do else call smei_skyd_scan(ra_min,ra_max,dec_min,dec_max, & d_eq_,n_eq_,nxeq_,nyeq_,cxeq_,cyeq_, & d_pl_,n_pl_,nxpl_,nypl_,cxpl_,cypl_, img,xmg) end if if (ibreak .lt. 0) then cStr = cSplit(:iSplit-1)//'sky'//cSplit(iSplit:jSplit)//cSkyType else i = Time2Str(SMEI__UT_FORMAT, tlastfrm , cTmp) cStr = cSplit(:iSplit-1)//'sky'//cSplit(iSplit:jSplit)//'_'//cTmp(:i)//cSkyType bExists = iSearch(11,cStr,cTmp) .ne. 0 end if ! Delete the file if it exists already if (bExists) then if (.not. bOverwrite .and. .not. bCheckVersion) call Say(cSayFts,'E','#'//cStr,'exists. Should not have happened') i = iOSDeleteFile(cStr) call Say(cSayFts,'I','#'//cStr,'deleted') end if ! Write FITS file iU = iGetLun(cStr) istat = 0 ! Initialize to zero ! Stays zero as long as there are no FITS errors ! Open a new Fits file. This will fail if the Fits file ! exists already. call Say(cSayFts,'I','#'//cStr,'is skymap file') call FTINIT(iU, cStr, 1, istat) if (istat .ne. 0) call say_fts(cSayFts,cSeverity,istat) ! Write equatorial map as primary data nAxes(1) = nxeq_ nAxes(2) = nyeq_ iBitPix = -32 ! Write primary header if (n_eq_ .ne. 0) then call FTPHPR(iU, bSimple, iBitPix, nAxis, nAxes, 0, 1, bExtend, istat) else call FTPHPR(iU, bSimple, iBitPix, 1, 0, 0, 1, bExtend, istat) endif if (istat .ne. 0) call say_fts(cSayFts,cSeverity,istat) ! Add to primary header ! Software version number indicator call FTPKYG(iU, 'SMEI_SKY', sky_version, 2, 'Software version number', istat) if (istat .ne. 0) call say_fts(cSayFts,cSeverity,istat) call FTPKYG(iU, 'SMEI_CAL', cal_version, 2, 'smei_cal version number', istat) if (istat .ne. 0) call say_fts(cSayFts,cSeverity,istat) call FTPKYG(iU, 'SMEI_ORB', orb_version, 2, 'smei_orb version number', istat) if (istat .ne. 0) call say_fts(cSayFts,cSeverity,istat) call Time2System(time) ! File creation time i = Time2Str('YYYY-MN-DD hh:mm',time,cTmp) call FTPKYS(iU, 'DATE', cTmp(:i), 'Creation UTC (CCCC-MM-DD) date of FITS header', istat) if (istat .ne. 0) call say_fts(cSayFts,cSeverity,istat) i = iSetFileSpec(cStr) i = iGetFileSpec(FIL__NAME,FIL__NAME,cTmp) ! File name call FTPKYS(iU, 'NAME', cTmp(:i), 'File name', istat) if (istat .ne. 0) call say_fts(cSayFts,cSeverity,istat) call FTPKYS(iU, 'IMG_TYPE', 'skymap_sky', 'Image type', istat) if (istat .ne. 0) call say_fts(cSayFts,cSeverity,istat) call FTPKYJ(iU, 'CAMERA', icam, 'Camera', istat)! Camera if (istat .ne. 0) call say_fts(cSayFts,cSeverity,istat) call FTPKYJ(iU, 'MODE', mode, 'Mode', istat) ! Mode if (istat .ne. 0) call say_fts(cSayFts,cSeverity,istat) i = Time2Str(SMEI__UT_FORMAT, tfirstfrm, cTmp) call FTPKYS(iU, 'STIME', cTmp(:i), 'First frame used', istat) if (istat .ne. 0) call say_fts(cSayFts,cSeverity,istat) i = Time2Str(SMEI__UT_FORMAT, tlastfrm , cTmp) call FTPKYS(iU, 'ETIME', cTmp(:i), 'Last frame used', istat) if (istat .ne. 0) call say_fts(cSayFts,cSeverity,istat) call FTPKYJ(iU, 'NFRAMES', nframe, 'Total nr frames', istat) if (istat .ne. 0) call say_fts(cSayFts,cSeverity,istat) call FTPKYJ(iU, 'IFRAMES', iframe, 'Nr frames used', istat) if (istat .ne. 0) call say_fts(cSayFts,cSeverity,istat) call FTPKYL(iU, 'CLNEDGE', .FALSE., 'Clean edge', istat) if (istat .ne. 0) call say_fts(cSayFts,cSeverity,istat) call FTPKYL(iU, 'ONEORBIT', .not. bOneSky, 'Orbital map', istat) if (istat .ne. 0) call say_fts(cSayFts,cSeverity,istat) if (bOneSky) then call Time2MJD(0,tfirstfrm,time) ! Obs start time in MJD call Time2Day8(0,time,dval) call FTPKYG(iU, 'MJD-OBS', dval, 7, 'First frame used (MJD)', istat) if (istat .ne. 0) call say_fts(cSayFts,cSeverity,istat) call Time2Delta(tfirstfrm,zero_gain_time,time) call Time2Day8(0,time,dval) dval = 1.0d0/inverse_zero_gain(icam)-0.01d0*dval/365.25d0 call FTPKYG(iU, 'GAIN', dval, 6, 'Camera gain at MJD-OBS', istat) else ! Orbit number call FTPKYJ(iU, 'ORBIT', thisorbit, 'Orbit number', istat) if (istat .ne. 0) call say_fts(cSayFts,cSeverity,istat) call FTPKYG(iU, 'LOWEXT', forbit_ext(1), 7, 'Start extended orbital fraction', istat) if (istat .ne. 0) call say_fts(cSayFts,cSeverity,istat) call FTPKYG(iU, 'HIGHEXT', forbit_ext(2), 7, 'End extended orbital fraction' , istat) if (istat .ne. 0) call say_fts(cSayFts,cSeverity,istat) call FTPKYG(iU, 'LOWLTD', forbit_ltd(1), 7, 'Start limited orbital fraction' , istat) if (istat .ne. 0) call say_fts(cSayFts,cSeverity,istat) call FTPKYG(iU, 'HIGHLTD', forbit_ltd(2), 7, 'End limited orbital fraction' , istat) if (istat .ne. 0) call say_fts(cSayFts,cSeverity,istat) call smei_orbit_time2(thisorbit,0.0d0,time) ! Exact orbit start time i = Time2Str(SMEI__UT_FORMAT//'.fff', time, cTmp) call FTPKYS(iU, 'TIME', cTmp(:i), 'Orbit start time', istat) if (istat .ne. 0) call say_fts(cSayFts,cSeverity,istat) ! Exact orbit period call FTPKYG(iU, 'PORBIT', thisperiod, 3, 'Orbit period in seconds', istat) if (istat .ne. 0) call say_fts(cSayFts,cSeverity,istat) call Time2MJD(0,time,time_) ! Orbit start time in MJD call Time2Day8(0,time_,dval) call FTPKYG(iU, 'MJD-OBS', dval, 7, 'Orbit start time (MJD)', istat) if (istat .ne. 0) call say_fts(cSayFts,cSeverity,istat) call Time2Delta(time,zero_gain_time,time_) call Time2Day8(0,time_,dval) dval = 1.0d0/inverse_zero_gain(icam)-0.01d0*dval/365.25d0 call FTPKYG(iU, 'GAIN', dval, 6, 'Camera gain at MJD-OBS', istat) end if ! Map description call FTPKYS(iU, 'MAP', 'Equatorial map RA=[0,360],DEC=[-60,+60]', 'Map description', istat) if (istat .ne. 0) call say_fts(cSayFts,cSeverity,istat) ! Begin WCS keywords call FTPKYS(iU, 'RADESYS', 'FK5', 'IAU 1984 system' , istat) if (istat .ne. 0) call say_fts(cSayFts,cSeverity,istat) call FTPKYG(iU, 'EQUINOX', 2000.0d0, 1, 'J2000 coordinates' , istat) if (istat .ne. 0) call say_fts(cSayFts,cSeverity,istat) call FTPKYS(iU, 'CTYPE1', 'RA', 'Axis type' , istat) if (istat .ne. 0) call say_fts(cSayFts,cSeverity,istat) call FTPKYS(iU, 'CTYPE2', 'DEC', 'Axis type' , istat) if (istat .ne. 0) call say_fts(cSayFts,cSeverity,istat) call FTPKYS(iU, 'CUNIT1', 'deg', 'Angles in degrees' , istat) if (istat .ne. 0) call say_fts(cSayFts,cSeverity,istat) call FTPKYS(iU, 'CUNIT2', 'deg', 'Angles in degrees' , istat) if (istat .ne. 0) call say_fts(cSayFts,cSeverity,istat) call FTPKYG(iU, 'CRPIX1', cxeq_, 2, 'Array index (base 1) for RA=0' , istat) if (istat .ne. 0) call say_fts(cSayFts,cSeverity,istat) call FTPKYG(iU, 'CRPIX2', cyeq_, 2, 'Array index (base 1) for DEC=0', istat) if (istat .ne. 0) call say_fts(cSayFts,cSeverity,istat) call FTPKYG(iU, 'CDELT1', 1.0d0/d_eq_, 2, 'Degrees per bin' , istat) if (istat .ne. 0) call say_fts(cSayFts,cSeverity,istat) call FTPKYG(iU, 'CDELT2', 1.0d0/d_eq_, 2, 'Degrees per bin' , istat) if (istat .ne. 0) call say_fts(cSayFts,cSeverity,istat) call FTPKYG(iU, 'CRVAL1', 0.0d0, 2, 'RA=0.0' , istat) if (istat .ne. 0) call say_fts(cSayFts,cSeverity,istat) call FTPKYG(iU, 'CRVAL2', 0.0d0, 2, 'DEC=0.0' , istat) if (istat .ne. 0) call say_fts(cSayFts,cSeverity,istat) ! End WCS keywords call FTPKYG(iU, 'BAD_DATA', badvalue, 2, 'Missing data flag', istat) if (istat .ne. 0) call say_fts(cSayFts,cSeverity,istat) ! Additions to primary header set in smei_frm_ok ! Modified nped_min and/or ndark_min parameters if ( hdrok(SMEI__HDROK_NPED_MIN) .ne. 0.0d0 .or. hdrok(SMEI__HDROK_NDARK_MIN) .ne. 0.0d0 ) then call FTPKYJ(iU, 'NPED_MIN',nint(hdrok(SMEI__HDROK_NPED_MIN )), 'Non-default nped_min threshold' , istat) if (istat .ne. 0) call say_fts(cSayFts,cSeverity,istat) call FTPKYJ(iU, 'NDRK_MIN',nint(hdrok(SMEI__HDROK_NDARK_MIN)), 'Non-default ndark_min threshold', istat) if (istat .ne. 0) call say_fts(cSayFts,cSeverity,istat) end if ! Skipped open-shutter frames every time shutter opens if ( hdrok(SMEI__HDROK_SHUTTER_SKIP) .ne. 0.0d0 ) then call FTPKYJ(iU, 'SHUT_SKIP' ,nint(hdrok(SMEI__HDROK_SHUTTER_SKIP )), 'Nr open-shutter frames skipped' , istat) if (istat .ne. 0) call say_fts(cSayFts,cSeverity,istat) end if ! Exclusion areas for Sun, Moon and Venus call ArrR8Copy(SMEI__HDROK_NBOX*SMEI__HDROK_NBODY,hdrok(SMEI__HDROK_BODYEXCL),bodyexcl) do i=1,SMEI__HDROK_NBODY if (bodyexcl(1,i) .ne. 0.0d0 .or. bodyexcl(2,i) .ne. 0.0d0 .or. & bodyexcl(3,i) .ne. 0.0d0 .or. bodyexcl(4,i) .ne. 0.0d0 ) then call FTPKYG(iU,'TX1'//SMEI__HDROK_BODIES(i)(:5),bodyexcl(1,i),3,'Min theta_x excl '//SMEI__HDROK_BODIES(i), istat) call FTPKYG(iU,'TX2'//SMEI__HDROK_BODIES(i)(:5),bodyexcl(2,i),3,'Max theta_x excl '//SMEI__HDROK_BODIES(i), istat) call FTPKYG(iU,'TY1'//SMEI__HDROK_BODIES(i)(:5),bodyexcl(3,i),3,'Min theta_y excl '//SMEI__HDROK_BODIES(i), istat) call FTPKYG(iU,'TY2'//SMEI__HDROK_BODIES(i)(:5),bodyexcl(4,i),3,'Max theta_y excl '//SMEI__HDROK_BODIES(i), istat) end if end do ! Data array call FTPPRE(iU, 0, 1, n_eq_, xmg, istat) ! Array if (istat .ne. 0) call say_fts(cSayFts,cSeverity,istat) ! Write polar maps as first and second extension iBitPix = -32 do i=1,2 nAxes(1) = nxpl_(i) nAxes(2) = nypl_(i) call FTCRHD(iU, istat) if (istat .ne. 0) call say_fts(cSayFts,cSeverity,istat) if (n_pl_(i) .ne. 0) then call FTPHPR(iU, bSimple, iBitPix, nAxis, nAxes, 0, 1, bExtend, istat) else call FTPHPR(iU, bSimple, iBitPix, 1, 0, 0, 1, bExtend, istat) end if if (istat .ne. 0) call say_fts(cSayFts,cSeverity,istat) ! Scaling constants for north/south pole map call FTPKYS(iU, 'MAP', 'Map of '//cPole(i)//' pole DEC=['//cSign(i)//'50,'//cSign(i)//'90]', & 'Map description', istat) if (istat .ne. 0) call say_fts(cSayFts,cSeverity,istat) ! Begin WCS keywords call FTPKYS(iU, 'RADESYS', 'FK5', 'IAU 1984 system' , istat) if (istat .ne. 0) call say_fts(cSayFts,cSeverity,istat) call FTPKYG(iU, 'EQUINOX', 2000.0d0, 1, 'J2000 coordinates' , istat) if (istat .ne. 0) call say_fts(cSayFts,cSeverity,istat) call FTPKYS(iU, 'CTYPE1', 'RA---ARC', 'Axis type: equidistant azimuthal', istat) if (istat .ne. 0) call say_fts(cSayFts,cSeverity,istat) call FTPKYS(iU, 'CTYPE2', 'DEC--ARC', 'Axis type: equidistant azimuthal', istat) if (istat .ne. 0) call say_fts(cSayFts,cSeverity,istat) call FTPKYS(iU, 'CUNIT1', 'deg', 'Angles in degrees' , istat) if (istat .ne. 0) call say_fts(cSayFts,cSeverity,istat) call FTPKYS(iU, 'CUNIT2', 'deg', 'Angles in degrees' , istat) if (istat .ne. 0) call say_fts(cSayFts,cSeverity,istat) call FTPKYG(iU, 'CRPIX1', cxpl_(i), 2, 'Array index (base 1) for '//cPole(i)//' pole', istat) if (istat .ne. 0) call say_fts(cSayFts,cSeverity,istat) call FTPKYG(iU, 'CRPIX2', cypl_(i), 2, 'Array index (base 1) for '//cPole(i)//' pole', istat) if (istat .ne. 0) call say_fts(cSayFts,cSeverity,istat) j = 3-2*i ! +1 = NP; -1 is SP call FTPKYG(iU, 'CDELT1', dsign(i)/d_pl_(i), 2, 'Degrees per bin' , istat) if (istat .ne. 0) call say_fts(cSayFts,cSeverity,istat) call FTPKYG(iU, 'CDELT2', 1.0d0/d_pl_(i), 2, 'Degrees per bin' , istat) if (istat .ne. 0) call say_fts(cSayFts,cSeverity,istat) call FTPKYG(iU, 'CRVAL1', dsign(i)*90.0d0, 2, 'RA of '//cPole(i)//' pole' , istat) if (istat .ne. 0) call say_fts(cSayFts,cSeverity,istat) call FTPKYG(iU, 'CRVAL2', dsign(i)*90.0d0, 2, 'DEC of '//cPole(i)//' pole', istat) if (istat .ne. 0) call say_fts(cSayFts,cSeverity,istat) call FTPKYG(iU, 'LONPOLE', 0.0d0, 2, 'Native longitude of celestial pole', istat) if (istat .ne. 0) call say_fts(cSayFts,cSeverity,istat) ! End WCS keywords call FTPKYG(iU, 'BAD_DATA', badvalue, 2, 'Missing data flag', istat) if (istat .ne. 0) call say_fts(cSayFts,cSeverity,istat) call FTPPRE(iU, 0, 1, n_pl_(i), xmg(n_eq_+(i-1)*n_pl_(1)+1), istat) ! Write array if (istat .ne. 0) call say_fts(cSayFts,cSeverity,istat) end do ! img(i) is at most 2^level do i=1,n_eq_+n_pl_(1)+n_pl_(2) jmg(i) = min( nint(img(i)), 127 ) end do ! Write equatorial counts as 3rd extension nAxes(1) = nxeq_ nAxes(2) = nyeq_ iBitPix = 8 call FTCRHD(iU, istat) if (istat .ne. 0) call say_fts(cSayFts,cSeverity,istat) call FTPHPR(iU, bSimple, iBitPix, nAxis, nAxes, 0, 1, bExtend, istat) if (istat .ne. 0) call say_fts(cSayFts,cSeverity,istat) call FTPKYS(iU, 'MAP', 'Counts for equatorial map RA=[0,360],DEC=[-60,+60]', 'Map description', istat) ! Begin WCS keywords call FTPKYS(iU, 'RADESYS', 'FK5', 'IAU 1984 system' , istat) if (istat .ne. 0) call say_fts(cSayFts,cSeverity,istat) call FTPKYG(iU, 'EQUINOX', 2000.0d0, 1, 'J2000 coordinates' , istat) if (istat .ne. 0) call say_fts(cSayFts,cSeverity,istat) call FTPKYS(iU, 'CTYPE1', 'RA', 'Axis type' , istat) if (istat .ne. 0) call say_fts(cSayFts,cSeverity,istat) call FTPKYS(iU, 'CTYPE2', 'DEC', 'Axis type' , istat) if (istat .ne. 0) call say_fts(cSayFts,cSeverity,istat) call FTPKYS(iU, 'CUNIT1', 'deg', 'Angles in degrees' , istat) if (istat .ne. 0) call say_fts(cSayFts,cSeverity,istat) call FTPKYS(iU, 'CUNIT2', 'deg', 'Angles in degrees' , istat) if (istat .ne. 0) call say_fts(cSayFts,cSeverity,istat) call FTPKYG(iU, 'CRPIX1', cxeq_, 2, 'Array index (base 1) for RA=0' , istat) if (istat .ne. 0) call say_fts(cSayFts,cSeverity,istat) call FTPKYG(iU, 'CRPIX2', cyeq_, 2, 'Array index (base 1) for DEC=0', istat) if (istat .ne. 0) call say_fts(cSayFts,cSeverity,istat) call FTPKYG(iU, 'CDELT1', 1.0d0/d_eq_, 2, 'Degrees per bin' , istat) if (istat .ne. 0) call say_fts(cSayFts,cSeverity,istat) call FTPKYG(iU, 'CDELT2', 1.0d0/d_eq_, 2, 'Degrees per bin' , istat) if (istat .ne. 0) call say_fts(cSayFts,cSeverity,istat) call FTPKYG(iU, 'CRVAL1', 0.0d0, 2, 'RA=0.0' , istat) if (istat .ne. 0) call say_fts(cSayFts,cSeverity,istat) call FTPKYG(iU, 'CRVAL2', 0.0d0, 2, 'DEC=0.0' , istat) if (istat .ne. 0) call say_fts(cSayFts,cSeverity,istat) ! End WCS keywords call FTPPRB(iU, 0, 1, n_eq_, jmg, istat) ! Array if (istat .ne. 0) call say_fts(cSayFts,cSeverity,istat) ! Write polar counts as 4th and 5th extension iBitPix = 8 do i=1,2 nAxes(1) = nxpl_(i) nAxes(2) = nypl_(i) call FTCRHD(iU, istat) if (istat .ne. 0) call say_fts(cSayFts,cSeverity,istat) call FTPHPR(iU, bSimple, iBitPix, nAxis, nAxes, 0, 1, bExtend, istat) if (istat .ne. 0) call say_fts(cSayFts,cSeverity,istat) call FTPKYS(iU, 'MAP', 'Counts for map of '//cPole(i)//' pole DEC=['//cSign(i)//'50,'//cSign(i)//'90]', & 'Map description', istat) if (istat .ne. 0) call say_fts(cSayFts,cSeverity,istat) ! Begin WCS keywords call FTPKYS(iU, 'RADESYS', 'FK5', 'IAU 1984 system' , istat) if (istat .ne. 0) call say_fts(cSayFts,cSeverity,istat) call FTPKYG(iU, 'EQUINOX', 2000.0d0, 1, 'J2000 coordinates' , istat) if (istat .ne. 0) call say_fts(cSayFts,cSeverity,istat) call FTPKYS(iU, 'CTYPE1', 'RA---ARC', 'Axis type: equidistant azimuthal', istat) if (istat .ne. 0) call say_fts(cSayFts,cSeverity,istat) call FTPKYS(iU, 'CTYPE2', 'DEC--ARC', 'Axis type: equidistant azimuthal', istat) if (istat .ne. 0) call say_fts(cSayFts,cSeverity,istat) call FTPKYS(iU, 'CUNIT1', 'deg', 'Angles in degrees' , istat) if (istat .ne. 0) call say_fts(cSayFts,cSeverity,istat) call FTPKYS(iU, 'CUNIT2', 'deg', 'Angles in degrees' , istat) if (istat .ne. 0) call say_fts(cSayFts,cSeverity,istat) call FTPKYG(iU, 'CRPIX1', cxpl_(i), 2, 'Array index (base 1) for '//cPole(i)//' pole', istat) if (istat .ne. 0) call say_fts(cSayFts,cSeverity,istat) call FTPKYG(iU, 'CRPIX2', cypl_(i), 2, 'Array index (base 1) for '//cPole(i)//' pole', istat) if (istat .ne. 0) call say_fts(cSayFts,cSeverity,istat) j = 3-2*i ! +1 = NP; -1 is SP call FTPKYG(iU, 'CDELT1', dsign(i)/d_pl_(i), 2, 'Degrees per bin' , istat) if (istat .ne. 0) call say_fts(cSayFts,cSeverity,istat) call FTPKYG(iU, 'CDELT2', 1.0d0/d_pl_(i), 2, 'Degrees per bin' , istat) if (istat .ne. 0) call say_fts(cSayFts,cSeverity,istat) call FTPKYG(iU, 'CRVAL1', dsign(i)*90.0d0, 2, 'RA of '//cPole(i)//' pole' , istat) if (istat .ne. 0) call say_fts(cSayFts,cSeverity,istat) call FTPKYG(iU, 'CRVAL2', dsign(i)*90.0d0, 2, 'DEC of '//cPole(i)//' pole', istat) if (istat .ne. 0) call say_fts(cSayFts,cSeverity,istat) call FTPKYG(iU, 'LONPOLE', 0.0d0, 2, 'Native longitude of celestial pole', istat) if (istat .ne. 0) call say_fts(cSayFts,cSeverity,istat) ! End WCS keywords call FTPPRB(iU, 0, 1, n_pl_(i), jmg(n_eq_+(i-1)*n_pl_(1)+1), istat) ! Write array if (istat .ne. 0) call say_fts(cSayFts,cSeverity,istat) end do nAxes(1) = SKYD__NXLO ! Write low resolution maps nAxes(2) = SKYD__NYLO iBitPix = -32 do i=1,SKYD__LO_NMAP j = itrim(cFitsID(i)) if (j .gt. 0) then call FTCRHD(iU, istat) if (istat .ne. 0) call say_fts(cSayFts,cSeverity,istat) call FTPHPR(iU, bSimple, iBitPix, nAxis, nAxes, 0, 1, bExtend, istat) if (istat .ne. 0) call say_fts(cSayFts,cSeverity,istat) call FTPKYS(iU, 'MAP', cFitsID(i)(:j)//' RA=[0,360],DEC=[-90,+90]', 'Map description', istat) if (istat .ne. 0) call say_fts(cSayFts,cSeverity,istat) if (cFitsID(i)(:j) .eq. 'Time (sec since TORIGIN)') then ! Add time origin for array k = Time2Str(SMEI__UT_FORMAT, torigin, cTmp) call FTPKYS(iU, 'TORIGIN', cTmp(:k), 'Time origin', istat) if (istat .ne. 0) call say_fts(cSayFts,cSeverity,istat) ! Exact orbit period call FTPKYG(iU, 'PORBIT', thisperiod, 3, 'Orbit period in seconds', istat) if (istat .ne. 0) call say_fts(cSayFts,cSeverity,istat) end if ! Begin WCS keywords call FTPKYS(iU, 'RADESYS', 'FK5', 'IAU 1984 system' , istat) if (istat .ne. 0) call say_fts(cSayFts,cSeverity,istat) call FTPKYG(iU, 'EQUINOX', 2000.0d0, 1, 'J2000 coordinates' , istat) if (istat .ne. 0) call say_fts(cSayFts,cSeverity,istat) call FTPKYS(iU, 'CTYPE1', 'RA', 'Axis type' , istat) if (istat .ne. 0) call say_fts(cSayFts,cSeverity,istat) call FTPKYS(iU, 'CTYPE2', 'DEC', 'Axis type' , istat) if (istat .ne. 0) call say_fts(cSayFts,cSeverity,istat) call FTPKYS(iU, 'CUNIT1', 'deg', 'Angles in degrees' , istat) if (istat .ne. 0) call say_fts(cSayFts,cSeverity,istat) call FTPKYS(iU, 'CUNIT2', 'deg', 'Angles in degrees' , istat) if (istat .ne. 0) call say_fts(cSayFts,cSeverity,istat) call FTPKYG(iU, 'CRPIX1',origin(1, 0), 2, 'Array index (base 1) for RA=0' , istat) if (istat .ne. 0) call say_fts(cSayFts,cSeverity,istat) call FTPKYG(iU, 'CRPIX2',origin(1,SKYD__NYLO), 2, 'Array index (base 1) for DEC=0', istat) if (istat .ne. 0) call say_fts(cSayFts,cSeverity,istat) call FTPKYG(iU, 'CDELT1', 1.0d0/SKYD__D_LO, 2, 'Degrees per bin' , istat) if (istat .ne. 0) call say_fts(cSayFts,cSeverity,istat) call FTPKYG(iU, 'CDELT2', 1.0d0/SKYD__D_LO, 2, 'Degrees per bin' , istat) if (istat .ne. 0) call say_fts(cSayFts,cSeverity,istat) call FTPKYG(iU, 'CRVAL1', 0.0d0, 2, 'RA=0.0' , istat) if (istat .ne. 0) call say_fts(cSayFts,cSeverity,istat) call FTPKYG(iU, 'CRVAL2', 0.0d0, 2, 'DEC=0.0' , istat) if (istat .ne. 0) call say_fts(cSayFts,cSeverity,istat) ! End WCS keywords call FTPKYG(iU, 'BAD_DATA', badvalue, 2, 'Missing data flag', istat) if (istat .ne. 0) call say_fts(cSayFts,cSeverity,istat) call FTPPRE(iU, 0, 1, SKYD__N_LO, lowres(1,1,i), istat) ! Write array if (istat .ne. 0) call say_fts(cSayFts,cSeverity,istat) end if end do nAxes(1) = nX ! Write bad pixel map nAxes(2) = nY iBitPix = 32 call FTCRHD(iU, istat) if (istat .ne. 0) call say_fts(cSayFts,cSeverity,istat) call FTPHPR(iU, bSimple, iBitPix, nAxis, nAxes, 0, 1, bExtend, istat) if (istat .ne. 0) call say_fts(cSayFts,cSeverity,istat) call FTPKYS(iU, 'MAP', 'Bad pixel map', 'Map description', istat) if (istat .ne. 0) call say_fts(cSayFts,cSeverity,istat) call FTPPRJ(iU, 0, 1, nX*nY, badpix, istat) ! Write array if (istat .ne. 0) call say_fts(cSayFts,cSeverity,istat) call FTCLOS(iU, istat) ! Close Fits file if (istat .ne. 0) call say_fts(cSayFts,cSeverity,istat) iU = iFreeLun(iU) return end