C+ C NAME: C rot_image_2pole.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,inum,NOBS,centx,centy,dmax,rotrmm,PDIST,ST,aimage,aimagexyr,aimagexyn,aimagexyo) 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) C indx = 1 indx 0 + shift index in y so pole is at the a given pixel location C indx = 2 indx 0 + 1 + warp the image relative to the y center pixel location C indx = 3 indx 0 + 1 + 2 + change the image to spherical coordinates C indx = 4 indx 0 + 1 + 2 + 3 + rotate the image around the pole according to the new sidereal time 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 November-2022, Bernard Jackson (UCSD) C- subroutine rot_image_2pole(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 image & PED (NOBS) ! The Pedestal during the image taking 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(inum2-inum2+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 ill = 5 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 print*,' ' write(*,'(3(A,I5))') 'Into rot_image_2pole.f, inum1 =',inum1,', inum2 =',inum2,', ill =',ill write(*,'(3(A,I5))') 'Into rot_image_2pole.f, inum3 =',inum3,', inum4 =',inum4 print*, ' ' YDIS = 0.0 ANYDIS = 0.0 KKKK = 0 ! Trigger to provide a sun-centered point in the analysis inum = -1 do KK=inum1,inum4 K = KK + inum + 1 distxci = centx(K) C distxc = centx(K) distyci = centy(K) do J=1,NXY do I=1,NXY aimagexyrk(I,J,K) = 0.0 imagexyrk (I,J,K) = 0 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(K.eq.inum1) then aimagexyr (I,J) = 0.0 aimagexyr1(I,J) = 0.0 aimagexys (I,J) = 0.0 imagexyr1 (I,J) = 0.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 C C Begin to composite the image C if(indx.ge.0) then do J=1,NXY do I=1,NXY distxc = distxci distyc = distyci AJYO = float(J)-distyc AIXO = float(I)-distxc DIS = sqrt(AIXO**2 + AJYO**2) if(DIS.le.dmax) then 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)) if(J.eq.nint(distycl).and.I.eq.nint(distxc)) then C print*, ' ' print*,I,J,'Image =',K C write(*,'(3(A,F8.1))'),'indx = 1, distycl= ',distycl,' PDIST= ',PDIST(K),' distyc = ',distyc end if 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 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 C End warp the image C C Begin background image subtraction ibak = 0 if(ibak.eq.1) then IXnew = nint(Xnew) JYnew = nint(Ynew) aimageb = aimagebak(IXnew,JYnew) distxyw = sqrt((float(IXnew)-distxc)**2 + (float(JYnew)-distyc)**2) if(distxyw.lt.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 aimage(IXnew,JYnew,K) = aimage(IXnew,JYnew,K) - aimageb aimage(IXnew,JYnew,K) = aimage(IXnew,JYnew,K) - aimageb + bakmin end if end if C End background image subtraction 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. IXnew = nint(Xnew) JYnew = nint(Ynew) itrim = 1 if(itrim.eq.1) then if(indx.le.2) then distxyw = sqrt((float(IXnew)-distxc)**2 + (float(JYnew)-distyc)**2) if(distxyw.ge.dmax) then aimage(IXnew,JYnew,K) = 0.0 end if end if end if C C Begin change the image to spherical coordinates with coordinate theta at the center of the image 3/31/2023 C DCONST = 20480.0/1800.0 RCONST = 20480.0/3600.0 if(indx.ge.3) then C Main image if(distxyw.lt.dmax.and.KKKK.eq.0) KKKK = 1 distycl = distyc F90 = 10.396*90.0 ! degrees/pixel x 90deg = normalized r; Pixel dist Pole to Eq. 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 from 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 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 C iecliptic = 0 iecliptic = 1 iheliographic = 0 C iheliographic = 1 if(iecliptic.eq.1) then Doy8 = Doys8(K) nYr = 2022 call SunNewcomb8(0,nYr,Doy8,dLngSun,dLatSun,dDisSun) ! Sun's position in ecliptic coordinates rDisSun = dDisSun ! Geocentric eclip. coord(Sun) rLngSun = dLngSun ! Geocentric eclip. lng(Sun) rLng = dLngSun rLat = dLatSun ! Geocentric ecliptic lat(Sun) C print*,'to here 2', rlng,rlat call ECLIPTIC_EQUATOR8(0,nYr,Doy8,rLng,rLat) dRA = rLng ! RA (Sun) dRA = RAD-dRA ! RA(P)-RA(Sun) if (abs(dRA) .gt. 180.) dRA = dRA-sign(1.,dRA)*360. C print*,'to here 3', dRA call ECLIPTIC_EQUATOR8(1,nYr,Doy8,RAD,DECD) ! Geocentric eclip. coord(P) RAD = RAD-rLngSun rLngP = RAD ! Geocentric ecliptic lng(P)-lng(Sun) rLatP = DECD ! Geocentric ecliptic lat(P) RAD = rLngSun+180.+RAD ! Heliocentric eclip. lng(P) if(KKKK.eq.1) then ! This provides a point at the location of the Sun. do IS=1,5 do JS=1,5 rLngP = rLngSun-2+IS rLatP = rLatP-2+JS aimagexyr(rLatP,rLngP) = 4000.0 end do end do end if if(iecliptic.eq.1.and.iheliographic.eq.0) then Xnew = rLngP ! Geocentric ecliptic lng(P)-lng(Sun) Ynew = rLatP ! Geocentric ecliptic lat(P) end if C print*,'to here 4', rLngP,rLatP if(iheliographic.eq.1) then call ECLIPTIC_HELIOGRAPHIC8(0,nYr,Doy8,RAD,DECD) xlon = RAD ! Heliographic long(P) xlat = DECD ! Heliographic lat(P) Xnew = xlon Ynew = xlat C print*,'to here 5', xlon,xLat end if end if ihammer = 0 C ihammer = 1 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)*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 end if 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(IXnew.ge.1.and.IXnew.le.NXY.and.JYnew.ge.1.and.JYnew.le.NXY) then ! Added protection if dmax is too big if(K.gt.4449.and.I.eq.999.and.J.eq.1000) then end if if(aimage(I,J,K).lt.aulim.and.aimage(I,J,K).ne.0.0) then imagexyrk (IXnew,JYnew,K) = imagexyrk (IXnew,JYnew,K) + 1 aimagexyrk(IXnew,JYnew,K) = aimagexyrk(IXnew,JYnew,K) + 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 if(K.eq.inum2) inum=inum3-inum2-2 if(K.eq.inum4) go to 1111 end do 1111 continue C stop C C End composite the image C inum = -1 do KK=inum1,inum4 K = KK + inum + 1 distxci = centx(K) C distxc = centx(K) distyci = centy(K) do J=1,NXY ! This combines all images into one do I=1,NXY aimagexyr1(I,J) = aimagexyr1(I,J) + aimagexyrk(I,J,K) imagexyr1 (I,J) = imagexyr1 (I,J) + imagexyrk (I,J,K) end do end do if(imagexyrk(I,J,K).ne.0) aimagexyrk(I,J,K) = aimagexyrk(I,J,K)/float(imagexyrk(I,J,K)) ! This just deals with the K image C print *,'To here 2a',inum1,inum2 if(K.eq.inum4) then C print *,'To here 2b',K,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,K,(aimagexyr1(I,J),I=1126,1128) C write(*,'(3I5,9I8)') J,I,K,(imagexyr1 (I,J),I=1126,1128) C end if end do end if if(K.eq.inum4.and.(inum4-inum1+1).ge.ill) then C print *,'To here 4',ill,inum4-inum1+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 iinum = -1 do III=inum1,inum4 II = iinum + III + 1 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(II.eq.inum2) iinum=inum3-inum2-2 if(II.eq.inum4) go to 2223 end do 2223 continue sdiv = sqrt(sigs/float(nsig-1)) nsig = 0 nout = 0 iinum = -1 do III=inum1,inum4 II = iinum + III + 1 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(II.eq.inum2) iinum=inum3-inum2-2 if(II.eq.inum4) go to 2224 else aimagexyo(I,J) = aimagexyr1(I,J) + aimagexyo(I,J) nout = nout + 1 end if 2224 continue 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 end if else num = 0 iinum = -1 do III=inum1,inum4 II = iinum + III + 1 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(II.eq.inum2) iinum=inum3-inum2-2 if(II.eq.inum4) go to 2225 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(K.eq.inum4) then write(*,'(A,I5,A)') 'There are only ', inum2-inum1+1,' images. Combine with all pixels' do J=1,NXY do I=1,NXY num = 0 iinum = -1 do III=inum1,inum4 II = iinum + III + 1 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(II.eq.inum2) iinum=inum3-inum2-2 if(II.eq.inum4) go to 2226 end do 2226 continue 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 do end do end if end if if(K.eq.inum2) inum=inum3-inum2-2 if(K.eq.inum4) go to 2222 end do ! end K=num1,num2 2222 continue C stop C C Begin background image subtraction ibak = 1 if(ibak.eq.1) then if(indx.eq.2) 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.lt.4) 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 C aimagexyr(I,J) = 0.0 C aimagexyo(I,J) = 0.0 C else 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_2pole.f' return end