C+ C NAME: C program ashi_HR_maker C PURPOSE: C To read the times from the ashi data sheet to determine the latitude an longitude of the ASHI balloon, C and then interpolate these values to obtain the values at the time of the ASHI image data. Once the C longitude values are obtained, they can be used to determine the sidereal time of each image taken. C The sidereal times and RA of each star will then give the Hour Angle values of each stellar measurement C determined, and these can then be used to locate and rotate each image to provide the same stellar map. C All stars with an RA and Dec known can then be located on this map. C C CALLING SEQUENCE: C ./ashi_hr_maker C INPUTS: C iyr integer year of observation C idoyst integer starting day of the observations C NLMAX integer maximum number of track values allowed C NOBSM integer maximum number of image observations C NSRC integer maximum numbers of RA and Dec source positions used C C OUTPUTS: C alat(I) real latitude of ASHI C alon(I) real longitude of ASHI C jd C C FUNCTIONS/SUBROUTINES: C read_ashi_track, poly_interpolate, sidereal, C C MODIFICATION HISTORY: C September 2022, Bernard Jackson (UCSD) C- program ashi_hr_maker parameter (NLMAX = 10000) ! Maximum number of Number track locations parameter (NSTAR = 30) ! Number of "Stars" tracked parameter (NOBSM = 5000) ! Maximum number of images with different times C to analyze parameter (NXYM = 40) parameter (NVAl = 41) character cfilei*15 /'ashiimages1.txt'/ character cfile101*24 character cfile1*27 /'Polaris_______centroids.txt'/, & cfile2*27 /'B_Capricorn___centroids.txt'/, ! New & cfile3*27 /'E_Pegasus_____centroids.txt'/, ! New & cfile4*27 /'G_Draco_______centroids.txt'/, & cfile5*27 /'Saturn________centroids.txt'/, & cfile6*27 /'D_Capricornus_centroids.txt'/ ! Bad character cfilen*15 /'722NTflight.rpt'/ real alat (NLMAX), & alon (NLMAX), & RA (NSTAR), & DEC (NSTAR), & vmag (NSTAR), & XS1 (NOBSM), & YS1 (NOBSM), & XS11 (NOBSM), & YS11 (NOBSM), & XS12 (NOBSM), & YS12 (NOBSM), & XS13 (NOBSM), & YS13 (NOBSM), & XS14 (NOBSM), & YS14 (NOBSM), & XS15 (NOBSM), & YS15 (NOBSM), & XS16 (NOBSM), & YS16 (NOBSM), & ST3 (NOBSM), & HA (NSTAR,NOBSM), & SDISTP (NSTAR,NOBSM), & SDISTM (NSTAR,NOBSM), & SDISTH (NSTAR,NOBSM), & SDISTMZ(NSTAR,NOBSM), & SDISTZ (NSTAR,NOBSM), & XYDISTZ(NSTAR,NOBSM), & DISDECC(NSTAR,NOBSM), & DISPXYC(NSTAR,NOBSM), & SMD (NSTAR,NOBSM), & ROTA (NSTAR,NOBSM), & AROTA (NSTAR,NOBSM), & xOPpos (NSTAR,NOBSM), & yOPpos (NSTAR,NOBSM), & ROTAP (NOBSM), & scalestah (NSTAR), & scalestav (NSTAR), & scalesta (NSTAR), & DISCIJRSAS (NSTAR), & DCENTRORSAS(NSTAR) real DZ (NSTAR,NOBSM), ! distance to pole from the star along the celestial meridian & DM (NSTAR,NOBSM), ! distance to the celestial meridian from the star & DZP (NSTAR,NOBSM), ! distance to pole from the star along the celestial meridian in pixels & DMP (NSTAR,NOBSM), ! distance to the celestial meridian from the star in pixels & DC (NSTAR,NOBSM), ! distance from the latitude of ASHI to the pole in degrees & DCP (NSTAR,NOBSM), ! distance from the latitude of ASHI to the pole in pixels & DCZX (NSTAR,NOBSM), ! location of X at the zenith in degrees. The image center & DCZY (NSTAR,NOBSM), ! location of Y at the zenith in degrees. The image center & DCZXP(NSTAR,NOBSM), ! location of X at the zenith in pixels. The image center & DCZYP(NSTAR,NOBSM), ! location of X at the zenith in pixels. The image center & XSBRO(NSTAR,NOBSM), ! X pole star position before rotation & YSBRO(NSTAR,NOBSM), ! Y pole star position before rotation & PSXY (NSTAR,NOBSM), ! pole star position from X, Y before rotation & XCENPO (NSTAR,NOBSM), ! location of the original X center of the image in degrees from the pole & YCENPO (NSTAR,NOBSM), ! location of the original Y center of the image in degrees from the pole & XCENPOP(NSTAR,NOBSM), ! location of the original X center of the image in degrees from the pole in pixels & YCENPOP(NSTAR,NOBSM), ! location of the original Y center of the image in degrees from the pole in pixels & DXCENPOP(NSTAR,NOBSM), ! difference in location of the original X center of the image in degrees from the pole in pixels & DYCENPOP(NSTAR,NOBSM), ! difference in location of the original Y center of the image in degrees from the pole in pixels & XCEN (NSTAR,NOBSM), ! location of the X center of the image in degrees & YCEN (NSTAR,NOBSM), ! location of the Y center of the image in degrees & DXCEN(NSTAR,NOBSM), ! change in location of the X center of the image in degrees from image 1 & DYCEN(NSTAR,NOBSM), ! change in location of the Y center of the image in degrees from image 1 & XCENP(NSTAR,NOBSM), ! location of the X center of the image in pixels & DZZ (NSTAR,NOBSM), ! theoretical Y distance from this star to the image center along the central meridian & DZZP (NSTAR,NOBSM), ! theoretical Y distance from this star to the image center along the central meridian in pixels & DZZPT(NSTAR,NOBSM), ! theoretical total distance from this star to the image center in pixels & YCENP(NSTAR,NOBSM), ! location of the Y center of the image in pixels & CMP(NSTAR,NOBSM), ! modeled total distance of the star to the original X, Y image center & CTP(NSTAR,NOBSM), ! theoretical total distance of the star to original X, Y image center & DXCENP(NSTAR,NOBSM), ! change in location of the X center of the image in pixels from image 1 & DYCENP(NSTAR,NOBSM), ! change in location of the Y center of the image in pixels from image 1 & XS1MCP(NSTAR,NOBSM), ! New X center location & YS1MCP(NSTAR,NOBSM), ! New Y center location & XS1MCOP(NSTAR,NOBSM), ! Old X center location & YS1MCOP(NSTAR,NOBSM), ! Old Y center location & XS1rot (NSTAR,NOBSM), ! Rotated X centroid location & YS1rot (NSTAR,NOBSM), ! Rotated Y centroid location & DCENTROR(NSTAR,NOBSM) real OPPLX(NSTAR,NOBSM), ! Original pole pixel X location before rotation & OPPLY(NSTAR,NOBSM), ! Original pole pixel X location before rotation & OCPX (NSTAR,NOBSM), ! Original position of the X Lat-Lon center & OCPY (NSTAR,NOBSM) ! Original position of the Y Lat-Lon center C & OPCXP(NSTAR,NOBSM), ! The X position of the lat-lon center on the image I C & OPCYP(NSTAR,NOBSM) ! The Y position of the lat-lon center on the image I real HAIJ (2000,2000), ! Hour angle of original XII, YJJ image pixel & RAIJ (2000,2000), ! RA of XII, YJJ image pixel & DECIJ (2000,2000), ! Dec of XII, YJJ image pixel & DISIJR (2000,2000), ! Distance to the Pole from the XII, YJJ image pixel & DISDECR (2000,2000), ! Distance to the center from the Pole & DISCIJR (2000,2000), ! Distance to the center from the XII, YJJ image pixel & XII (2000), & YJJ (2000) real XSMAX (NXYM,NVAL), & YSMAX (NXYM,NVAL), & XSMIN (NXYM,NVAL), & YSMIN (NXYM,NVAL) integer NIMAX (NXYM,NVAL), & NIMIX (NXYM,NVAL), & NIMAY (NXYM,NVAL), & NIMIY (NXYM,NVAL) real distxcmax(NOBSM), & distxcmin(NOBSM), & distycmax(NOBSM), & distycmin(NOBSM) real alatyy(NOBSM), & alonxx(NOBSM) integer NII(NSTAR), & NJJ(NSTAR) real XXII (NSTAR,500,500), & YYII (NSTAR,500,500), & RAIJS (NSTAR,500,500), & DECIJS (NSTAR,500,500), & DISIJRS (NSTAR,500,500), & DISDECRS (NSTAR,500,500), & DISCIJRS (NSTAR,500,500), & DCENTRORS (NSTAR,500,500) real*8 JD (NLMAX), ! Julian Dates of each tracked ASHI position & JDX (NOBSM), ! Julian date of the image from the image file & JDXS1 (NOBSM), ! Julian date of the image from the image file that includes stellar x y centroid positions & JDXSD (NOBSM), ! Julian date of the image from the image file that includes the dummy stellar x y centroid positions & JDbeg /2459817.58333333333d0/, ! Julian date at 8 pm local Ft Sumner time Aug 25 or 2 UTC Aug 26 & doyutveq8,JDio,JEpoch,doyveq8,stbeg real*8 doyck8,JDioo ! Sidereal time tests character cMon*3 ! Sidereal time tests iyr = 2022 ! year of ASHI observation idoyst = 237 ! starting DOY of ASHI observation xcentpx = 1126.5 ! xcenter in pixels as determined by Polaris rotation ycentpy = 994.5 ! ycenter in pixels as determined by Polaris rotation distcx = 1126.5 ! center of the image in X as determined by Matthew distcy = 1000.5 ! center of the image in Y as determined by Matthew xscale = 568.5/54.73 ! xscale (10.387) of image in pixels from declination of Polaris yscale = 569.5/54.73 ! yscale (10.406) of image in pixels from declination of Polaris xscaleu = 10.396 ! xscale of image used until we know better yscaleu = 10.396 ! yscale of image used until we know better call read_ashi_track(cfilen,NLMAX,NL,iyr,idoyst,JD,alat,alon) ! Determine Julian dates and alat and alon values along the ASHI track call readimagetimes(cfilei,NOBSM,NOBS,JDX) ! read an ashi file to determine the output JD UTC of each image. iyrveq = 2022 ! year of the vernal equinox doyutveq8 = 79.0d0 + 15.0d0/24.0d0 + 33.0d0/1440.0d0 ! doy of the vernal equinox (2022 March 20 at 15:33 UTC) C C Stars: C from https://en.wikipedia.org/w/index.php?title=Category:Bright_Star_Catalogue_objects&pagefrom=8217%0A35+Vulpeculae#mw-pages C (Coordinates in year 2000 values) C Saturn: C from https://in-the-sky.org/ephemeris.php?ird=1&objtype=1&objpl=Saturn&objtxt=Saturn&tz=1&startday=24&startmonth=8&startyear=2022&interval=1&rows=25 C (Coordinates in year 2000 values) C RA(1) = 2.0*15.0 + 31.0/4.0 + 49.09/240.0 ! Polaris RA in degrees DEC(1) = 89.0 + 15.0/60.0 + 50.8/3600.0 ! Polaris DEC in degrees scalestah(1) = 10.396 ! Polaris star horizontal scale in pixels/deg scalestav(1) = 10.396 ! Polaris star vertical scale in pixels/deg C scalesta (1) = 10.396 ! Polaris star scale in pixels/deg (before limits were 10.378 to 10.406) C scalesta (1) = 10.4218 ! Polaris star scale in pixels/deg (attempts to make both NS and EW asymmetry equal) C scalesta (1) = 10.47739 (most work with this) ! Polaris star scale in pixels/deg (makes NS asymetry equal at 1.421 pixels scalesta (1) = 10.5466 ! Polaris star scale in pixels/deg (makes NS asymetry equal at 1.421 pixels to north) C With 10.4198 the east -west discrepancy is east by 0.7005 pixels, but the scale is too big by about 1.715 pixels vmag(1) = 1.98 ! variable (1.86-2.13) RA(2) = 20.0*15.0 + 21.0/4.0 + 00.7/240.0 ! Beta Capricorn RA in degrees ! New DEC(2) = -14.0 - 46.0/60.0 - 53.00/3600.0 ! Beta Capricorn Dec in degrees C scalesta(2) = 10.396 ! Beta Capricorn scale in pixels/deg scalesta (2) = 10.5466 ! Beta Capricorn scale in pixels/deg vmag(2) = 2.05 ! variable (2.01 - 2.10) RA(3) = 21.0*15.0 + 44.0/4.0 + 11.15614/240.0 ! Epsilon Pegasus RA in degrees ! New DEC(3) = 09.0 + 52.0/60.0 + 30.0311/3600.0 ! Epsilon Pegasus Dec in degrees C scalesta(3) = 10.396 ! Epsilon Pegasus scale in pixels/deg scalesta (1) = 10.5466 ! Epsilon Pegasus scale in pixels/deg vmag(3) = 2.357 ! not variable RA(4) = 17.0*15.0 + 56.0/4.0 + 36.36988/240.0 ! Gamma Draco RA in degrees DEC(4) = 51.0 + 29.0/60.0 - 20.0242/3600.0 ! Gamma Draco Dec in degrees C scalesta(5) = 10.396 ! Gamma Draco scale in pixels/deg scalesta (4) = 10.5466 ! Gamma Draco scale in pixels/deg vmag(4) = 2.23 ! not variable K5 RA(5) = 21.0*15.0 + 34.0/4.0 + 12.0/240.0 ! Saturn RA in degrees (within a second or so at midnight) DEC(5) = -15.0 - 49.0/60.0 - 5.0/3600.0 ! Saturn Dec in degrees (within a few seconds at midnight) C scalesta(6) = 10.396 ! Saturn scale in pixels/deg scalesta (5) = 10.5466 ! Saturn scale in pixels/deg vmag(5) = 0.3 ! approximate RA(6) = 5.0*15.0 + 16.0/4.0 + 41.35871/240.0 ! Capella RA in degrees DEC(6) = 45.0 + 59.0/60.0 + 52.7693/3600.0 ! Capella Dec in degrees Dec in degrees C scalesta(6) = 10.396 ! Capella scale in pixels/deg scalesta (6) = 10.5466 ! Capella scale in pixels/deg vmag(6) = -0.3 ! OK RA(7) = 0.0*15.0 + 43.0/4.0 + 35.37090/240.0 ! Beta Cetus RA in degrees DEC(7) = -17.0 - 59.0/60.0 - 11.7827/3600.0 ! Beta Cetus RA in degrees Dec in degrees C scalesta(7) = 10.396 ! Beta Cetus scale in pixels/deg scalesta (7) = 10.5466 ! Beta Cetus scale in pixels/deg vmag(7) = 1.01 ! B-V Mag RA(10) = 21.0*15.0 + 47.0/4.0 + 2.44424/240.0 ! Delta Capricornus RA in degrees ! Bad DEC(10) = -16.0 - 07.0/60.0 - 38.2335/3600.0 ! Delta Capricornus Dec in degrees C scalesta(3) = 10.396 ! Delta Capricornus scale in pixels/deg scalesta (10) = 10.5466 ! Delta Capricornus scale in pixels/deg vmag(10) = 2.81 ! eclipsing binary 1.022 day variable (-2.90 -3.05) RA (11) = 320.59 DEC(11) = -16.72 scalesta (11) = 10.5466 vmag(11) = 0.0 RA (12) = 310.35 DEC(12) = 45.28 scalesta (12) = 10.5466 vmag(12) = 0.0 RA (13) = 297.69 DEC(13) = 8.86 scalesta (13) = 10.5466 vmag(13) = 0.0 RA (14) = 305.55 DEC(14) = 40.25 scalesta (14) = 10.5466 vmag(14) = 0.0 RA (15) = 279.23 DEC(15) = 38.78 scalesta (15) = 10.5466 vmag(15) = 0.0 RA (16) = 40.41 DEC(16) = 89.26 scalesta (16) = 10.5466 vmag(16) = 0.0 RA (17) = 269.14 DEC(17) = 51.48 scalesta (17) = 10.5466 vmag(17) = 0.0 RA (18) = 222.67 DEC(18) = 74.15 scalesta (18) = 10.5466 vmag(18) = 0.0 RA (19) = 247.35 DEC(19) = -26.43 scalesta (19) = 10.5466 vmag(19) = 0.0 RA (20) = 165.92 DEC(20) = 61.75 scalesta (20) = 10.5466 vmag(20) = 0.0 RA (21) = 193.50 DEC(21) = 55.95 scalesta (21) = 10.5466 vmag(21) = 0.0 RA (22) = 206.88 DEC(22) = 49.31 scalesta (22) = 10.5466 vmag(22) = 0.0 RA (23) = 213.91 DEC(23) = 19.18 scalesta (23) = 10.5466 vmag(23) = 0.0 call interpolate_latlon(JDbeg,JD,alat,alon,NL,alatx,alonb) ! interpolate the first beginning JDbeg track value to determine the ASHI beginning longitude ialonb = alonb alonbm = alonb*60.0 - float(ialonb*60) ilonbm = alonbm alonbs = alonbm*60.0 - float(ilonbm*60) write(*,'(A,F16.7,A,F10.5,A,2I3,F6.2)') 'Beginning JDbeg = ',JDbeg,' has a longitude of',alonb,' or ',ialonb,ilonbm,alonbs C stbeg = 113.5417415874d0 stbeg = 255.0737582540d0 ! Starts at Julian date at 8 pm local Ft Sumner time Aug 25 or 2 UTC Aug 26 write(*,'(A)') 'The start true sidereal time at this location from:' write(*,'(A)') 'http://neoprogrammics.com/sidereal_time_calculator/index.php is:' write(*,'(F12.7,A)') stbeg,' degrees' cfile101 = 'all_13_centroids_500.txt' ! Matthew's 13 star centroid list for 13 stars of rotated image 500 Nstarb = 11 ! beginning of the image 500 star list Nstart = 23 ! end of the image 500 star list ixoffset = 121 ! Matthew's X offset of the image 500 star list placed at 122 jyoffset = 6 ! Matthew's Y offset of the image 500 star list placed at 7 call read_stellar_centroids(cfile101,NSTAR,NOBSM,Nstarb,Nstart,ixoffset,jyoffset,imagen,XS1rot,YS1rot) C do J=1,NSTAR ! 6 stars do J=1,1 ! Only Polaris do K=1,NOBM XS1(K) = -999.999 YS1(K) = -999.999 end do if(J.eq.1) call read_image_t_x_y(cfile1,NOBSM,NOBS,JDXS1,XS1,YS1) ! read an ashi file that gives time and x, y of a star if(J.eq.2) call read_image_t_x_y(cfile2,NOBSM,NOBS,JDXS1,XS1,YS1) ! read an ashi file that gives time and x, y of a star if(J.eq.3) call read_image_t_x_y(cfile3,NOBSM,NOBS,JDXS1,XS1,YS1) ! read an ashi file that gives time and x, y of a star if(J.eq.4) call read_image_t_x_y(cfile4,NOBSM,NOBS,JDXS1,XS1,YS1) ! read an ashi file that gives time and x, y of a star if(J.eq.5) call read_image_t_x_y(cfile5,NOBSM,NOBS,JDXS1,XS1,YS1) ! read an ashi file that gives time and x, y of a star if(J.eq.6) call read_image_t_x_y(cfile6,NOBSM,NOBS,JDXS1,XS1,YS1) ! read an ashi file that gives time and x, y of a star C call find_min_max(XS1,YS1,NOBS,XSMAX,NIMAX,XSMIN,NIMIX,YSMAX,NIMAY,YSMIN,NIMIY,NXYM,NVAL,distcx,distcy) C Nbeg = 2499 C Nend = 3552 Nbeg = 1 do I=1,NOBS C do I=Nbeg,Nend C print*,' ' print*,' ' print*,' Begin I loop: I =', I call interpolate_latlon(JDX(I),JD,alat,alon,NL,alaty,alonx) ! Interpolate at the JDX(I) value to find the latitude longitude of ASHI iprint_lat_lon = 0 ip_star_pole_st_dist = 1 if(ip_lat_lon.eq.1.or.ip_star_pole_st_dist.eq.1) then alatyy(I) = alaty alonxx(I) = alonx end if C call Julian8(0,iyrveq,doyutveq8,JDio,JEpoch) C doyveq8 = JDio C call sidereal(JDX(I),doyveq8,alonx,st) call sidereal_t(JDbeg,JDbeg,stbeg,alonb,alonb,st) ! Find the sidereal time at the JDX(I) time value, and longitude alonx of ASHI print*,' ' ist = st/15.0 aist = st/15.0 - float(ist) istm = aist*60.0 asts = (aist*60.0 - float(istm))*60.0 write(*,'(A,F10.4,I3,I3,F5.1)') 'I,Sidereal Time at ASHI beginning = degrees, hr, mi, sec',st,ist,istm,asts call sidereal_t(JDX(I),JDbeg,stbeg,alonx,alonb,st) ! Find the sidereal time at the JDX(I) time value, and longitude alonx of ASHI if(I.eq.1) st1 = st ideglo = alonb aminlo = (alonb - ideglo)*60.0 iminlo = aminlo seclo = (aminlo - iminlo)*60.0 write(*,'(A,I5,F10.3,5X,2I3,F6.2)'),'I,beginning long,ideg,imin,sec', I, alonb, ideglo, iminlo, seclo AJDsthu = (JDbeg - JDbeg)*24 + 2.0 IAJDsthu = AJDsthu thum = (AJDsthu - IAJDsthu)*60.0 ithum = thum tsecu = (thum - ithum)*60.0 IAJDsthl = IAJDsthu - 6 if(IAJDsthl.lt.0) IAJDsthl = IAJDsthl + 12 write(*,'(A,I5,F17.8,I3,I3,F5.1)') 'I,JDbeg - hr, mi, sec ',I,JDbeg,IAJDsthu,ithum,tsecu print *,' ' ist = st/15.0 aist = st/15.0 - float(ist) istm = aist*60.0 asts = (aist*60.0 - float(istm))*60.0 write(*,'(A,I5,F10.4,I3,I3,F5.1)') 'I, Sidereal Time at ASHI = degrees, hr, mi, sec',I,st,ist,istm,asts write(*,'(A,I5,F10.4)') 'I, Latitude (Declination) at ASHI in degrees = ',I,alaty C if(I.eq.I1.or.I.eq.I2) print*, 'JDX(I), st, RA(I)',JDX(I),st, RA(J) if(ST.lt.260.0) ST = ST + 360.0 XCEN(J,I) = ST YCEN(J,I) = alaty AJDsthu = (JDX(I) - JDbeg)*24 + 2.0 IAJDsthu = AJDsthu thum = (AJDsthu - IAJDsthu)*60.0 ithum = thum tsecu = (thum - ithum)*60.0 IAJDsthl = IAJDsthu - 6 if(IAJDsthl.lt.0) IAJDsthl = IAJDsthl + 12 write(*,'(A,I5,F17.8,I3,I3,F5.1)') 'I,UT time, hr, mi, sec',I,JDX(I),IAJDsthu,ithum,tsecu C print*,'here 00',JDbeg,JDX(I),stbeg,'JDbeg,JDX(I),stbeg' C print*,'here 0 ',ST,alonx,alaty,'ST,alonx,alat' ideglo = alonx aminlo = (alonx - ideglo)*60.0 iminlo = aminlo seclo = (aminlo - iminlo)*60.0 write(*,'(A,I5,F10.3,5X,2I3,F6.2)'),' I,longitude,ideg,imin,sec', I, alonx, ideglo, iminlo, seclo write(*,'(A,I5,F17.8,5X,2I3,F6.2)'),' I,JDX(I), Image UT = ', I, JDX(I), IAJDsthu, ithum, tsecu write(*,'(A,I5,F17.8,5X,2I3,F6.2)'),' I,JDX(I), Image LOCAL = ', I, JDX(I), IAJDsthl, ithum, tsecu DXCEN(J,I) = ST - ST1 DYCEN(J,I) = alaty - alaty1 XCENP(J,I) = ST*scalesta(J) YCENP(J,I) = alaty*scalesta(J) DXCENP(J,I) = XCENP(J,I) - ST1*scalesta(J) DYCENP(J,I) = YCENP(J,I) - alaty1*scalesta(J) HA(J,I) = st - RA(J) ! hour angle of the star from the central meridian in angle print*,'here 1a',HA(J,I),st,RA(J) if(HA(J,I).gt.180.0) HA(J,I) = HA(J,I) - 360.0 ! hour angle towards the east - if(HA(L,I).lt.-180.0) HA(J,I) = HA(J,I) + 360.0 ! hour angle towards the west + DZ(J,I) = (90.0 - DEC (J))*cosd(HA(J,I)) ! distance to pole from the star along the celestial meridian print*,'here 1a',HA(J,I),st,RA(J) print*,'here 1b',HA(J,I),(90.0 - DEC (J)),(90.0 - DEC (J))*cosd(HA(J,I)) DM(J,I) = (90.0 - DEC (J))*sind(HA(J,I)) ! distance to the celestial meridian from the star DZP(J,I) = DZ(J,I)*scalesta(J) ! distance to pole from the star along the celestial meridian in pixels DMP(J,I) = DM(J,I)*scalesta(J) ! distance to the celestial meridian from the star in pixels DZZ(J,I) = 90.0 - DZ(J,I) ! theoretical Y distance from this star to the north pole along the central meridian DZZP(J,I) = (DZZ(J,I) - alaty)*scalesta(J) ! theoretical Y distance from this star to the image center along the central meridian in pixels DZZPT(J,I) = sqrt(DZZP(J,I)**2 + DMP(J,I)**2) ! theoretical total distance from this star to the image center in pixels write(*,'(A,I2,I5,F10.3,4F10.6)')' here 2 ',J,I,HA(J,I),DZ(J,I),DZP(J,I),DM(J,I),DMP(J,I) C print*,'I = ',I,'DZ = ',DZ(J,I),' 2*DZ = ', 2.0*DZ(J,I) C if(I.eq.500) stop C if(I.eq.600) stop C if(I.eq.1000) stop C if(I.eq.1250) stop C if(I.eq.1500) stop C if(I.eq.2000) stop C if(I.eq.2500) stop C if(I.eq.3000) stop C if(I.eq.3500) stop C if(I.eq.4000) stop C if(I.eq.4200) stop C if(I.eq.4500) stop C print*,'here 2b', J,I,DZZ(J,I) C print*,'here 2c', J,I,DZZP(J,I),DMP(J,I),DZZPT(J,I) DC(J,I) = (90.0 - alaty) ! distance between the pole and the latitude of ASHI in degrees DCP(J,I) = (90.0 - alaty)*scalesta(J) ! theoretical distance to the center of the image from the pole in pixels DCZX(J,I) = cosd(ST)/(90.0 - alaty) ! location of X at the zenith in pixels from the pole. The image center DCZY(J,I) = sind(ST)/(90.0 - alaty) ! location of Y at the zenith in pixels from the pole. The image center DCZXP(J,I) = DCZX(J,I)*scalesta(J) DCZYP(J,I) = DCZY(J,I)*scalesta(J) C if(J.eq.1) then XCN = distcx ! = distcx: pixel center of the image in X according to Matthew YCN = distcy ! = distcy: pixel center of the image in Y according to Matthew XS1MCOP(J,I) = XS1(I) - distcx YS1MCOP(J,I) = YS1(I) - distcy XS1MCO = XS1MCOP(J,I) YS1MCO = YS1MCOP(J,I) CSOXYP = SQRT(XS1MCOP(J,I)**2 + YS1MCOP(J,I)**2) C The below intends to find the extra rotation needed to get to the Pole location from the current Rotated Pole position for each star, and C in addition we also want to find the new center location for each image as derived from the min -max positions. XSBRO(J,I) = -DMP(J,I) + distcx ! X pole star position before rotation YSBRO(J,I) = DZP(J,I) + CSOXYP + distcy ! Y pole star position relative to star and y before rotation PSXY(J,I) = SQRT((-DMP(J,I))**2 + (DZP(J,I) + CSOXYP)**2) ! pole star distance from X, Y before rotation AXY = asind(-DMP(J,I)/PSXY(J,I)) ! This extra rotation angle to the pole from the star is subject to the scale value used ! in the analysis, and so is best determined from Polaris but should work for other stars. CSOYP = PSXY(J,I) ! pole star distance from X, Y before rotation print*,'here 2d XSBRO, YSBRO, PSXY',XSBRO(J,I), YSBRO(J,I), PSXY(J,I) print*,'here 2e', J,I,AXY,CSOYP distcymydCS = (distcy - YS1(I))/CSOXYP if(distcymydCS.gt.1.0) distcymydCS = 1.0 ! if this condition is not met, the rotation will give a NAN if(distcymydCS.lt.-1.0) distcymydCS = -1.0 ! if this condition is not met, the rotation will give a NAN ROTAP(I) = asind(distcymydCS) ! Rotation counterclockwise to get pole position at top C if(I.eq.1000.and.J.eq.2) print*, XS1(I),distcx,YS1(I),distcy XCN = XCN + DXCENP(J,I) ! However the new center is at a different x location in pixels now YCN = YCN + DYCENP(J,I) ! However the new center is at a different Y location in pixels now XS1MC = XS1(I) - XCN YS1MC = YS1(I) - YCN XS1MCP(J,I) = XS1MC YS1MCP(J,I) = YS1MC CXS1 = SQRT(XS1MC*XS1MC + YS1MC*YS1MC) + DZP(J,I) ! One theoretical distance of pole to image center C All this below is only approximate if(J.eq.1) then ! OK for Polaris xOPpos(J,I) = distcx - DMP(J,I) ! ~X Distance of the Pole location derived from the star using the original X center yOPpos(J,I) = distcy - CSOYP ! ~Y Distance of the Pole location derived from the star using the original Y center CMP(J,I) = CSOYP ! ~Rotated Distance to the Pole derived from the star using the original Y center CTP(J,I) = DCP(J,I) ! ~Theoretical Rotated Dist to the Pole derived from the star using the original Y center C CMP(J,I) = CXS1 else ! Better for other stars or Saturn C CTP(J,I) = DCP(J,I) ! worked pretty well in average CTP(J,I) = DZZPT(J,I) ! should be accurate CMP(J,I) = sqrt(XS1MCOP(J,I)**2+YS1MCOP(J,I)**2) ! Best for non-pole star: observed distance - star to zenith end if print*,'Here 3a',xOPpos(J,I),HA(J,I),DZP(J,I),distcx print*,'Here 3b',yOPpos(J,I),HA(J,I),DMP(J,I),distcy print*,'Here 3c',CTP(J,I),XS1MCOP(J,I),YS1MCOP(J,I),CMP(J,I) if(XS1MC.gt.CXS1) XS1MC = CXS1 if(-XS1MC.gt.CXS1) XS1MC = -CXS1 C ROTAP(I) = asind(XS1MC/CXS1) + atand(DMP(J,I)/CXS1) ! rotation angle to place the pole in the same location print*, 'XS1(I),YS1(I),XS1MC,YS1MC',XS1(I),YS1(I),XS1MC,YS1MC print*, 'CXS1,XS1MC',CXS1,XS1MC print*, 'J,I,atand(DMP(J,I)/CXS1)',J,I,atand(DMP(J,I)/CXS1) print*,'asind(XS1MC/CXS1),ROTAP(I)',asind(XS1MC/CXS1),ROTAP(I) print*,'XCENP(J,I),DXCENP(J,I),YCENP(J,I),DYCENP(J,I)', XCENP(J,I),DXCENP(J,I),YCENP(J,I),DYCENP(J,I) print*, 'YS1MCO,XS1MCO,ROTAP(I) original',YS1MCO,XS1MCO,ROTAP(I) print*, 'YS1MCO,XS1MCO,ROTAP(I) modified',YS1MCO,XS1MCO,ROTAP(I) if(YS1MCO.gt.0.0.and.XS1MCO.gt.0.0) ROTAP(I) = 90.0 + ROTAP(I) if(YS1MCO.le.0.0.and.XS1MCO.ge.0.0) ROTAP(I) = 90.0 + ROTAP(I) if(YS1MCO.gt.0.0.and.XS1MCO.le.0.0) ROTAP(I) = 270.0 - ROTAP(I) if(YS1MCO.le.0.0.and.XS1MCO.le.0.0) ROTAP(I) = 270.0 - ROTAP(I) print*, 'YS1MCO,XS1MCO,ROTAP(I) "fixed" ',YS1MCO,XS1MCO,ROTAP(I) ROTA(J,I) = ROTAP(I) C Below for subrouine nsew_max.f that allows an Excel file produced to provide min-max values for scale determination min_max = 1 if(min_max.eq.1) then C XCENPOP is the X pole distance from the original X center of the image from the pole in pixels C YCENPOP is the Y pole distance from the original Y center of the image from the pole in pixels if(ROTA(J,I).gt.0.0.and.ROTA(J,I).lt.180.0) XCENPOP(J,I) = XS1(I) + DZP(J,I) ! west if(ROTA(J,I).lt.360.0.and.ROTA(J,I).gt.180.0) XCENPOP(J,I) = XS1(I) - DZP(J,I) ! east if(ROTA(J,I).gt.270.0.or.ROTA(J,I).lt.90.0) YCENPOP(J,I) = YS1(I) + DZP(J,I) ! north if(ROTA(J,I).gt.90.0.and.ROTA(J,I).lt.270.0) YCENPOP(J,I) = YS1(I) - DZP(J,I) ! south DXCENPOP(J,I) = XCENPO(J,I) - distcx DYCENPOP(J,I) = YCENPO(J,I) - distcy XCENP(J,I) = DC(J,I) YCENP(J,I) = DC(J,I) print*, ' ' if(ROTA(J,I).gt.0.0.and.ROTA(J,I).lt.180.0) write(*,'(A,I2,I5,F9.4,3F7.1,F8.3,F7.2)') & 'west ',J,I,HA(J,I),XCENPOP(J,I),XS1(I),DZP(J,I),XCENP(J,I),ROTA(J,I) if(ROTA(J,I).lt.360.0.and.ROTA(J,I).gt.180.0) write(*,'(A,I2,I5,F9.4,3F7.1,F8.3,F7.2)') & 'east ',J,I,HA(J,I),XCENPOP(J,I),XS1(I),DZP(J,I),XCENP(J,I),ROTA(J,I) if(ROTA(J,I).gt.270.0.or.ROTA(J,I).lt.90.0) write(*,'(A,I2,I5,F9.4,3F7.1,F8.3,F7.2)') & 'north',J,I,HA(J,I),YCENPOP(J,I),YS1(I),DZP(J,I),YCENP(J,I),ROTA(J,I) if(ROTA(J,I).gt.90.0.and.ROTA(J,I).lt.270.0) write(*,'(A,I2,I5,F9.4,3F7.1,F8.3,F7.2)') & 'south',J,I,HA(J,I),YCENPOP(J,I),YS1(I),DZP(J,I),YCENP(J,I),ROTA(J,I) C if(I.eq.4038) stop C if(I.eq.4335) stop C if(I.eq.4172) stop C if(I.eq.4438) stop C if(I.eq.625) stop C if(I.eq.833) stop C if(I.eq.713) stop C if(I.eq.1200) stop end if C End below for subroutine nsew_max.f and check tests that allow an Excel file produced to provide min-max values for scale determination C Below for subroutine find_pole_pos.f and check tests that allow a determination of the pole pixel position for a star ifind_pp = 1 if(ifind_pp.eq.1) then C The below intends to find the extra rotation needed to get to the Pole location from the current Rotated Pole position for each star, and C in addition we also want to find the new center location for each image as derived from the min -max positions. print*,'AXY',AXY DDTANP = -DMP(J,I)/tand(AXY) ! The distance along the celestial meridian out to the perpendicular to the meridian through x, y print*,'Compare DDTANP, CSOYP ~SAME?', DDTANP, CSOYP DCPMDDTANP = DCP(J,I) - DDTANP ! Distance in pixels different between the latitude location and DDTANP XRO = atand(DCPMDDTANP/-DMP(J,I)) + (90.0 - AXY) ! The extra angle that the diameter needs to rotate through greater than ROTA and AXY DIAM = SQRT(DCPMDDTANP**2 + DMP(J,I)**2) ! The diameter beyond the x,y position TRXO = - XRO - ROTA(J,I) - AXY ! The total rotation angle to get to the end ofdiameter TRO = TRXO + 90.0 ! just the increment to the perpendicular to the box with the diameter print*,' ' print*,'TRXO, TRO',TRXO, TRO if(TRO.lt.-360.0) TRO = TRO + 360.0 C From X, Y center the lat Lon center is at: AROTAXC = sind(90.0-TRO)*DIAM ! The X position of the image latitude, longitude center AROTAYC = cosd(90.0-TRO)*DIAM print*,' ' print*,'DDTANP, DCPMDDTANP, XRO',DDTANP, DCPMDDTANP, XRO print*,'DMP, DIAM',DMP(J,I),DIAM print*,'TRXO, TRO',TRXO, TRO print*,'AROTAXC, AROTAYC', AROTAXC, AROTAYC print*,' ' print*,'CSOXYP DZP',CSOXYP,DZP(J,I) print*,'DMP DCP CMP',DMP(J,I), DCP(J,I), CMP(J,I) AXYP = SQRT(DMP(J,I)**2 + DCPMCMP**2) ! Distance in pixels for the diameter to get to the new Pole image center on the rotated image AXYR = atand(-DMP(J,I)/DCPMCMP) AXYRO = 180.0 + atand(-DMP(J,I)/DCPMCMP) ! Approximate angle from the X,Y position to the Apex location DOCPX = AROTAXC ! the difference from the X center lat-lon position to the X center on the original I image DOCPY = AROTAYC ! the difference from the Y center lat-lon position to the Y center on the original I image DOPPLX = -sind(-ROTA(J,I) - AXY)*CMP(J,I) ! It seems this value is always negative DOPPLY = cosd(-ROTA(J,I) - AXY)*CMP(J,I) ! It seems this value is always positive AROTA(J,I) = -ROTA(J,I) - AXY if(AROTA(J,I).lt.-360.0) AROTA(J,I) = AROTA(J,I) + 360.0 XCENPOP(J,I) = XS1(I) + DZP(J,I) OCPX(J,I) = DOCPX + distcx ! the X center lat-lon location position on the original I image OCPY(J,I) = DOCPY + distcy ! the Y center lat-lon location position on the original I image OPPLX(J,I) = DOPPLX + distcx ! the X pole location position on the original I image OPPLY(J,I) = DOPPLY + distcy ! the Y pole location position on the original I image print*,' ' write(*,'(A,I5,A,F7.2,A,F7.2)')'Image =',I,', Rotation (degrees) = ',ROTA(J,I) ,', HA(J,I) (degrees) = ',HA(J,I) write(*,'(A,I5,A,F7.2,A,F7.2)')'Image =',I,', Pole Ro frm N (deg) = ',AXY ,',Min Rot+AXY(AROTA)deg = ',AROTA(J,I) write(*,'(A,I5,A,F7.2,A,F7.2)')'Image =',I,', Extra Rot frm Pole = ',AXYR ,', Rot of Diam frm Pole = ',AXYRO write(*,'(A,I5,A,F7.2,A,F7.2)')'Image =',I,', distcx = ',distcx ,', distcy = ',distcy write(*,'(A,I5,A,F7.2,A,F7.2)')'Image =',I,', Or Pol diam frm X,Y = ',AXYP ,', Minus Rot-AROTAC(deg)= ',AROTAC write(*,'(A,I5,A,F7.2,A,F7.2)')'Image =',I,', Pole XY dif frm Orig = ',AXYP write(*,'(A,I5,A,F7.2,A,F7.2)')'Image =',I,', Polaris frm X,Y = ',CSOXYP ,', Extra Dist to Pole = ',DZP(J,I) write(*,'(A,I5,A,F7.2,A,F7.2)')'Image =',I,', Dist to pole frm X,Y = ',CSOYP ,', Pole dis to la,lo=DCP= ',DCP(J,I) print*,' Values at AROTA, AROTAC' write(*,'(A,I5,A,F7.2,A,F7.2)')'Image =',I,', Lon C dif from X Cen = ',DOCPX ,', Lat C dif from Y Cen = ',DOCPY write(*,'(A,I5,A,F7.2,A,F7.2)')'Image =',I,', X la,lo position = ',OCPX(J,I) ,', Y la,lo position = ',OCPY(J,I) write(*,'(A,I5,A,F7.2,A,F7.2)')'Image =',I,', Pol X dif from X = ',DOPPLX ,', Pol Y dif from Y = ',DOPPLY write(*,'(A,I5,A,F7.2,A,F7.2)')'Image =',I,', X Pol from X pos = ',OPPLX(J,I),', Y Pole from Y pos = ',OPPLY(J,I) C if(I.eq.500) stop ! Matthew's image C if(I.eq.625) stop ! ROTA = 0.57 AROTA = -1.35 DIAM = 12.7632 C if(I.eq.833) stop ! ROTA = 179.66 AROTA = -180.42 DIAM = 7.7545 C if(I.eq.4175) stop ! ROTA = 358.12 AROTA = -358.49 C if(I.eq.4170) stop ! ROTA = 1.89 AROTA = -2.13 DIAM = 10.6938 C if(I.eq.4042) stop ! ROTA = 88.35 AROTA = -88.62 C if(I.eq.4036) stop ! ROTA = 91.90 AROTA = -92.27 C if(I.eq.4438) stop ! ROTA = 179.54 AROTA = -179.70 C if(I.eq.4430) stop ! ROTA = 184 C if(I.eq.4335) stop ! ROTA = 271.08 AROTA = -271.28 C if(I.eq.2826) stop ! ROTA = 1.27 AROTA = -1.83 DIAM = 10.9487 C if(I.eq.3526) stop ! ROTA = 180.26 AROTA = -180.66 DIAM = 4.1243 ST2 = ST if(ST2.gt.360.0) ST2 = ST2 - 360.0 ! Sidereal time that does not go over 360 degrees C Here begin the image search for each star. Although each Star can be searched, The RA, and Decs are only currently determined by Polaris print*,' ' print*,' ' print*,'Now begin the determination of RA and Dec on an image and the image search ' print*,' for a few known stellar positions at their designated RA and DEC' DISPXYC(J,I) = SQRT((OPPLX(J,I)-distcx)**2 + (OPPLY(J,I)-distcy)**2) print*,'DISDEC,DCP(J,I)',DCP(J,I),DISPXYC(J,I) ! DISDEC and DCP should be the same, and is. DISPXYC(J,I) determined here C if(I.eq.4170) stop C ipij = 0 C ipij = 1 C IMC = 0 ! 0 if BVJ's center is used IMC = 1 ! 1 if Matthew's center is used lsource1 = 11 ! first source attempted lsource2 = 23 ! last source attempted C lsource1 = 19 ! first source attempted C lsource2 = 19 ! last source attempted if(I.eq.5000) then do JJJ=lsource1,lsource2 NII(JJJ) = 0 NJJ(JJJ) = 0 if(JJJ.eq.1) then call read_image_t_x_y(cfile1,NOBSM,NOBS,JDXSD,XS11,YS11) ! read an ashi file that gives time and x, y of this star write(*,'(A,I5,A,F9.3,A,F9.3)')'This star image',I,' is at XS11 =',XS11(I),' and a YS11 =',YS11(I) Write(*,'(A,I5,A,F9.3,A,F9.3)')'Polaris is at: ',I,' has an RA =',RA(JJJ),' and a Dec =',DEC(JJJ) end if if(JJJ.eq.2) then call read_image_t_x_y(cfile2,NOBSM,NOBS,JDXSD,XS12,YS12) ! read an ashi file that gives time and x, y of this star write(*,'(A,I5,A,F9.3,A,F9.3)')'This star image',I,' is at XS12 =',XS12(I),' and a YS12 =',YS12(I) Write(*,'(A,I5,A,F9.3,A,F9.3)')'This star image',I,' has an RA =',RA(JJJ),' and a Dec =',DEC(JJJ) end if if(JJJ.eq.3) then call read_image_t_x_y(cfile3,NOBSM,NOBS,JDXSD,XS13,YS13) ! read an ashi file that gives time and x, y of this star write(*,'(A,I5,A,F9.3,A,F9.3)')'This star image',I,' is at XS13 =',XS13(I),' and a YS13 =',YS13(I) Write(*,'(A,I5,A,F9.3,A,F9.3)')'This star image',I,' has an RA =',RA(JJJ),' and a Dec =',DEC(JJJ) end if if(JJJ.eq.4) then call read_image_t_x_y(cfile4,NOBSM,NOBS,JDXSD,XS14,YS14) ! read an ashi file that gives time and x, y of this star write(*,'(A,I5,A,F9.3,A,F9.3)')'This star image',I,' is at XS14 =',XS14(I),' and a YS14 =',YS14(I) Write(*,'(A,I5,A,F9.3,A,F9.3)')'This star image',I,' has an RA =',RA(JJJ),' and a Dec =',DEC(JJJ) end if if(JJJ.eq.5) then call read_image_t_x_y(cfile5,NOBSM,NOBS,JDXSD,XS15,YS15) ! read an ashi file that gives time and x, y of this star write(*,'(A,I5,A,F9.3,A,F9.3)')'This star image',I,' is at XS15 =',XS15(I),' and a YS15 =',YS15(I) Write(*,'(A,I5,A,F9.3,A,F9.3)')'Saturn is at: ',I,' has an RA =',RA(JJJ),' and a Dec =',DEC(JJJ) end if if(JJJ.eq.6) then irotimage = 1 ! 1, a rotated image is used Write(*,'(A,I5,A,F9.3,A,F9.3)')'Capella is at: ',I,' has an RA =',RA(JJJ),' and a Dec =',DEC(JJJ) end if if(JJJ.eq.7) then irotimage = 1 ! 1, a rotated image is used Write(*,'(A,I5,A,F9.3,A,F9.3)')'B Cetus is at: ',I,' has an RA =',RA(JJJ),' and a Dec =',DEC(JJJ) end if if(JJJ.ge.11.or.JJJ.le.23) then irotimage = 1 ! 1, a rotated image is used Write(*,'(2(A,I5),A,F9.3,A,F9.3)') & 'Rotated Star:',JJJ,' Image',I,' has an RA =',RA(JJJ),' and a Dec =',DEC(JJJ) end if print*,'I, IMC = ',I,IMC print*,'I,lsource1,lsource2,irotimage',I,lsource1,lsource2,irotimage do II=1,2000 ! X pixels do JJ=1,2000 ! Y pixels XII(II) = float(II) ! Float X value of original image pixel YJJ(JJ) = float(JJ) ! Float Y value of original image pixel C if(II.eq.1028.and.JJ.eq.991) then ! center RA, 100 pixels east Dec - OK C if(II.eq.1108.and.JJ.eq.991) then ! center RA, 20 pixels east Dec - OK C if(II.eq.1128.and.JJ.eq.997) then ! center RA, Dec - up ~ 5 deg pixels OK C if(II.eq.1123.and.JJ.eq.991) then if(JJJ.eq.1) then if(II.ge.1143.and.II.le.1146.and.JJ.ge.1558.and.JJ.le.1560) then ! RA, Dec position for Polaris on image 4170 ipij = 1 write(*,'(A,2I5,2F8.3)')'Got here for Polaris at II, JJ, RA, DEC',II,JJ,RA(JJJ),DEC(JJJ) else ipij = 0 end if end if if(JJJ.eq.2) then if(II.ge.1604.and.II.le.1608.and.JJ.ge.526.and.JJ.le.530) then ! RA, Dec position for B Capricorn image 4170 ipij = 1 write(*,'(A,2I5,2F8.3)')'Got here for Beta Capricorn at II, JJ, RA, DEC',II,JJ,RA(JJJ),DEC(JJJ) else ipij = 0 end if end if if(JJJ.eq.6) then if(RAIJ(II-1,JJ-1).ge.(RA(JJJ)-1.0).and.RAIJ(II-1,JJ-1).le.(RA(JJJ)+1.0).and. ! RA, Dec position for Capella, image 4170 & DECIJ(II-1,JJ-1).ge.(DEC(JJJ)-1.0).and.DECIJ(II-1,JJ-1).le.(DEC(JJJ)+1.0)) then ipij = 1 ! Switch on printing write(*,'(A,2I5,2F8.3)')'Got here for Capella at II, JJ, RA, DEC',II,JJ,RA(JJJ),DEC(JJJ) else ipij = 0 end if end if if(JJJ.eq.7) then if(RAIJ(II-1,JJ).ge.(RA(JJJ)-1.0).and.RAIJ(II-1,JJ).le.(RA(JJJ)+1.0).and. ! RA, Dec position for Beta Cetus, image 4170 & DECIJ(II-1,JJ).ge.(DEC(JJJ)-1.0).and.DECIJ(II-1,JJ).le.(DEC(JJJ)+1.0)) then ipij = 1 write(*,'(A,2I5,2F8.3)')'Got here for Beta Cetus at II, JJ, RA, DEC',II,JJ,RA(JJJ),DEC(JJJ) else ipij = 0 end if end if if(JJJ.ge.11) then if(RAIJ(II-1,JJ).ge.(RA(JJJ)-1.0).and.RAIJ(II-1,JJ).le.(RA(JJJ)+1.0).and. ! RA, Dec position for stars greater than 11 & DECIJ(II-1,JJ).ge.(DEC(JJJ)-1.0).and.DECIJ(II-1,JJ).le.(DEC(JJJ)+1.0)) then ipij = 1 write(*,'(A,I2,A,2I5,2F8.3)')'Got here for Star ',JJJ,' at II, JJ, RA, DEC',II,JJ,RA(JJJ),DEC(JJJ) else ipij = 0 end if end if C if(II.eq.1133.and.JJ.eq.991) then C if(II.eq.1428.and.JJ.eq.500) then ! center RA, Dec - down ~.5 deg pixels OK C if(II.eq.1148.and.JJ.eq.991) then ! center RA, 20 pixels west Dec - OK C if(II.eq.1228.and.JJ.eq.991) then ! center RA, 100 pixels west Dec - OK if(irotimage.ne.1) then ! Non-rotated stars DISCIJ = SQRT((OCPX(J,I)-XII(II))**2 + (OCPY(J,I)-YJJ(JJ))**2) ! Dist XII, YJJ to lat-lon center DISDEC = SQRT((OPPLX(J,I)-OCPX(J,I))**2 + (OPPLY(J,I)-OCPY(J,I))**2) ! Dist of Pole to lat-lon center DISIJ = SQRT((OPPLX(J,I)-XII(II))**2 + (OPPLY(J,I)-YJJ(JJ))**2) ! Dist XII, YJJ to the Pole OLXDMC = (OCPX(J,I)-XII(II))/DISCIJ if(OLXDMC.gt.1.0) OLXDMC = 1.0 if(ipij.eq.1) then print*,' ' print *,'OLXDMC', OLXDMC end if TOPANG = asind(-OLXDMC) ASIGN = -1.0 if(ipij.eq.1) then print *,'TOPANG', TOPANG print*,'XII(II), OCPX(J,I)',XII(II), OCPX(J,I) print*,'YJJ(II), OCPY(J,I)',YJJ(JJ), OCPY(J,I) end if if(YJJ(JJ).ge.OCPY(J,I)) then TANG = 1.0 ! West Stars if(TOPANG.ge.AROTA(J,I)) TANG = -1.0 ! East stars if(ipij.eq.1) then if(TANG.eq.-1.0) write(*,'(A,4F7.1)')'upper west TOPANG,AROTA,ASIGN,TANG',TOPANG,AROTA(J,I),ASIGN,TANG if(TANG.eq.1.0) write(*,'(A,4F7.1)')'upper east TOPANG,AROTA,ASIGN,TANG',TOPANG,AROTA(J,I),ASIGN,TANG end if else TANG = 1.0 ! East Stars if(TOPANG.ge.AROTA(J,I)) TANG = -1.0 ! West Stars if(ipij.eq.1) then if(TANG.eq.-1.0) write(*,'(A,4F7.1)')'lower east TOPANG,AROTA,ASIGN,TANG',TOPANG,AROTA(J,I),ASIGN,TANG if(TANG.eq.1.0) write(*,'(A,4F7.1)')'lower west TOPANG,AROTA,ASIGN,TANG',TOPANG,AROTA(J,I),ASIGN,TANG end if end if else ! Rotated stars ************ if(IMC.ne.1) then ! The distance to BVJ center used to the pole DISCIJ = SQRT((1134.4-XII(II))**2 + (990.8-YJJ(JJ))**2) ! Dist XII, YJJ to lat-lon center (image 500) C DISCIJ = SQRT((OCPX(J,I)-XII(II))**2 + (OCPY(J,I)-YJJ(JJ))**2) ! Dist XII, YJJ to lat-lon center DISDEC = 1570.5 - 990.8 ! Dist of Pole to lat-lon center (image 500) DISIJ = SQRT((1134.4-XII(II))**2 + (1570.5-YJJ(JJ))**2) ! Dist XII, YJJ to the Pole (image 500) C DISIJ = SQRT((1134.4-XII(II))**2 + (990.8+DCP(J,I)-YJJ(JJ))**2) ! Dist XII, YJJ to the Pole C DISIJ = SQRT((OCPX(J,I)-XII(II))**2 + (OCPY(J,I)+DCP(J,I)-YJJ(JJ))**2) ! Dist XII, YJJ to the Pole DCENTRO = SQRT((1134.4-XS1rot(JJJ,I))**2 + (990.8-YS1rot(JJJ,I))**2) ! distance from BVJ's center to the star centroid TANG = 1.0 if(XII(II).lt.1134.4) then ! If the distance to BVJ's center is used to the pole ASIGN = -1.0 if(ipij.eq.1) write(*,'(A,4F7.1)')'east AROTA,ASIGN',AROTA(J,I),ASIGN else ASIGN = 1.0 if(ipij.eq.1) write(*,'(A,4F7.1)')'west AROTA,ASIGN',AROTA(J,I),ASIGN end if end if if(IMC.eq.1) then ! The distance to Mathew's center is used to the pole DISCIJ = SQRT((distcx-XII(II))**2 + (distcy-YJJ(JJ))**2) ! Dist XII, YJJ to Matthew's X.Y center DISDEC = SQRT((1134.4-distcx)**2 + (distcy-1570.5)**2) ! Distance of Matthew's X,Y center to the pole DISIJ = SQRT(((distcx+(distcx-1134.4))-XII(II))**2 + (1570.5-YJJ(JJ))**2) ! Dist XII, YJJ to the Pole DCENTRO = SQRT(((distcx+(distcx-1134.4))-XS1rot(JJJ,I))**2 + (distcy-YS1rot(JJJ,I))**2) ! distance from Matthew's centroid to cent TANG = 1.0 if(XII(II).lt.distcx) then ASIGN = -1.0 if(ipij.eq.1) write(*,'(A,4F7.1)')'east AROTA,ASIGN',AROTA(J,I),ASIGN else ASIGN = 1.0 if(ipij.eq.1) write(*,'(A,4F7.1)')'west AROTA,ASIGN',AROTA(J,I),ASIGN end if end if end if ANUM = DISIJ**2 + DISDEC**2 - DISCIJ**2 HAIJ(II,JJ) = acosd(ANUM/(2.0*DISIJ*DISDEC)) DISIJR (II,JJ) = DISIJ ! Dist Pole to XII, YJJ DISDECR(II,JJ) = DISDEC ! Dist Pole to the center DISCIJR(II,JJ) = DISCIJ ! Dist center to XII, YJJ DCENTROR(JJJ,I) = DCENTRO ! Distance from the center to the star centroid if(ipij.eq.1) then print*,'XII(II),YJJ(JJ) =',XII(II),YJJ(JJ) print*,'Distances =',DISIJ,DISDEC,DISCIJ print*,'HAIJ(II,JJ) =',ASIGN*HAIJ(II,JJ) ! Hour angle of star print*,'ST2 =',ST2 print*,'Ingrediants =',TANG,ST2,-ASIGN*HAIJ(II,JJ),-AROTA(J,I) if(irotimage.ne.1) print*,'RAIJ(II,JJ) =',TANG*(ST2 - ASIGN*HAIJ(II,JJ) - AROTA(J,I)), irotimage if(irotimage.eq.1) print*,'RAIJ(II,JJ),irot=',TANG*(ST2 - ASIGN*HAIJ(II,JJ)), irotimage print*,'DECIJ(II,JJ) =', 90.0 - DISIJ/scalesta(J) end if if(irotimage.ne.1) RAIJ (II,JJ) = TANG*(ST2 - ASIGN*HAIJ(II,JJ) + AROTA(J,I)) ! RA of XII, YJJ if(irotimage.eq.1) RAIJ (II,JJ) = TANG*(ST2 - ASIGN*HAIJ(II,JJ)) ! RA of XII, YJJ if(RAIJ (II,JJ).gt.360.0) RAIJ(II,JJ) = RAIJ(II,JJ) - 360.0 if(RAIJ (II,JJ).lt. 0.0) RAIJ(II,JJ) = RAIJ(II,JJ) + 360.0 DECIJ(II,JJ) = 90.0 - DISIJ/scalesta(J) ! Dec of XII, YJJ if(ipij.eq.1) then print*,' ' write(*,'(A,I5,A,F7.2,A,F7.2)')'Image number = ',I,', Rotation =', ROTA(J,I),', HA(J,I) = ',HA(J,I) write(*,'(A,I5,A,F7.2,A,F7.2)')'Image number = ',I,', AXY =', AXY if(JJJ.eq.1) then write(*,'(A,I5,A,F7.2,A,F7.2)')'Image number = ',I,', XS11 =', XS11(I),', YS11 = ',YS11(I) end if if(JJJ.eq.2) then write(*,'(A,I5,A,F7.2,A,F7.2)')'Image number = ',I,', XS12 =', XS12(I),', YS12 = ',YS12(I) end if if(JJJ.eq.3) then write(*,'(A,I5,A,F7.2,A,F7.2)')'Image number = ',I,', XS13 =', XS13(I),', YS13 = ',YS13(I) end if if(JJJ.eq.4) then write(*,'(A,I5,A,F7.2,A,F7.2)')'Image number = ',I,', XS14 =', XS14(I),', YS14 = ',YS14(I) end if if(JJJ.eq.5) then write(*,'(A,I5,A,F7.2,A,F7.2)')'Image number = ',I,', XS15 =', XS15(I),', YS15 = ',YS15(I) end if write(*,'(A,F7.1,A,F7.1)') 'Original OCPX(J,I) =',OCPX(J,I),', OCPY(J,I) =',OCPY(J,I) write(*,'(A,F7.1,A,F7.1,A,F7.1)')'Original OPPLX(J,I) =',OPPLX(J,I),', OPPLY(J,I) =',OPPLY(J,I) write(*,'(A,F7.1,A,F7.1,A,F9.3)')'XII =',XII(II),', YJJ =',YJJ(JJ), ', ST =',ST2 write(*,'(A,F7.1,A,F7.1,A,F7.1)')'DISDEC =',DISDEC,', DISCIJ =',DISCIJ,', DISIJ =',DISIJ write(*,'(A,F7.3,A,F7.3)') 'HAIJ(II,JJ) = ',ASIGN*HAIJ(II,JJ) write(*,'(A,F7.3,A,F7.3)') 'RAIJ(II,JJ) = ',RAIJ(II,JJ) ,', DECIJ(II,JJ)= ',DECIJ(II,JJ) end if end do end do end do C stop do JJJ=lsource1,lsource2 do II=1,2000 ! X pixels do JJ=1,2000 ! Y pixels C FAC = 90.0/(90.0 - DEC(JJJ)) C RAC = FAC*SQRT(FAC) C DAC = 0.5*SQRT(RAC) RAC = 0.5 DAC = 0.5 if(RAIJ(II,JJ).gt.(RA(JJJ)-1.0*RAC).and.RAIJ(II,JJ).lt.(RA(JJJ)+1.0*RAC)) then if(DECIJ(II,JJ).lt.0.0) then C print*,'to here 7', RAIJ(II,JJ),RA(JJJ),DECIJ(II,JJ),DEC(JJJ) C print*,'to here 7b',XII(II),YJJ(JJ) C print*,'to here 7c',HAIJ(II,JJ),ST,ST2 end if if(DECIJ(II,JJ).gt.(DEC(JJJ)-1.0*DAC).and.DECIJ(II,JJ).lt.(DEC(JJJ)+1.0*DAC)) then print*,'to here 8',RAIJ(II,JJ),RA(JJJ),DECIJ(II,JJ),DEC(JJJ) NII(JJJ) = NII(JJJ) + 1 NJJ(JJJ) = NJJ(JJJ) + 1 if(NII(JJJ).gt.500) go to 1111 print*,'to here 8b',NII(JJJ),NJJ(JJJ) XXII (JJJ,NII(JJJ),NJJ(JJJ)) = XII(II) YYII (JJJ,NII(JJJ),NJJ(JJJ)) = YJJ(JJ) RAIJS (JJJ,NII(JJJ),NJJ(JJJ)) = RAIJ(II,JJ) DECIJS (JJJ,NII(JJJ),NJJ(JJJ)) = DECIJ(II,JJ) DISIJRS (JJJ,NII(JJJ),NJJ(JJJ)) = DISIJR(II,JJ) ! Dist Pole to XII, YJJ DISDECRS(JJJ,NII(JJJ),NJJ(JJJ)) = DISDECR(II,JJ) DISCIJRS(JJJ,NII(JJJ),NJJ(JJJ)) = DISCIJR(II,JJ) ! Dist center to XII, YJJ DCENTRORS(JJJ,NII(JJJ),NJJ(JJJ)) = DCENTROR(JJJ,I) print*,'to here 9a',DCENTRORS(JJJ,NII(JJJ),NJJ(JJJ)) write(*,'(A,I5,4F7.1,2F9.3)')'to here 9', & JJJ,XII(II),YJJ(JJ),XXII(JJJ,NII(JJJ),NJJ(JJJ)),YYII(JJJ,NII(JJJ),NJJ(JJJ)), & RAIJS (JJJ,NII(JJJ),NJJ(JJJ)),DECIJS(JJJ,NII(JJJ),NJJ(JJJ)) 1111 continue end if end if end do C stop end do write(*,'(A,2I4)')'NII(JJJ),NJJ(JJJ)', NII(JJJ),NJJ(JJJ) if(NII(JJJ).gt.500) then NII(JJJ) = 500 NJJ(JJJ) = 500 end if end do do IXYS=1,NOBS XS16(IXYS) = 0.0 YS16(IXYS) = 0.0 end do do JJJ=lsource1,lsource2 print*,' ' print*,' ' if(irotimage.ne.1) print*, ' Star Observed R,D Calculated Actual Determined' if(irotimage.eq.1) print*, ' Star # X,Y Location R,D Actual R,D Determine Dist P Dist C' if(irotimage.ne.1) then ! Non-rotated stars do II=1,NII(JJJ) JJ=II if(JJJ.eq.1) write(*,'(3I5,4F7.1,4F8.3)') JJJ,II,JJ,XS11(I),YS11(I),XXII(JJJ,II,JJ),YYII(JJJ,II,JJ), & RA(JJJ),DEC(JJJ),RAIJS(JJJ,II,JJ),DECIJS(JJJ,II,JJ) if(JJJ.eq.2) write(*,'(3I5,4F7.1,4F8.3)') JJJ,II,JJ,XS12(I),YS12(I),XXII(JJJ,II,JJ),YYII(JJJ,II,JJ), & RA(JJJ),DEC(JJJ),RAIJS(JJJ,II,JJ),DECIJS(JJJ,II,JJ) if(JJJ.eq.3) write(*,'(3I5,4F7.1,4F8.3)') JJJ,II,JJ,XS13(I),YS13(I),XXII(JJJ,II,JJ),YYII(JJJ,II,JJ), & RA(JJJ),DEC(JJJ),RAIJS(JJJ,II,JJ),DECIJS(JJJ,II,JJ) if(JJJ.eq.4) write(*,'(3I5,4F7.1,4F8.3)') JJJ,II,JJ,XS14(I),YS14(I),XXII(JJJ,II,JJ),YYII(JJJ,II,JJ), & RA(JJJ),DEC(JJJ),RAIJS(JJJ,II,JJ),DECIJS(JJJ,II,JJ) if(JJJ.eq.5) write(*,'(3I5,4F7.1,4F8.3)') JJJ,II,JJ,XS15(I),YS15(I),XXII(JJJ,II,JJ),YYII(JJJ,II,JJ), & RA(JJJ),DEC(JJJ),RAIJS(JJJ,II,JJ),DECIJS(JJJ,II,JJ) if(JJJ.gt.5) write(*,'(3I5,4F7.1,4F8.3)') JJJ,II,JJ,XS16(I),YS16(I),XXII(JJJ,II,JJ),YYII(JJJ,II,JJ), & RA(JJJ),DEC(JJJ),RAIJS(JJJ,II,JJ),DECIJS(JJJ,II,JJ) end do else ! Rotated Stars ANII = 0.0 XXIIA = 0.0 YYIIA = 0.0 RAIJSA = 0.0 DECIJSA = 0.0 DISIJRSA = 0.0 DISCIJRSA = 0.0 do II=1,NII(JJJ) JJ = II if(IMC.ne.1) write(*,'(A,I3,I4,2F7.1,2(F8.3,F7.3),2F7.1)') 'BVJs Center ',JJJ,II,XXII(JJJ,II,JJ), & YYII(JJJ,II,JJ),RA(JJJ),DEC(JJJ),RAIJS(JJJ,II,JJ),DECIJS(JJJ,II,JJ),DISIJRS(JJJ,II,JJ),DISCIJRS(JJJ,II,JJ) if(IMC.eq.1) write(*,'(A,I3,I4,2F7.1,2(F8.3,F7.3),2F7.1)') 'Matts Center',JJJ,II,XXII(JJJ,II,JJ), & YYII(JJJ,II,JJ),RA(JJJ),DEC(JJJ),RAIJS(JJJ,II,JJ),DECIJS(JJJ,II,JJ),DISIJRS(JJJ,II,JJ),DISCIJRS(JJJ,II,JJ) ANII = ANII + 1.0 XXIIA = XXIIA + XXII (JJJ,II,JJ) YYIIA = YYIIA + YYII (JJJ,II,JJ) RAIJSA = RAIJSA + RAIJS (JJJ,II,JJ) DECIJSA = DECIJSA + DECIJS (JJJ,II,JJ) DISIJRSA = DISIJRSA + DISIJRS (JJJ,II,JJ) DISCIJRSA = DISCIJRSA + DISCIJRS(JJJ,II,JJ) DCENTRORSA= DCENTRORS (JJJ,II,JJ) end do print*,' ' if(irotimage.eq.1) print*, ' Star # X,Y Location R,D Actual R,D Determine Dist P Dist C' if(ANII.ne.0.0) then if(IMC.ne.1) write(*,'(A,I3,I4,2F7.1,2(F8.3,F7.3),2F7.1)') 'BVJs Cen Ave',JJJ,II-1,XXIIA/ANII, & YYIIA/ANII,RA(JJJ),DEC(JJJ),RAIJSA/ANII,DECIJSA/ANII,DISIJRSA/ANII,DISCIJRSA/ANII if(IMC.eq.1) write(*,'(A,I3,I4,2F7.1,2(F8.3,F7.3),2F7.1)') 'Matt Cen Ave',JJJ,II-1,XXIIA/ANII, & YYIIA/ANII,RA(JJJ),DEC(JJJ),RAIJSA/ANII,DECIJSA/ANII,DISIJRSA/ANII,DISCIJRSA/ANII DISCIJRSAS (JJJ) = DISCIJRSA/ANII DCENTRORSAS(JJJ) = DCENTRORSA end if end if end do icentroids = 1 if(icentroids.eq.1) then write(*,'(A)') ' Star DIST C DIST Centroid DIFF DIS ' do JJJ=lsource1,lsource2 if(IMC.ne.1) write(*,'(A,I5,2F12.1,F13.5)') 'BVJs Cen DIF',JJJ,DISCIJRSAS(JJJ),DCENTRORSAS(JJJ), & DISCIJRSAS(JJJ)/DCENTRORSAS(JJJ) if(IMC.eq.1) write(*,'(A,I5,2F12.1,F13.5)') 'Mats Cen DIF',JJJ,DISCIJRSAS(JJJ),DCENTRORSAS(JJJ), & DISCIJRSAS(JJJ)/DCENTRORSAS(JJJ) end do end if C stop end if C if(I.eq.500) stop C if(I.eq.625) stop ! ROTA = 0 C if(I.eq.1200) stop ! ROTA = 90 C if(I.eq.833) stop ! ROTA = 180 C if(I.eq.713) stop ! ROTA = 270 C if(I.eq.4172) stop ! ROTA = .4 C if(I.eq.4175) stop ! ROTA = 358 C if(I.eq.4038) stop ! ROTA = 90.8 C if(I.eq.4042) stop ! ROTA = 88.3 C if(I.eq.4438) stop ! ROTA = 179 C if(I.eq.4430) stop ! ROTA = 184 C if(I.eq.4335) stop ! ROTA = 271.08 STN = ST if(ST.lt.0.0) STN = STN + 360.0 if(ST.gt.360.0) STN = STN - 360.0 ST3(I) = STN end if end do end do C stop ip_min_max = 0 ip_lat_lon = 0 ip_star_pole = 0 ip_star_pole_st = 0 ip_star_pole_st_dist = 0 ip_t_lat_lon_2rot_dt = 0 ip_t_lat_lon_2r_dtrm = 0 ip_star_pole_xy_rm = 1 Iave = 3 do J=1,1 ! Polaris Only C do J=1,NSTAR C call print_cent_rot(J,NSTAR,NOBS,xOPpos,distcx,yOPpos,distcy,XCENP,DXCENP,YCENP,DYCENP,ROTA) C call print_cent_rot(J,NSTAR,NOBS,XCENPOP,distcx,YCENPOP,distcy,XCENP,DXCENP,YCENP,DYCENP,ROTA) C call print_distcdiff(J,NSTAR,NOBS,NXYM,NVAL,distcx,distcy,XYDISTZ,XSMAX,NIMAX,XSMIN,NIMIX,YSMAX,NIMAY, C & YSMIN,NIMIY,distxcmax,distxcmin,distycmax,distycmin,ROTA) if(ip_min_max.eq.1) call nsew_max(J,NSTAR,NOBS,XCENPOP,YCENPOP,XCENP,YCENP,ROTA) if(ip_lat_lon.eq.1) call print_lat_lon(NOBS,JDX,alatyy,alonxx) C call print_max_determine(J,NSTAR,NOBS,XS1,YS1,CMP,CTP,XS1MCOP,YS1MCOP,XS1MCP,YS1MCP,ROTA) if(ip_star_pole.eq.1) call print_star_pole_cen_xy_rot(J,NSTAR,NOBS,OPPLX,OPPLY,OCPX,OCPY,XS1,YS1,ROTA,AROTA) if(ip_star_pole_st.eq.1) call print_star_pole_cen_xy_st_rot(J,NSTAR,NOBS,OPPLX,OPPLY,OCPX,OCPY,AROTA,ST3) if(ip_star_pole_st_dist.eq.1)call pr_star_pole_cen_xy_dist_st_rot(J,NSTAR,NOBS,OPPLX,OPPLY,DCP,DISPXYC,alatyy,AROTA,ST3) if(ip_t_lat_lon_2rot_dt.eq.1) call pr_t_lat_lon_2rot_delt(NOBS,NSTAR,J,JDbeg,JDX,ST3,alatyy,alonxx,AROTA) if(ip_t_lat_lon_2r_dtrm.eq.1) call pr_t_lat_lon_2rot_delt_rm(NOBS,NSTAR,J,Iave,JDbeg,JDX,ST3,alatyy,alonxx,AROTA) if(ip_star_pole_xy_rm.eq.1) call pr_star_pole_cen_xy_rm(J,NSTAR,NOBS,Iave,OPPLX,OPPLY,OCPX,OCPY,XS1,YS1) end do stop end