C+ C NAME: C smei_skyd_scan C PURPOSE: C Extracts part of skymap inside box bounded in right ascension C and declination C CATEGORY: C camera/for/sky C CALLING SEQUENCE: subroutine 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) C INPUTS: C ra_min double precision minimum right ascension (deg) in [0,360] C ra_max double precision maximum right ascension (deg) in [0,360] C If ra_min > ra_max then RA=[ra_min,360] C and RA[0,ra_max] are extracted C dec_min double precision minimum declination C dec_max double precision maximum declination C xmg(*) real sky map C img(*) real sky map counts C OUTPUTS: C For the eqautorial map: C C d_eq_ double precision bin size (skybins per degree) C n_eq_ integer nr of sky bins (=nxeq_*nyeq_) C nxeq_ integer nr of sky bins in RA direction C nyeq_ integer nr of sky bins in dec direction C cxeq_ double precision array index for RA=0 C cyeq_ double precision array index for dec=0 C C For the polar maps: C C d_pl_ double precision radial bin size (skybins per degree) C n_pl_(2) integer nr of sky bins in polar maps (=nxpl_*nypl_) C nxpl_(2) integer nr of sky bins in x-direction of polar map C nypl_(2) integer nr of sky bins in y-direction of polar map C cxpl_(2) double precision horizontal array index of pole C cypl_(2) double precision vertical array index of pole C C xmg(*) real updated sky maps C img(*) real updated sky map counts C INCLUDE: include 'smei_skyd_dim.h' C PROCEDURE: C MODIFICATION HISTORY: C MAR-2006, Paul Hick (UCSD/CASS; pphick@ucsd.edu) C- double precision ra_min double precision ra_max double precision dec_min double precision dec_max double precision d_eq_ integer n_eq_ integer nxeq_ integer nyeq_ double precision cxeq_ double precision cyeq_ real img(*) real xmg(*) real xrow(SKYD__NXEQ) integer irow(SKYD__NXEQ) equivalence (xrow,irow) double precision cxpl double precision cypl double precision d_pl_(2) integer n_pl_(2) integer nxpl_(2) integer nypl_(2) double precision cxpl_(2) double precision cypl_(2) double precision px double precision py double precision dmin double precision dmax double precision PLMaxRad include 'smei_skyd_fnc.h' ip_eq_(i,j ) = (j-1)*nxeq_ +i ip_pl_(i,j,k) = n_eq_+(k-1)*n_pl_(1)+(j-1)*nxpl_(k)+i d_eq_ = SKYD__D_EQ cxeq_ = origin(1,0) cyeq_ = origin(1,SKYD__NYEQ) ilow = nint(origin(1,0)+ra_min *SKYD__D_EQ+0.5d0) ihigh = nint(origin(1,0)+ra_max *SKYD__D_EQ-0.5d0) jlow = nint(origin(1,SKYD__NYEQ)+dec_min*SKYD__D_EQ+0.5d0) jhigh = nint(origin(1,SKYD__NYEQ)+dec_max*SKYD__D_EQ-0.5d0) if (ilow .gt. ihigh) then ! The requested RA range bridges the 0/360 deg meridian ! Shuffle the array to put ilow in the first column n = SKYD__NXEQ-ilow+1 ! # grid points between ra_min and 360 deg do j=1,SKYD__NYEQ call ArrR4Copy(ilow-1,xmg(ip_eq(1,1,j)),xrow) call ArrR4Copy(n,xmg(ip_eq(1,ilow,j)),xmg(ip_eq(1,1,j))) call ArrR4Copy(ilow-1,xrow,xmg(ip_eq(1,n+1,j))) call ArrR4Copy(ilow-1,img(ip_eq(1,1,j)),irow) call ArrR4Copy(n,img(ip_eq(1,ilow,j)),img(ip_eq(1,1,j))) call ArrR4Copy(ilow-1,irow,img(ip_eq(1,n+1,j))) end do cxeq_ = cxeq_-ilow+1 ilow = 1 ihigh = ihigh+n end if nxeq_ = ihigh-ilow+1 nyeq_ = jhigh-jlow+1 n_eq_ = nxeq_*nyeq_ if (n_eq_ .eq. 0) then nxeq_ = 0 nyeq_ = 0 cxeq_ = 0.0d0 cyeq_ = 0.0d0 d_eq_ = 0.0d0 else if (n_eq_ .ne. SKYD__N_EQ) then do j=1,nyeq_ do i=1,nxeq_ xmg(ip_eq_(i,j)) = xmg(ip_eq(1,ilow+i-1,jlow+j-1)) img(ip_eq_(i,j)) = img(ip_eq(1,ilow+i-1,jlow+j-1)) end do end do cxeq_ = cxeq_-ilow+1 cyeq_ = cyeq_-jlow+1 end if cxpl = origin(1,SKYD__NXPL) cypl = origin(1,SKYD__NYPL) px = cxpl-0.5d0 py = cypl-0.5d0 PLMaxRad = dsqrt(px*px+py*py)/SKYD__D_PL do k=1,2 ipm = 3-2*k ! +1 at North Pole, -1 at South Pole dmin = SKYD__D_PL*(90.0d0-ipm*dec_min) dmax = SKYD__D_PL*(90.0d0-ipm*dec_max) i1 = nint(cxpl+dmin*dcosd(ra_min)) i2 = nint(cxpl+dmin*dcosd(ra_max)) i3 = nint(cxpl+dmax*dcosd(ra_min)) i4 = nint(cxpl+dmax*dcosd(ra_max)) j1 = nint(cypl+dmin*dsind(ra_min)) j2 = nint(cypl+dmin*dsind(ra_max)) j3 = nint(cypl+dmax*dsind(ra_min)) j4 = nint(cypl+dmax*dsind(ra_max)) ilow = min(i1,i2,i3,i4) ihigh = max(i1,i2,i3,i4) jlow = min(j1,j2,j3,j4) jhigh = max(j1,j2,j3,j4) if ( (ra_min .le. ra_max .and. (ra_min .le. 0.0d0 .and. 0.0d0 .le. ra_max)) .or. & (ra_min .gt. ra_max .and. (ra_min .le. 0.0d0 .or. 0.0d0 .le. ra_max)) ) & ihigh = (2-k)*nint(cxpl+dmin)+(k-1)*nint(cxpl+dmax) if ( (ra_min .le. ra_max .and. (ra_min .le. 90.0d0 .and. 90.0d0 .le. ra_max)) .or. & (ra_min .gt. ra_max .and. (ra_min .le. 90.0d0 .or. 90.0d0 .le. ra_max)) ) & jhigh = (2-k)*nint(cypl+dmin)+(k-1)*nint(cypl+dmax) if ( (ra_min .le. ra_max .and. (ra_min .le. 180.0d0 .and. 180.0d0 .le. ra_max)) .or. & (ra_min .gt. ra_max .and. (ra_min .le. 180.0d0 .or. 180.0d0 .le. ra_max)) ) & ilow = (2-k)*nint(cxpl-dmin)+(k-1)*nint(cxpl-dmax) if ( (ra_min .le. ra_max .and. (ra_min .le. 270.0d0 .and. 270.0d0 .le. ra_max)) .or. & (ra_min .gt. ra_max .and. (ra_min .le. 270.0d0 .or. 270.0d0 .le. ra_max)) ) & jlow = (2-k)*nint(cypl-dmin)+(k-1)*nint(cypl-dmax) ilow = max(ilow ,1) ihigh = max( min(ihigh,SKYD__NXPL), ilow-1 ) jlow = max(jlow ,1) jhigh = max( min(jhigh,SKYD__NYPL), jlow-1 ) d_pl_(k) = SKYD__D_PL cxpl_(k) = cxpl cypl_(k) = cypl nxpl_(k) = ihigh-ilow+1 nypl_(k) = jhigh-jlow+1 n_pl_(k) = nxpl_(k)*nypl_(k) if (n_pl_(k) .eq. 0) then nxpl_(k) = 0 nypl_(k) = 0 cxpl_(k) = 0.0d0 cypl_(k) = 0.0d0 d_pl_(k) = 0.0d0 else if (n_pl_(k) .ne. SKYD__N_PL) then do j=1,nypl_(k) do i=1,nxpl_(k) xmg(ip_pl_(i,j,k)) = xmg(ip_pl(1,ilow+i-1,jlow+j-1,k)) img(ip_pl_(i,j,k)) = img(ip_pl(1,ilow+i-1,jlow+j-1,k)) end do end do cxpl_(k) = cxpl_(k)-ilow+1 cypl_(k) = cypl_(k)-jlow+1 end if end do return end