C+ C NAME: C rot_image_2ra_dec.f C PURPOSE: C Read the sidereal time and the mean rotation determined from Polaris and Matthew's center to place the Pole at the map top C CATEGORY: C I C CALLING SEQUENCE: C call rot_image_2pole(indx,NXY,inum1,inum2,inum3,inum4,NOBS,NRA,NDEC,centx,centy,dmax,rotrmm,PDIST,Doys8,ST, C & aimagebak,aimage,bakmin,aimagexyr,aimagexyn,aimagexyo,cimradec,ncimradec,ocimradec,XDIST,YDIST) C C INPUTS: C indx integer index of what to do C indx = 0 Rotate the image so that the pole is at the top (always do this) and make a composite C indx = 1 indx 0 + shift index in y so pole is at the a given pixel location and make a composite C indx = 2 indx 0 + 1 + warp the image relative to the y center pixel location and make a composite C indx = 3 indx 0 + 1 + 2 + subtract an image and make a composite C INDX = 4 indx 4 read a rotated image, skip the rotate and make an RA and Dec composite C indx = 5 indx 5 read a rotated image, skip the rotate and make an RA and Dec composite and subtract an image C C NXY integer Maximum number of x or y values in each grd image C inum integer Image number C NOBS integer Number of images C centx(NOBS) real Image X centers C centy(NOBS) real Image Y centers C dmax real Distance of image from center to cut off analysis C PDIST(NOBS) real Unrotated distance of the pole from the image center C ST(NOBS) real Sidereal time C rotrmm(NOBS) real Mean rotation from Matthew's center to get the pole to the top of the image C aimagexy(NXY,NXY) real Original image brightness values in ADU C OUTPUTS: C aimagexyr(NXY,NXY)real Combined, rotated, and refined image brightness values in ADU C aimagexyn(NXY,NXY)real Number of image brightness hits per pixel in the combined, rotated, and refined images C aimagexyo(NXY,NXY)real Rejected average image brightness values of the combined, rotated, and refined images C C FUNCTIONS/SUBROUTINES: C C MODIFICATION HISTORY: C May-2023, Bernard Jackson (UCSD) C- subroutine rot_image_2ra_dec(indx,NXY,inum1,inum2,inum3,inum4,NOBS,NRA,NDEC,centx,centy,dmax,rotrmm,PDIST,Doys8,ST, & aimagebak,aimage,bakmin,aimagexyr,aimagexyn,aimagexyo,cimradec,ncimradec,ocimradec,XDIST,YDIST,BACKM,PED) real centx (NOBS), ! Image X centers & centy (NOBS), ! Image Y centers & rotrmm(NOBS), ! Mean rotation from Matthew's center to get the pole to the top of the image & PDIST (NOBS), ! Mean distance to the pole from the centx and centy location in pixels in original images & ST (NOBS), ! Sidereal time in degrees & XDIST (NOBS), ! distxc & YDIST (NOBS), ! final value of distyc & BACKM (NOBS), ! Background minimum near the center of the background & PED (NOBS) ! A pedestal determined by an image portion distant from the image center integer IVAL (NOBS) ! The num3 image at the center of the composite real aimagexyr(NXY,NXY), ! Rotated combined image brightness values in ADU & aimagexyn(NXY,NXY), ! Number of image brightness hits per pixel in the rotated image & aimagexyo(NXY,NXY), ! Brightness of pixels removed & aimagebak(NXY,NXY) ! Background image to subtract real aimagexyr1(NXY,NXY), ! Dummy rotated image brightness values in ADU & aimagexys (NXY,NXY) ! Dummy combined rotated image brightness values in ADU real aimage (NXY,NXY,NOBS), ! Original image set to rotate and combine & aimagexyrk(NXY,NXY,NOBS) ! Rotated image for each image real sigv(inum4-inum1+1) real DEC (NXY,NXY), & HA (NXY,NXY), & RA (NXY,NXY), & DEC1 (NXY,NXY), & HA1 (NXY,NXY), & RA1 (NXY,NXY) real cimradec (NRA,NDEC), ! composite images in RA, DEC & ocimradec(NRA,NDEC) ! thrown out composite images in RA, DEC integer ncimradec(NRA,NDEC) ! Number of hits in composite images in RA, DEC C real DECC (NXY,NXY,NOBS), C & HAA (NXY,NXY,NOBS), C & RAA (NXY,NXY,NOBS), C & DECCC(NXY,NXY,NOBS), C & HAAA (NXY,NXY,NOBS), C & RAAA (NXY,NXY,NOBS) real*8 Doys8(NOBS),Doy8,dLngSun,dLatSun,dDisSun ! Used in SunNewcom8 integer imagexyrk(NXY,NXY,NOBS) ! Total rotated number for each image integer imagexyr1(NXY,NXY) ! Dummy total rotated number for the combined image C ill = 5 ill = 6000 sig3 = 0.002700 ! Three standard deviations sig4 = 0.00006334 ! Four standard deviations allim = 100.0 ! minimum brightness below which image data are not accepted as brightness in combined images aulim = 65534.0 ! maximum brightness above which image data are not accepted as brightness in combined images PDIST1 = 576.0 inumi = 20 inumid2=inumi/2 inumhi = 9 RAD = 0.0 DECD = 0.0 print*,' ' write(*,'(3(A,I5))') 'Into rot_image_2ra_dec.f, inum1 =',inum1,', inum2 =',inum2,', ill =',ill write(*,'(5(A,I5))') 'Into rot_image_2ra_dec.f, inum3 =',inum3,', inum4 =',inum4,', inumh1 =',inumhi,' inumi =',inumi print*, ' ' do J=1,NXY do I=1,NXY aimagexyrk(I,J,1) = 0.0 imagexyrk (I,J,1) = 0 end do end do DCON = 2048.0/360.0 C iecliptic = 0 iecliptic = 1 0 iheliographic = 0 C iheliographic = 1 ihammer = 1 if(ihammer.eq.1) then if(iecliptic.eq.1) then ! Provide a mark at the Sun's location at image center KKK = 1 do IS=1,9 ! This provides a spot at the location of the Sun. do JS=1,15 rLatSun = 0.0 rLatP = ((rLatSun+(-7.0+float(JS))/7.0)+90.0)*DCON rLngP = ((0.0 +(-7.0+float(IS))/3.5)+180)*DCON aimagexyrk(nint(rLngP),nint(rLatP),KKK) = 4000.0 ! make this point bright in the KKK image imagexyrk (nint(rLngP),nint(rLatP),KKK) = 1 ! make sure this point is counted in the KKK image C print*,aimagexyrk(nint(rLngP),nint(rLatP),KKK),KKK,nint(rLngP),nint(rLatP) end do end do end if end if YDIS = 0.0 ANYDIS = 0.0 inum = -1 KKK = 0 do KK=inum1+inumhi,inum4,inumi K = KK + inum + 1 KKK = KKK + 1 distxci = centx(K) C distxc = centx(K) distyci = centy(K) Ifirst = 0 Jfirst = 0 do J=1,NXY do I=1,NXY if(KKK.gt.1) aimagexyrk(I,J,KKK) = 0.0 ! The first of these are zeroed at the first image set if(KKK.gt.1) imagexyrk (I,J,KKK) = 0 ! The first of these are zeroed at the first image set C DECCC (I,J,K) = 0.0 C HAAA (I,J,K) = 0.0 C RAAA (I,J,K) = 0.0 C DECC (I,J,K) = 0.0 C HAA (I,J,K) = 0.0 C RAA (I,J,K) = 0.0 if(KKK.eq.1) then aimagexyr (I,J) = 0.0 aimagexyr1(I,J) = 0.0 aimagexys (I,J) = 0.0 aimagexyn (I,J) = 0.0 aimagexyo (I,J) = 0.0 imagexyr1 (I,J) = 0 DEC1 (I,J) = 0.0 HA1 (I,J) = 0.0 RA1 (I,J) = 0.0 DEC (I,J) = 0.0 HA (I,J) = 0.0 RA (I,J) = 0.0 end if end do end do if(I.eq.1.and.J.eq.1) PED(K) = 0.0 C C Begin to composite the image C if(indx.ge.0) then do J=1,NXY do I=1,NXY if(I.eq.1.and.J.eq.1.and.KKK.eq.1) then ! Determine the Sun's location at the first image DCON = 2048.0/360.0 print*, 'to here K, aimage ',K,aimage(I,J,K),RAD, DECD Doy8 = Doys8(K) nYr = 2022 call SunNewcomb8(0,nYr,Doy8,dLngSun,dLatSun,dDisSun) ! Sun's position in ecliptic coordinates rDisSun = dDisSun ! Distance of the Earthfrom Sun rLngSun = dLngSun ! Geocentric ecliptic lng(Sun) rLatSun = dLatSun ! Geocentric ecliptic lat(Sun) = 0.0 deg write(*,'(A,2F10.3)') 'The ecliptic longitude, and latitude of the Sun is:',rLngSun,rLatSun if(iecliptic.eq.1.or.iheliographic.eq.1) then print*,'Before converted to ecliptic',RAD,DECD,rLngSun,rLatSun call ECLIPTIC_EQUATOR8(0,nYr,Doy8,rLngSun,rLatSun) write(*,'(A,2F10.3)') 'The RA, and Dec of the Sun is:', rLngSun,rLatSun rLngecl = rLngSun rLatecl = rLatSun print*,'As above the Sun Location in Geocentric',rLngecl,rLatecl call ECLIPTIC_EQUATOR8(1,nYr,Doy8,rLngecl,rLatecl) ! Geocentric Earth Equatorial to Eclipic coord(P) print*,'Sun Location in Ecliptic',rLngecl,rLatecl dRSunE = 180.0 - rLngecl RSun = rLngecl + dRSunE ! This is RA of the Sun for ecliptic coordinates RSunE = RSunE dDSun = 0.0 if(iheliographic.eq.1) then call ECLIPTIC_HELIOGRAPHIC8(0,nYr,Doy8,rLngecl,rLatecl) ! Ecliptic to Heliographic coord(P) print*,'Sun Location in Heliographic',rLngecl,rLatecl dRSunH = 180.0 - rLngecl Rsun = rLngecl + dRSunH ! overwrite RSun if iheliographic = 1 for below dDSun = 0.0 - rLatecl DSun = rLatecl + dDSun end if if(ieciptic.eq.1) print*,'Sun placed in Ecliptic at 180 deg and at Zero latitude',dRSunE,RSunE if(iheliographic.eq.1) write(*,'(A,2F9.4)')'Sun placed in Heliographic at 180 deg longitude',dRSunH,RSun if(iheliographic.eq.1) write(*,'(A,2F9.4)')'Sun placed in Heliographic at zero deg latitude',dDSun,DSun end if end if C C Before rotate or stretch images iped = 1 if(iped.eq.1.and.indx.eq.0) then ! only the unrotated and non-zeroed image will work indx = 0 if(J.eq.1.and.I.eq.1) then PEDDD = 0.0 ANPEDD = 0.0 end if if(J.le.200.and.I.le.200) then II = I JJ = J PEDDD = PEDDD + aimage(I,J,K) ANPEDD = ANPEDD + 1.0 end if if(J.eq.1824.and.I.eq.1824) then PEDD = PEDDD/ANPEDD PEDDD = 0.0 ANPEDD = 0.0 sigp = 1.0 sigm = 1.0 if(J.le.(NXY-199).and.I.le.(NXY-199)) then II = I - 1824 JJ = J - 1824 end if if(aimage(II,JJ,K).lt.(PEDD+sigp).and.aimage(II,JJ,K).gt.(PEDD-sigm)) then PEDDD = PEDDD + aimage(II,JJ,K) ANPEDD = ANPEDD + 1.0 end if if(I.eq.NXY.and.J.eq.NXY) then PEDD = PEDDD/ANPEDD PED(K) = PEDD print*, 'PEDDD, ANPEDD, PEDD(K) ', PEDDD, ANPEDD, PED(K) stop end if end if end if C End before rotate or stretch images C if(I.eq.1000.and.J.eq.1000) write(*,'(A,2I6,3F12.4)'), 'Before indx= 0 ',K,indx,aimage(I,J,K),distxc,distyc C idistok = 1 if(idistok.eq.1) then distxc = XDIST(K) distyc = YDIST(K) else distxc = distxci distyc = distyci end if AJYO = float(J)-distyc AIXO = float(I)-distxc DIS = sqrt(AIXO**2 + AJYO**2) if(DIS.le.dmax) then if(indx.eq.3.or.indx.eq.4.or.indx.eq.5) go to 999 ! No need to rotate or stretch images AJYD = AJYO/DIS if(AJYD.gt.1.0) AJYD = 1.0 ago = asind(AJYD) agoprot = ago+rotrmm(K) if(agoprot.gt.360.0) agoprot = agoprot - 360.0 if(agoprot.lt. 0.0) agoprot = agoprot + 360.0 CD = 1.0 SD = 1.0 ! OK if(float(I).lt.distxc.and.float(J).gt.distyc) then agoprot = ago-rotrmm(K) if(agoprot.gt.360.0) agoprot = agoprot - 360.0 if(agoprot.lt. 0.0) agoprot = agoprot + 360.0 CD = -1.0 SD = 1.0 end if if(float(I).le.distxc.and.float(J).le.distyc) then ! OK agoprot = ago-rotrmm(K) if(agoprot.gt.360.0) agoprot = agoprot - 360.0 if(agoprot.lt. 0.0) agoprot = agoprot + 360.0 CD = -1.0 SD = 1.0 end if Xnew = CD*DIS*cosd(agoprot) + distxc Ynew = SD*DIS*sind(agoprot) + distyc C C Begin shift the image to the pole if(indx.ge.1) then distycl = distyc distyc = distyc + (PDIST1 - PDIST(K)) end if C End shift the image to the pole C C Begin warp the image if(indx.ge.2) then distycl = distyc distyc = PDIST(K) - factor(PDIST(K))*PDIST(K) + distyc C if(J.eq.nint(distycl).and.I.eq.nint(distxc)) C & write(*,'(3(A,F8.1))'),'indx = 2, distycl= ',distycl,' PDIST= ',PDIST(K),' distyc = ',distyc,Ynewmd Xnew = CD*factor(DIS)*DIS*cosd(agoprot) + distxc Ynew = SD*factor(DIS)*DIS*sind(agoprot) + distyc C if(J.eq.nint(distycl).and.I.eq.nint(distxc)) then C write(*,'(2(A,F8.1),A,F8.5)'),'indx = 2, PDIS(K = ',PDIST(K),' FDIS= ',factor(PDIST(K))*PDIST(K) C write(*,'(2(A,F8.1),A,F8.5)'),'indx = 2, distyc = ',distyc,' PDIS= ',PDIST(K),' f(DIS)*DIS = ', C & factor(DIS)*DIS C end if end if ineedcent = 1 ! Save the values of distxc, distyc if(ineedcent.eq.1) then C write(*,'(A,2I5,2F10.3)'), 'K, inum2, distxc, distyc',K, inum2, distxc, distyc if(I.eq.1000.and.J.eq.1000) then XDIST(K) = distxc YDIST(K) = distyc C write(*,'(A,I5,2F10.3)'), 'K, YDIST(K)',XDIST(K) end if end if C if(I.eq.1000.and.J.eq.1000) write(*,'(A,2I6,3F12.4)'), 'Before Spheric1',K,indx,aimage(I,J,K),distxc,distyc C End warp the image C 999 Continue C C************************************************************************************************************************** C C Begin change the image to spherical coordinates with coordinate theta at the center of the image 3/31/2023 C if(indx.ge.3) then C C Trim the image so that it does not extend into the corral if a sidereal composite is made. Farther gets into the corral. C C Make sure distxc and distyc are the ones that have been given after the data have been rotated and stretched. if(XDIST(K).lt.0.1) then print*, 'No x or y center values have been input, STOP!!' stop else distxc = XDIST(K) distyc = YDIST(K) end if if(I.eq.1000.and.J.eq.1000.and.iecliptic) write(*,'(A,2I6,3F12.4,A)'), & 'Before Spheric2',K,indx,aimage(I,J,K),distxc,distyc,' Ecliptic' if(I.eq.1000.and.J.eq.1000.and.iheliographic) write(*,'(A,2I6,3F12.4,A)'), & 'Before Spheric2',K,indx,aimage(I,J,K),distxc,distyc,' Heliographic' if(distxc.eq.0.or.distyc.eq.0) then ! These are the new input centers print*,'There are no values for this K =',K,' stop!!' stop end if if(indx.eq.3.or.indx.eq.4) then ! This is new for the images input xnew = float(I) ynew = float(J) end if IXnew = nint(Xnew) ! Trim the new images to dmax from center JYnew = nint(Ynew) itrim = 1 if(itrim.eq.1) then distxyw = sqrt((float(IXnew)-distxc)**2 + (float(JYnew)-distyc)**2) if(distxyw.ge.dmax) then aimage(IXnew,JYnew,K) = 0.0 end if if(Ifirst.eq.0.and.Jfirst.eq.0.and.aimage(IXnew,JYnew,K).ge.0.0) then Ifirst = I Jfirst = J C print*,'Ifirst, Jfirst',Ifirst, Jfirst end if end if C DCONST = 20480.0/1800.0 RCONST = 20480.0/3600.0 C Main image distycl = distyc F90 = 10.396*90.0 ! degrees/pixel x 90deg = normalized r; Pixel dist Pole to Equator Xnewb = Xnew Ynewb = Ynew Xnewmd = Xnew-distxc Ynewmd = Ynew-distyc DISXY = SQRT((Ynewmd)**2 + (Xnewmd)**2) ! Pixel distance from the image center DISXY9 = DISXY/F90 ! Distance from the image center relative to unit distance PHI = DISXY9*90.0 ! Bourke's phi value THET = atan2d(Ynewmd,Xnewmd) ! Bourke's theta value C The two below are not needed ALONG = THET + 90.0 ! Longitude Value relative to the image center according to Bourke ALAT = atan2d(1.0,DISXY9) ! Latitude Value relative to the image center C PDEG = PDIST(K)/10.396 ! degrees to the north pole from distyc (distyc is modified in 2) C C The below is the critical step glossed over in Bourke. It converts the values of his phi and theta values C into unit vectors C AX = sind(PHI)*cosd(THET) ! The normalized X-axis (east-west) in our fisheye skymap AY = sind(PHI)*sind(THET) ! The normalized Y-axis (north-south) in our fisheye skymap AZ = cosd(PHI) ! The normalized Z-axis in our fisheye skymap Beta = PDEG ! Angle to rotate from zenith pointing to north pole (see above) C C Now rotate to pole via Euler transform C Euler rotation of the Z axis (the second Euler rotation) AXV = AY*cosd(Beta) - AZ*sind(Beta) AYV = AX ! By the way our rotation in the Euler formulation is around X rather than Y AZV = AY*sind(Beta) + AZ*cosd(Beta) C C now we have our values of Hour Angle and Declination HAD = atan2d(AXV,AYV) + 90.0 ! Hour angle of the point in the Kth image (degrees) DECD = atan2d(AZV,SQRT(AXV**2 + AYV**2)) ! Declination angle of the point in the Kth image (degrees) RAD = ST(K) - HAD ! Right Ascension of the point in the Kth image (degrees) C RA can get no less than 0 or larger than 360 degrees if(RAD.lt.0.0) RAD = RAD + 360.0 if(RAD.gt.360.0) RAD = RAD - 360.0 C if(I.eq.Ifirst.and.J.eq.Jfirst.and.KKK.eq.1) then ! Provide a mark at the Sun's location at the first image DCON = 2048.0/360.0 print*, 'to here K, aimage ',K,aimage(I,J,K),RAD, DECD Doy8 = Doys8(K) nYr = 2022 call SunNewcomb8(0,nYr,Doy8,dLngSun,dLatSun,dDisSun) ! Sun's position in ecliptic coordinates rDisSun = dDisSun ! Distance of the Earthfrom Sun rLngSun = dLngSun ! Geocentric ecliptic lng(Sun) rLatSun = dLatSun ! Geocentric ecliptic lat(Sun) = 0.0 deg write(*,'(A,2F10.3)') 'The ecliptic longitude, and latitude of the Sun mark is:',rLngSun,rLatSun end if C The below establishes Ynew and Xnew to locate pixel values into the original 2048 x 2048 pixel map - other ways are better. irect = 0 C irect = 1 if(irect.eq.1) then Ynew = (DECD+90.0)*DCONST if(Ynew.le.0.5) Ynew = 1.0 if(Ynew.ge.2048.5) Ynew = 2048.0 Xnew = RAD*RCONST if(Xnew.le.0.05) Xnew = 1.0 if(Xnew.ge.2048.5) Xnew = 2048.0 end if if(iecliptic.eq.0.and.iheliographic.eq.0) then ! this is a sidereal map (Earth Equatorial Coordinates) if(I.eq.Ifirst.and.J.eq.Jfirst) then ! Provide a mark at the Sun's location at in the first image call ECLIPTIC_EQUATOR8(0,nYr,Doy8,rLngSun,rLatSun) write(*,'(A,2F10.3)') 'The RA, and Dec of the Sun is:', rLngSun,rLatSun do IS=1,9 ! This provides a spot at the location of the Sun. do JS=1,15 rLatP = ((rLatSun+(-7.0+float(JS))/7.0)+90.0)*DCON rLngP = (rLngSun +(-7.0+float(IS))/3.5)*DCON aimagexyrk(nint(rLngP),nint(rLatP),KKK) = 4000.0 ! make this point bright in the KKK image imagexyrk (nint(rLngP),nint(rLatP),KKK) = 1 ! make sure this point is counted in the KKK image print*,aimagexyrk(nint(rLngP),nint(rLatP),KKK),KKK,nint(rLngP),nint(rLatP) end do end do end if else ! this is an sun-centered ecliptic or heliographic map if(iheliographic.eq.1) then ! this is a heliographic map (Earth Heliographic Coordinates) if(I.eq.Ifirst.and.J.eq.Jfirst) then ! Provide a mark at the Sun's location in each image do IS=1,9 ! This provides a spot at the location of the Sun. do JS=1,15 rLatP = (((-7.0+float(JS))/7.0)+90.0-dDSun)*DCON ! The sun is not at zero rLngP = (180.0 +(-7.0+float(IS))/3.5)*DCON ! The sun is assumed at the image center aimagexyrk(nint(rLngP),nint(rLatP),KKK) = 4000.0 ! make this point bright in the KKK image imagexyrk (nint(rLngP),nint(rLatP),KKK) = 1 ! make sure this point is counted in the KKK image C print*,'a point',aimagexyrk(nint(rLngP),nint(rLatP),KKK),KKK,nint(rLngP),nint(rLatP),90.0-dDSun end do end do end if end if if(iecliptic.eq.1) then ! this is an ecliptic map (Earth Ecliptic Coordinates) if(I.eq.Ifirst.and.J.eq.Jfirst) then ! Provide a mark at the Sun's location in each image do IS=1,9 ! This provides a spot at the location of the Sun. do JS=1,15 rLatP = (((-7.0+float(JS))/7.0)+90.0)*DCON ! The sun is at zero rLngP = (180.0 +(-7.0+float(IS))/3.5)*DCON ! The sun is assumed at the image center aimagexyrk(nint(rLngP),nint(rLatP),KKK) = 4000.0 ! make this point bright in the KKK image imagexyrk (nint(rLngP),nint(rLatP),KKK) = 1 ! make sure this point is counted in the KKK image C print*,'a point',aimagexyrk(nint(rLngP),nint(rLatP),KKK),KKK,nint(rLngP),nint(rLatP),90.0 end do end do end if end if call ECLIPTIC_EQUATOR8(1,nYr,Doy8,RAD,DECD) ! Geocentric Earth Equatorial to eclip. coord(P) RADDD = RAD + dRSunE DECDD = DECD if(iheliographic.eq.1) call ECLIPTIC_HELIOGRAPHIC8(0,nYr,Doy8,RAD,DECD) ! Ecliptic to Heliographic coord(P) C if(iecliptic.eq.1) RAD = RAD + dRSunE ! Shift RAD to ecliptic Sun Center C if(iheliographic.eq.1) RAD = RAD + dRSunH ! Shift RAD to heliographic Sun Center C if(RAD.gt.360.0) RAD = RAD - 360.0 if(RAD.lt.0.0) RAD = RAD + 360.0 if(I.eq.Ifirst.and.J.eq.Jfirst.and.KKK.eq.1) print*,'I,J 1000 Ecliptic RA, Dec ',RADDD,DECDD if(I.eq.Ifirst.and.J.eq.Jfirst.and.KKK.eq.1.and.iheliographic.eq.1) & print*,'I,J 1000 Heliographic RA, Dec',RAD,DECD C if(I.eq.Ifirst.and.J.eq.Jfirst.and.KKK.eq.1) stop end if ihammer = 0 C ihammer = 1 if(I.eq.Ifirst.and.J.eq.Jfirst) write(*,'(A,I2,I5,I6,4F12.4)'), 'Before ihammer', & ihammer,K,indx,aimage(I,J,K),RAD,DECD,90.0-dDSun if(ihammer.eq.1) then sqtopcc = SQRT(1.0 + cosd(DECD)*cosd(RADPM/2.0)) RADPM = RAD - 180.0 Ynew = ((((sind(DECD))/sqtopcc)*90.0) + 90.0 - dDSun)*RCONST if(Ynew.le.0.5) Ynew = 1.0 if(Ynew.ge.2048.5) Ynew = 2048.0 Xnew = ((((cosd(DECD)*sind(RADPM/2.0))/sqtopcc)*180.0) + 180.0)*RCONST if(Xnew.le.0.05) Xnew = 1.0 if(Xnew.ge.2048.5) Xnew = 2048.0 else Ynew = (DECD + 90.0)*RCONST if(Ynew.le.0.5) Ynew = 1.0 if(Ynew.ge.2048.5) Ynew = 2048.0 Xnew = RAD*RCONST if(Xnew.le.0.05) Xnew = 1.0 if(Xnew.ge.2048.5) Xnew = 2048.0 if(I.eq.Ifirst.and.J.eq.Jfirst) write(*,'(A,I2,I5,I6,7F12.4)'), 'No ihammer ', & ihammer,K,indx,aimage(I,J,K),RAD,DECD,90.0-dDSun,Xnew,Ynew,RCONST end if if(I.eq.Ifirst.and.J.eq.Jfirst) print*,' ' iprint = 0 C iprint = 1 if(iprint.eq.1) then C if(K.eq.2803) stop if(K.ge.2802.and.K.lt.2803.and.nint(Xnewmd).eq.0.and.distxyw.lt.dmax) then print*,'Mark 1 I J distxc distyc Xnewb Ynewb' write(*,'(A,2I5,2F10.3,2F10.3)'),' ',I,J,distxc,distyc,Xnewb,Ynewb write(*,'(3(A,F9.3))'), 'indx = 3a, PDIST = ',PDIST(K),' Ynewmd= ',Ynewmd,' Xnewmd= ',Xnewmd write(*,'(3(A,F9.3))'), 'indx = 3b, F90 = ',F90 ,' DISXY = ',DISXY ,' DISXY9= ',DISXY9 write(*,'(3(A,F9.3))'), 'indx = 3c, PHI = ',PHI ,' THET = ',THET ,' ALAT = ',ALAT write(*,'(3(A,F9.3))'), 'indx = 3d, AY = ',AY ,' AX = ',AX ,' AZ = ',AZ print*,' ' write(*,'(3(A,F9.3))'), 'indx = 3e, PDEG = ',PDEG,' AYCOS = ',AY*cosd(Beta),' AZSIN = ',AZ*sind(Beta) write(*,'(3(A,F9.3))'), 'indx = 3f, Beta = ',Beta,' AYSIN = ',AY*sind(Beta),' AZCOS = ',AZ*cosd(Beta) write(*,'(3(A,F9.3))'), 'indx = 3g, AYV = ',AYV ,' AXV = ',AXV ,' AZV = ',AZV write(*,'(3(A,F9.3))'), 'indx = 3h, DEC = ',DECD ,' HA = ',HAD ,' RA = ',RAD write(*,'(3(A,F9.1))'), 'indx = 3i, ST = ',ST(K) ,' Ynew = ',Ynew ,' Xnew = ',Xnew print*, ' ' print*, ' ' end if end if end if if(DIS.le.Dmax.and.K.eq.2812.and.I.eq.1000) then if(DECD.gt.DECDB) DECDB = DECD if(DECD.lt.DECDS) DECDS = DECD if(Ynew.gt.YnewB) YnewB = Ynew if(Ynew.lt.YnewS) YnewS = Ynew end if C C End change the image to spherical coordinates C IXnew = nint(Xnew) JYnew = nint(Ynew) if(ihammer.eq.0) then ! this accumulates non-hammer images into the bottom half of the 2048 image NXYM = NXY ! and this works and makes sure the images don't go beyond 1024 pixels else NXYM = NXY end if if(IXnew.ge.1.and.IXnew.le.NXY.and.JYnew.ge.1.and.JYnew.le.NXYM) then ! Added protection if dmax is too big <1 , >NYX if(IXnew.eq.875.and.JYnew.eq.569) print*, 'Point?',aimagexyrk(IXnew,JYnew,KKK),IXnew,JYnew,KKK C C Translate images to a new array and make images continuous C if(aimage(I,J,K).lt.aulim.and.aimage(I,J,K).ne.0.0) then imagexyrk (IXnew,JYnew,KKK) = imagexyrk (IXnew,JYnew,KKK) + 1 aimagexyrk(IXnew,JYnew,KKK) = aimagexyrk(IXnew,JYnew,KKK) + aimage(I,J,K) C DECC (IXnew,JYnew,K) = DECC (IXnew,JYnew,K) + DECCC (I,J,K) C HAA (IXnew,JYnew,K) = HAA (IXnew,JYnew,K) + HAAA (I,J,K) C RAA (IXnew,JYnew,K) = RAA (IXnew,JYnew,K) + RAAA (I,J,K) C if(IXnew.gt.500.and.IXnew.lt.1000) then C if(aimage(I,J,K).gt.2000.0) write(*,'(5I5,10F8.1)') J,I,K,JYnew,IXnew,aimage(I,J,K),aimagexyrk(IXnew,JYnew,K),Ynew,Xnew C end if end if end if if(JYnew.ge.998+500.and.JYnew.le.1002+500.and.IXnew.ge.1126+1000.and.IXnew.le.1128+1000) then C write(*,'(5I5,9F8.1)') J,I,K,JYnew,IXnew,ago,rotrmm(K),agoprot,aimagexyrk(IXnew,JYnew,K) C write(*,'(5I5,9I8)') J,I,K,JYnew,IXnew,imagexyrk(IXnew,JYnew,K) end if end if end do end do end if C write(*,'(A,8I6)')'Value of increments A',K, inum,inum1,inum2,inum3,inum4,inum1+inumhi,inum3-inum2-2-inumid2 if(K.eq.inum2) inum=inum3-inum2-2-inumid2 C write(*,'(A,8I6)')'Value of increments B',K, inum,inum1,inum2,inum3,inum4,inum1+inumhi,inum3-inum2-2-inumid2 if(K.ge.inum4) go to 1111 end do 1111 continue C write(*,'(A,8I6)')'Value of increments C',K, inum,inum1,inum2,inum3,inum4,inum1+inumhi,inum3-inum2-2-inumid2 C stop C C End composite the image C C do J=1,NXY C do I=1,NXY C aimagexyr(I,J) = 0.0 C end do C end do C Make image sequence continuous C inum = -1 C do KK=inum1,inum4 do KK=1,KKK do J=1,NXY ! This combines all images into one do I=1,NXY if(I.eq.1000.and.J.eq.1000.and.KK.eq.10) write(*,'(A,2I6,3F12.4)'), 'Before Image C1',KK,indx,aimagexyrk(I,J,KK), & float(imagexyrk (I,J,KK)),aimagexyr1(I,J) aimagexyr1(I,J) = aimagexyr1(I,J) + aimagexyrk(I,J,KK) imagexyr1 (I,J) = imagexyr1 (I,J) + imagexyrk (I,J,KK) if(I.eq.1000.and.J.eq.1000.and.KK.eq.10) write(*,'(A,2I6,3F12.4)'), 'Before Image C2',KK,indx,aimagexyrk(I,J,KK), & float(imagexyr1 (I,J)),aimagexyr1(I,J) end do end do if(imagexyrk(I,J,KK).ne.0) aimagexyrk(I,J,KK) = aimagexyrk(I,J,KK)/float(imagexyrk(I,J,KK)) ! This just deals with the KK image C print *,'To here 2a',inum1,inum2 if(KK.eq.KKK) then C print *,'To here 2b',KK,inum1,inum2 do J=1,NXY do I=1,NXY if(imagexyr1(I,J).ne.0) then aimagexyr1(I,J) = aimagexyr1(I,J)/float(imagexyr1(I,J)) end if end do C if(J.ge.998.and.J.le.1002) then C write(*,'(3I5,9F8.1)') J,I,KK,(aimagexyr1(I,J),I=1126,1128) C write(*,'(3I5,9I8)') J,I,KK,(imagexyr1 (I,J),I=1126,1128) C end if end do end if IT = 1000 JT = 1000 write(*,'(A,2I6,3F12.4)'), 'After image C ',KK,indx,aimagexyr1(IT,JT),float(imagexyr1(IT,JT)) if(KK.eq.KKK.and.KKK.ge.ill) then C print *,'To here 4',ill,inum4-inum1+inumhi+1 do J=1,NXY do I=1,NXY if(imagexyr1(I,J).ge.3) then ! Calculate sigma of values if more than 3 points exist sigs = 0.0 nsig = 0 do II=1,KKK sigv(II) = 0.0 if(aimagexyrk(I,J,II).gt.allim.and.aimagexyrk(I,J,II).lt.aulim) then sigv(II) = aimagexyrk(I,J,II) - aimagexyr1(I,J) sigs = sigs + sigv(II)**2 nsig = nsig +1 end if if(I.eq.1000.and.J.eq.1000.and.II.eq.10) write(*,'(A,2I6,3F12.4)'), 'After Sigma 1 ',II,indx,aimagexyrk(I,J,II),sigv(II) end do sdiv = sqrt(sigs/float(nsig-1)) nsig = 0 nout = 0 do II=1,KKK C if((sigv(II)/sdiv).ge.sig3) then if(sigv(II).ne.0.and.(sigv(II)/sdiv).ge.sig4) then aimagexys(I,J) = aimagexyrk(I,J,II) + aimagexys(I,J) aimagexyn(I,J) = imagexyrk (I,J,II) + aimagexyn(I,J) nsig = nsig + 1 if(I.eq.1000.and.J.eq.1000.and.II.eq.10) write(*,'(A,2I6,3F12.4)'), 'After Sigma 2 ',II,indx,aimagexyrk(I,J,II),sigv(II) else aimagexyo(I,J) = aimagexyr1(I,J) + aimagexyo(I,J) nout = nout + 1 if(I.eq.1000.and.J.eq.1000.and.II.eq.10) write(*,'(A,2I6,3F12.4)'), 'After Sigma 3 ',II,indx,aimagexyrk(I,J,II),sigv(II) end if end do if(nsig.ne.0) then aimagexyr(I,J) = aimagexys(I,J)/float(nsig) DEC (I,J) = DEC (I,J)/float(nsig) HA (I,J) = HA (I,J)/float(nsig) RA (I,J) = RA (I,J)/float(nsig) else aimagexyr(I,J) = 0.0 end if if(nout.ne.0) then aimagexyo(I,J) = aimagexyo(I,J)/float(nout) else aimagexyo(I,J) = 0.0 if(I.eq.1000.and.J.eq.1000.and.II.eq.10) write(*,'(A,2I6,3F12.4)'), 'After Sigma 4 ',II,indx,aimagexyo(I,J),sigv(II) end if else num = 0 do II=1,KKK if(aimagexyrk(I,J,II).gt.allim.and.aimagexyrk(I,J,II).lt.aulim) then aimagexyr(I,J) = aimagexyrk(I,J,II) + aimagexyr(I,J) aimagexyn(I,J) = imagexyrk (I,J,II) + aimagexyn(I,J) num = num + 1 end if if(I.eq.1000.and.J.eq.1000.and.III.eq.10) write(*,'(A,2I6,3F12.4)'), 'After Sigma 5 ',II,indx,aimagexyrk(I,J,II) if(I.eq.1000.and.J.eq.1000.and.III.eq.10) print*,' ' end do 2225 continue C aimagexyn(I,J) = imagexyrk (I,J,K) if(num.ne.0) aimagexyr(I,J) = aimagexyr(I,J)/float(num) aimagexyo(I,J) = 0.0 end if end do end do else if(KK.eq.KKK) then write(*,'(A,I5,A,I5,A)') 'There are only ', KKK,' images of ill =',ill, ' . Combine with all pixels.' do J=1,NXY do I=1,NXY num = 0 do II=1,KKK if(aimagexyrk(I,J,II).gt.allim.and.aimagexyrk(I,J,II).lt.aulim) then aimagexyr(I,J) = aimagexyrk(I,J,II) + aimagexyr(I,J) aimagexyn(I,J) = imagexyrk (I,J,II) + aimagexyn(I,J) num = num + 1 end if end do if(num.ne.0) aimagexyr(I,J) = aimagexyr(I,J)/float(num) aimagexyo(I,J) = 0.0 ! no sigma values used if(I.eq.1000.and.J.eq.1000.and.II.eq.10) write(*,'(A,2I6,3F12.4)'), 'After Sigma 6 ',II,indx,aimagexyrk(I,J,II) if(I.eq.1000.and.J.eq.1000.and.II.eq.10) print*,' ' end do end do end if end if end do C stop C C Begin background image subtraction ibak = 0 if(ibak.eq.1) then if(indx.le.3) then ! Safety since the backgrounds are only at index 2 do J=1,NXY do I=1,NXY aimageb = aimagebak(I,J) distxyw = sqrt((float(I)-distxc)**2 + (float(J)-distyc)**2) if(distxyw.le.dmax) then C if(IXNew.eq.1050.and.JYnew.lt.1005.and.JYnew.gt.995) print*, I,J,aimageb, distxyw, bakmin if(aimageb.lt.bakmin/10.0) aimageb = bakmin C aimagexyr(I,J) = aimagexyr(I,J) - aimageb aimagexyr(I,J) = aimagexyr(I,J) - aimageb + bakmin end if end do end do end if end if C End background image subtraction C C print*, 'to here 5' if(indx.eq.3) then aimagel = 0.0 aimagnl = 0.0 do J=1,NXY do I=1,NXY distxy = sqrt((float(I)-distxc)**2 + (float(J)-distyc)**2) if(distxy.le.dmax) then if(aimagexyr(I,J).gt.aimagel) aimagel = aimagexyr(I,J) if(float(imagexyr1(I,J)).gt.aimagnl) aimagnl = float(imagexyr1(I,J)) end if end do end do print*,' ' write(*,'(A,F10.1,A,F6.1)') ' The largest image value is', aimagel,'; The greatest radius mapped =',dmax write(*,'(A,F10.1,A,F6.1)') ' The largest img # value is', aimagnl end if print *, ' ' print*, 'Out of rot_image_2ra_dec.f' print *, ' ' return end