C C PURPOSE: C This will return the location of the Parker Solar Probe, C BepiColombo space crafts, and Solar Orbiter on a C given year and (fractional) day of year. C C INPUTS: C itp = integer, selects which space craft to iterpolate. C ITP = 17 => PSP C ITP = 18 => BC C ITP = 19 => SO C ITP = 20 => EARTH C iYr = integer, year C Doy = double precision, day of year C npoints = integer, number of points read from trajectory.txt C NMAX = integer, size(dates) C dates = real, dim(2,NMAX), dates for the location of PSP, C dates(1,:) = year, and C dates(2,:) = day of year C coords = real, dim(3,NMAX), coordinates of PSP, C coords(1,:) = radius (AU), C coords(2,:) = SE_Latitude, C coords(3,:) = SE_Longitude C C OUTPUTS: C RVect = real, dim(3), ecliptic coordinate, where C RVect(1) = radius (AU), C RVect(2) = latitude, C RVect(3) = longitude C C DEPENDENCIES: C poly_interpolate.f, C jd_from_doy.f C C NOTES: C Note the input date must be in modified julian form. C That is year, day of year. The location is determined C by 4 point interpolation of the data in the file C '/home/soft/dat/parker/psp_trajectory.txt'. So, if this file C does not exist in the above location it can be downloaded from C 'https://omniweb.gsfc.nasa.gov/coho/helios/heli.html'. C C The coordinate at index icur of the array coords(:,icur) C correspond to the location of the PSP at the date given in the C array dates(:,icur). C C VARIABLES: C MAX_ORDER = integer, parameter, maximum order expexcted C for interpolation. C order = integer, the order of interpolation c (order <= MAX_ORDER) C jd_start = integer, the julian day of the start of the C PSP data (year = 2018, doy = 225) C month = integer, number of month of year C dom = integer, number of day of the month C ijd = integer, julian day on which to interpolate C d_index = integer, index that corresponds to ijd in the C dates and coords arrays C jd = real, julina day on which to interpolate. This C will include the fractional part of the C day. C c_nodes = real, dim(order), coordinates of interpolation c nodes C d_nodes = real, dim(order), dates of interpolation nodes C error = real, dim(3), error estimate for each C coordinate: error(1) = error C estimate for radius, error(2) = C error estimate for latitude, C error(3) = error estimate for C longitude. C C 1/12/2022, Luke Cota: Fixed problem with year rollover. Now when C the year rolls over it continues counting C the day of year, so Jan 1 of the rolled C over year will be day of year 366 (or 367) C C DATE: C February 20, 2020 @ 4:28 PM C C AUTHOR: C Luke Cota C subroutine get_orbit(itp, iYr, Doy, npoints, NMAX, dates, coords, & RVect, error) implicit none C integer, parameter :: MAX_ORDER = 10 C integer, intent(in) :: iYr, NMAX, npoints, itp real, intent(in) :: Doy real, dimension(2,NMAX), intent(in) :: dates real, dimension(3,NMAX), intent(in) :: coords C real, dimension(3), intent(out) :: RVect real, dimension(3), intent(out) :: error C character(len=200) :: fmt1 integer :: date_index, mon, dom integer :: ijd, jd_start integer :: d_index, icur, jcur integer :: date_start, date_end integer :: order, sindex, prev_order logical :: leap real :: Vect, day_cur real, dimension(MAX_ORDER) :: d_nodes real, dimension(MAX_ORDER + 1) :: c_nodes C order = 2 C C Getting the index of the date array that correspond to the C dates that will be used in the interpolation C call get_sindex(NMAX, order, iYr, Doy, dates, sindex) call get_dnodes(NMAX, order, sindex, dates, d_nodes) C C Checking for a year rollover. Luke Cota, 1/12/2022 C call is_leap_year(iYr, leap) if (leap .and. (Doy .ge. 367.0)) then day_cur = Doy - 366.0 else if (.not. leap .and. (Doy .ge. 366.0)) then day_cur = Doy - 365.0 else day_cur = Doy end if C C Interpolating radius and latitude C call poly_interpolate(order, d_nodes, ! radius & coords(1,sindex:sindex + order), & day_cur, error(1), RVect(1)) C call poly_interpolate(order, d_nodes, ! latitude & coords(2,sindex:sindex + order), & day_cur, error(2), RVect(2)) C C Populating the c_node array. This bit is for interpolating the C longitude. We have to account for the fact that the longitude will C reach 360 degrees and then start over at 0 degres. Also, This C corredinate may be better served with a different order interpolation. C Note, everytime we change the interpolation order we need to check if C we have to add more nodes to the d_nodes array C prev_order = order order = 3 if (order .gt. prev_order) then call get_dnodes(NMAX, order, sindex, dates, d_nodes) end if C C Getting the longitude coordinates that will serve as nodes in the C interpolation polynomial C do icur = 1, order c_nodes(icur) = coords(3, sindex + icur - 1) end do C C Scaning interpolation coordinate values to ensure continuity. Looking C to catch when c_nodes(i) = 360 (or there abouts) and c_nodes(i+1) = 1 C (or there abouts) C do icur = 2, order if (c_nodes(icur) < c_nodes(icur - 1)) then c_nodes(icur) = c_nodes(icur) + 360.0 end if end do C C Interpolating the longitude C call poly_interpolate(order, d_nodes, c_nodes(1:order), & day_cur, error(3), RVect(3)) C Need to check if the result of the interpolation is more than 360 C degrees C if (RVect(3) > 360.0) then RVect(3) = RVect(3) - 360.0 end if C for testing only open(unit = 50, name = 'radial_interp_vals.txt', & action ='write', position='append') write(50,"(f8.4,x,f8.4,x,f8.3,x,f9.4)") Doy, RVect(1), RVect(2), & RVect(3) close(50) C Need to reorder RVect Vect = RVect(1) RVect(1) = RVect(3) RVect(3) = Vect return end subroutine get_orbit C C PURPOSE: C This finds the date in the dates array that is closest to, but C no after doy_cur C C INPUTS: C NMAX = int, Max number of dates C order = int, Order of interpolation C iYr = int, Current year C doy_cur = C dates = real, dimensions(2:NMAX), Array of dates: C dates(1:) = year C dates(2:) = doy C C OUTPUTS: C sindex = int, Index of the date that is closest to but not after C the date to be interpolated C C DEPENDENCIES: C None C C NOTES: C 1/12/2022, Luke Cota: Fixed problem with year rollover. Now when C the year rolls over it continues counting C the day of year, so Jan 1 of the rolled C over year will be day of year 366 (or 367) C subroutine get_sindex(NMAX, order, iYr, doy_cur, dates, sindex) implicit none C integer, intent(in) :: NMAX, order, iYr real, intent(in) :: doy_cur real, dimension(2,NMAX), intent(in) :: dates C integer, intent(out) :: sindex C logical :: leap integer :: icur, yr_cur, iiYr real :: ddoy_cur C ddoy_cur = doy_cur iiYr = iYr C C Checking to see if there is a year rollover for doy_cur. Luke Cota, C 1/12/2022 C call is_leap_year(iYr, leap) if (leap .and. (doy_cur .gt. 366.0)) then !if (leap .and. (doy_cur .ge. 367.0)) then ddoy_cur = doy_cur - 366.0 iiYr = iYr + 1 else if (doy_cur .gt. 365.0) then !else if (doy_cur .ge. 366.0) then ddoy_cur = doy_cur - 365.0 iiYr = iYr + 1 end if C do icur = 1, NMAX yr_cur = int(dates(1,icur)) if (ddoy_cur <= dates(2,icur)) then if (iiYr == yr_cur) then exit end if end if end do C sindex = icur - 1 print *, 'ddoy_cur: ', ddoy_cur, ', doy_cur: ', doy_cur, & ', sindex: ', sindex C return end subroutine get_sindex C C PURPOSE: C Populate an array with dates that will be used as nodes in the C date interpolation C C INPUTS: C NMAX = int, Max number of elements in dates C array C order = int, Order of interpolation C sindex = int, Starting index of date nodes in the C dates array C dates = real, dim(2:NMAX), Array of dates C C OUTPUTS: C d_nodes = real, dim(order), Array of day of years that will be C the date nodes in the interpolation. C Note, if there happens to be a year C roll over the next element will just C be 366 or 367 (depending on leap C year) C C NOTES: C subroutine get_dnodes(NMAX, order, sindex, dates, d_nodes) implicit none C integer, intent(in) :: NMAX, order, sindex real, dimension(2,NMAX), intent(in) :: dates C real, dimension(order) :: d_nodes C integer :: icur, yr_initial, ycur logical :: leap real :: dinyr C yr_initial = int(dates(1,sindex)) d_nodes(1) = dates(2,sindex) C do icur = 2, order ycur = int(dates(1, sindex + icur - 1)) C if (ycur > yr_initial) then ! We have a change of year call is_leap_year(yr_initial, leap) dinyr = 365.0 if (leap) dinyr = 366.0 C d_nodes(icur) = dates(2, sindex + icur - 1) + dinyr else d_nodes(icur) = dates(2, sindex + icur - 1) end if end do C return end subroutine get_dnodes C C PURPOSE: C Print error message if the order of interpolation is such that C more data points are needed before the date of interpolation C than are available. C C INPUTS: C iYr = integer, year of date to be interpolated C Doy = real, day of year to be interpolated C order = integer, order of interpolation C C OUTPUTS: C None C C DEPENDENCIES: C None C C DATE: C March 14, 2020 @ 3:40 PM C C AUTHOR: C Luke Cota C subroutine date_range_error(iYr, Doy, order) implicit none integer, intent(in) :: iYr, order real, intent(in) :: Doy C character(len=18) :: fmt1 C print *, 'Error: Cannot interpolate the requested date' fmt1 = "(8X,A7,I4,A8,F5.1)" print fmt1, 'year = ', iYr,', Doy = ', Doy fmt1 = "(8X,A32)" print fmt1, 'There is not enough data for the' fmt1 = "(8X,A13,I2)" print fmt1, 'given order: ',order stop C return end subroutine date_range_error