
>>>>>>>>>>>>>>>>>>> FILE: ../od/od_converge_mod.f90 <<<<<<<<<<<<<<<<<<<

module od_converge_mod

use constants_mod, only: wp, deg_per_rad
use od_mod, only: od_solve, od_statistics, od_prt_sum, od_prt_res, od_msg
use estpar_mod, only: n_param, n_elems, estpar_name, &
   dvx_est, dvx_idx, &
   dvy_est, dvy_idx, &
   dvz_est, dvz_idx, &
   a1_est, a1_idx, &
   a2_est, a2_idx, &
   a3_est, a3_idx, &
   dt_est, dt_idx, &
   aj1_est, aj1_idx, &
   aj2_est, aj2_idx, &
   alf_est, alf_idx, &
   del_est, del_idx, &
   et1_est, et1_idx, &
   et2_est, et2_idx, &
   dth_est, dth_idx, &
   lgk_est, lgk_idx, &
   rho_est, rho_idx, &
   s0_est, s0_idx, &
   n_est_sb, sb_navio_idx, sb_idx
use derivs_mod, only: ng_a1, ng_a2, ng_a3, ng_dt, &
   ng_aj1, ng_aj2, ng_alf, ng_del, ng_et1, ng_et2, ng_dth, &
   yrk_lgk, yrk_rho, sb_mass
use intpar_mod, only: dvx, dvy, dvz

implicit none
save
private

! Public procedures
public :: od_converge, od_converge_def, od_converge_par
! Public variables
public :: converge_msg

!==== Public variables =========================================================
character (len=256) :: converge_msg

!==== Private variables ========================================================
! Verbose control flag
logical :: verbose = .false.

! Default convergence controls
real (kind=wp) :: norm_corr_conv1 = 1.0e-4_wp
real (kind=wp) :: norm_corr_conv2 = 1.0e-2_wp
real (kind=wp) :: delta_rms_conv1 = 1.0e-5_wp
real (kind=wp) :: delta_rms_conv2 = 1.0e-3_wp
real (kind=wp) :: divergence_rate = 1.1
real (kind=wp) :: convergence_rate = 0.9
real (kind=wp) :: max_rms_res = 1.0e5_wp
integer        :: converge_max_iter = 10
integer        :: max_nonconvergent = 4
integer        :: max_partial_step = 3

contains

!########################### OD_CONVERGE #######################################
subroutine od_converge(asteroid, epoch, elements, &
   src_matrix, last_correction, converge_code, iterate, strict_conv, &
   silent, nocorr, new_start, elem_type, full_conv, rms_resid)

logical,                     intent(in)    :: asteroid
real (kind=wp),              intent(in)    :: epoch
real (kind=wp),              intent(inout) :: elements(:)
real (kind=wp),              intent(out)   :: src_matrix(:)
real (kind=wp),              intent(out)   :: last_correction(:)
integer,                     intent(out)   :: converge_code
logical,           optional, intent(in)    :: iterate
logical,           optional, intent(in)    :: strict_conv
logical,           optional, intent(in)    :: silent
logical,           optional, intent(in)    :: nocorr
logical,           optional, intent(in)    :: new_start
character (len=3), optional, intent(in)    :: elem_type
logical,           optional, intent(out)   :: full_conv
real (kind=wp),    optional, intent(out)   :: rms_resid

!===============================================================================
! od_converge computes a new orbit, typically with iterated differential 
!   corrections.
!
! Inputs:
!    asteroid - logical to indicate whether hyperbolic orbit should be 
!       considered a failure
!    epoch - the epoch of the initial_elements (JD)
! Input/Output:
!    elements - the list of estimated parameters. Input as the first guess, 
!       output as the converged solution.
! Output:
!    src_matrix - the sqrt covariance matrix
!    last_correction - the final (hopefully small) corrections applied.
!    converge_code - error code. Zero indicates success, see note below.
! Optional Input:
!    iterate - .true. to converge, .false. to do a single step (Default: .true.)
!    strict_conv - .true. to converge fully, .false. to stop upon 
!       partial convergence. (Default: .true.)
!    silent - .true. to avoid outputting short iteration report 
!       (Default: .false.)
!    nocorr - .true. to skip correcting the elements, including nongravs, etc.
!       (Default: .false.)
!    new_start - .true. to forget convergence information from previous entry.
!       Otherwise, use last RMS, etc. to allow convergence in one iteration.
!       (Default: .false.)
!    elem_type - Type of element set, e.g., CEL or CAR. (Default: CEL)
!
! Optional Output:
!    full_conv - Indicates whether a convergence was full or partial. 
!       (This only makes much sense when strict_conv = .false.)
!    rms_resid - The rms of residuals reported by od_statistics
!
! Notes:
!   - The converge_code has the following interpretation:
!        < 0  - failure
!        = 0  - unqualified success
!        > 0  - qualified success
!     The present convergence codes are 
!        +2  - Not converged, but only single step requested.
!        +1  - RMS paralyzed. (i.e., not changing, but elements still moving.)
!         0  - Converged.
!        -1  - Diverging.
!        -2  - Unreasonably large RMS.
!        -3  - Negative eccentricity. Too many consecutive partial steps.
!        -4  - Negative perihelion. Too many consecutive partial steps.
!        -5  - Hyperbolic asteroid. Too many consecutive partial steps.
!        -6  - Exceeded max_iterations.
!        -7  - Error in lower level routine.
!   - The module variable converge_msg contains diagnostic information.
!
!-------------------------------------------------------------------------------
! Written       Sep 19, 2000       Steve Chesley
! Modified      Feb 22, 2001       Steve Chesley
!      - Allow convergence based on paralyzed RMS.
! Modified      Aug 14, 2001       Steve Chesley
!      - Added silent optional arg.
! Modified      Sep 12, 2001       Steve Chesley
!      - Corrected call to deraup.
! Modified      Oct 04, 2001      Steve Chesley
!      - Added alt_jac argument.
! Modified      Mar 06, 2002       Steve Chesley
!   - Added "only:" restrictions to use statements.
!   - Substantial reworking to support F90 propagator upgrade.
! Modified      May 31, 2002       Steve Chesley
!   - Bug fix: added diff_avail condition to delta_rms convergence test.
! Modified     Jul 02, 2002       Steve Chesley
!   - Added noise to verbose output.
! Modified     Mar 06, 2003       Steve Chesley
!   - Added Yarko stuff.
! Modified     Feb 18, 2004       Steve Chesley
!   - Added nocorr argument to prevent updating of model parameters in derivs
! Modified     Jan 05, 2005       Steve Chesley
!   - Tweaked treatment of full_conv, ensuring that it is set, 
!     even for iter= .false.
! Modified     Apr 26, 2008       Steve Chesley
!   - Forbid negative RHO
! Modified     Nov 24, 2008       Steve Chesley
!   - Add new_start option. Allows second call to od_converge that forgets 
!     earlier convergence info.
!
! $Id: od_converge_mod.f90,v 1.18 2008/11/24 19:33:21 chesley Exp $
!$==============================================================================

!^########################## OD_CONVERGE_DEF ###################################
subroutine od_converge_def

!===============================================================================
! od_converge_def VLDEF's the convergence control parameters for od_converge_mod
!-------------------------------------------------------------------------------
! Written       Dec 14, 2000       Steve Chesley
!$==============================================================================

!^########################## OD_CONVERGE_PAR ###################################
subroutine od_converge_par(verbose_arg)

logical, intent(in) :: verbose_arg

!===============================================================================
! od_converge_par sets the verbose flag for the od_converge_mod, and reads the 
!    various convergence control parameters from the VARLIST into the module.
!-------------------------------------------------------------------------------
! Written       Dec 14, 2000       Steve Chesley
!$==============================================================================

>>>>>>>>>>>>>>>>>>> FILE: ../od/od_mod.f90 <<<<<<<<<<<<<<<<<<<

module od_mod

use constants_mod, only: wp, sin_eps, cos_eps, gm_sun, v_light, two_pi, pi, &
   sec_per_rad, au, rad_per_deg
use obs_mod, only: obs_get, obs_put, obs_get_coords, put_resid, &
   obsmsg, num_obs, max_obs, &
   len_reject, len_obs_type, len_obs_code, len_obs_code_trx, len_bounce_pt, &
   len_time_str, len_frame, len_ra, len_dec, len_observable
use estpar_mod, only: n_param, n_param_p1, n_elems, n_dyn_param, &
   n_con, n_con_dp, n_est_dp, n_est_ndp, n_col, n_utm, n_utm_rhs, &
   last_est_col, estpar_name, apr_sigs, apr_vals, s0_est, s0_idx, &
   n_est_sb, n_est_pl, n_est_gs, sb_navio_idx, sb_idx, pl_idx, gsthry_idx,gs_idx
use get_unit_mod, only: get_unit
use matvec_util_mod, only: vec_mag, vsmat_vec_mul, vsmat_mat_mul, vec_unit, &
   vec_cross, mat_inv
use gram_schmidt_mod, only: gram_schmidt
use stn_pos_mod, only: stn_pos_msg, stn_pos_ini, stn_pos
use eop_mod, only: eoplbl, eoptim
use obscode_mod, only: sites, obscode_index

implicit none
save
private

! Public procedures
public :: od_def
public :: od_ini
public :: od_epoch
public :: od_constraint
public :: od_solve
public :: od_pass
public :: od_filter
public :: od_attrib
public :: od_to_obs
public :: od_get
public :: od_put
public :: od_statistics
public :: od_stat_summ
public :: od_prt_sum
public :: od_prt_res
public :: od_clear

! Public variables
public :: od_msg
public :: od_num_obs
public :: od_t_first, od_t_last, od_use_sim_obs

! Public module variables
character (len=256) :: od_msg           ! Diagnostic message
integer :: od_num_obs                   ! Total number of observations loaded
real (kind=wp) :: od_t_first, od_t_last ! Arc of undeleted obs (JD)
logical :: od_use_sim_obs = .false.     ! Flag to indicate use of simulated obs

!===============================================================================
! This is the od_mod module. It really is.
!
! Public routines: 
!    od_def        - Defines od_mod VARLIST variables.
!    od_ini        - Loads the observational data, among other things.
!    od_epoch      - Computes the weighted midpoint of the observed arc, 
!                    rounded to the nearest day.
!    od_solve      - Computes the differential correction.
!    od_to_obs     - Dumps relevant od_mod data back into the obs_mod.
!    od_get        - Retrieves various arrays from the od_mod for use outside.
!    od_put        - Loads various arrays from outside into the od_mod.
!    od_statistics - Computes miscellaneous statistical information after 
!                    an iteration.
!    od_stat_summ  - Retrieves misc. statistical information for use outside 
!                    the od_mod.
!    od_prt_sum    - Prints a short report on the fit statistics.
!    od_prt_res    - Prints the residual report.
!
! Public variables:
!    od_msg     - Diagnostic message in case of error.
!    od_num_obs - The number of observations loaded by od_ini.
!    od_t_first - The JD of the first accepted observation.
!    od_t_last  - The JD of the last accepted observation.
!    od_use_sim_obs - Flag to indicate use of simulated obs. 
!
! Notes:
!     - The relevant data for radar obs is stored in the right ascension areas.
!     - A '*' in column 80 below indicates data that is typically written back 
!       to the obs file.
!
!-------------------------------------------------------------------------------
! Written       Sep 19, 2000       Steve Chesley
! Modified      Mar 02, 2001       Steve Chesley
!      - Added ACCEPT_ALL option.
! Modified      Jul 30, 2001      Steve Chesley
!     - Added arrays for Earth receive time delay for S/C obs.
! Modified      Sep 21, 2001      Steve Chesley
!     - Added a priori constraint module variables and VARLIST options.
! Modified      Mar 06, 2002       Steve Chesley
!   - Added "only:" restrictions to use statements.
!   - Modified to work with F90 propagator routines.
! Modified      Sep 24, 2003       Steve Chesley
!     - Add simulated observations.
! Modified      Nov 19, 2003       Steve Chesley
!     - Add chi_square
! Modified      Mar 12, 2004       Steve Chesley
!     - New routines: od_constraint, od_pass, od_filter, od_attrib.
!     - And ancillary module variables.
! Modified      Jan 06, 2005       Steve Chesley
!     - Removed maximum residual checking as this is done in od_prt_sum
!
! $Id: od_mod.f90,v 1.13 2007/06/14 22:54:06 chesley Exp $
!$==============================================================================

!^==============================================================================
subroutine od_def

!-------------------------------------------------------------------------------
! od_def:  Define all VARLIST variables unique to the od_mod.
!$------------------------------------------------------------------------------

!^ Files included in the od_mod
include "od_ini.f90"
include "od_epoch.f90"
include "od_solve.f90"
include "od_constraint.f90"
include "od_pass.f90"
include "od_filter.f90"
include "od_attrib.f90"
include "od_to_obs.f90"
include "od_get.f90"
include "od_put.f90"
include "od_statistics.f90"
include "od_stat_summ.f90"
include "od_prt_sum.f90"
include "od_prt_res.f90"
include "od_clear.f90"
!$

>>>>>>>>>>>>>>>>>>> FILE: ../od/od_ini.f90 <<<<<<<<<<<<<<<<<<<

subroutine od_ini(j2000_arg, eph_start, eph_stop, time_ref, &
   solve_verbose, radius_body, errcod, file_num)

logical,           intent(in)  :: j2000_arg
real (kind=wp),    intent(in)  :: eph_start
real (kind=wp),    intent(in)  :: eph_stop
real (kind=wp),    intent(in)  :: time_ref
logical,           intent(in)  :: solve_verbose
real (kind=wp),    intent(out) :: radius_body
integer,           intent(out) :: errcod
integer, optional, intent(in)  :: file_num

!===============================================================================
! od_ini imports the observational data using the specified file index. 
!    This must be called only after a call to obs_ini has loaded an 
!    observation file into the obs_mod.
!
!    The routine does a good deal of housekeeping:
!
!     a) Parses the parameter list to determine what and how many parameters
!        are to be estimated.
!     b) Initializes the center of light offset model, after reading the 
!        relevant VARLIST variables.
!     c) Initializes the Earth orientation routines (DELINI, RDRSYS, STNSYS).
!     d) Computes the radius body. This  is generally taken from the VARLIST 
!        (RADB) if present, otherwise it is set to zero. However, if there 
!        are PP (peak power) radar observations and radius_body is found to 
!        be zero, then it will compute the value based on the VARLIST H value 
!        and an assumed albedo.
!     e) Initializes the radar observation processing routine (RDRINI).
!     f) Counts the number of relevant observational data from the obs_mod, 
!        and allocates the requisite arrays accordingly. Observations 
!        marked "ignore" (#), and observations outside of the OBS_START and 
!        OBS_STOP time span are ignored.
!     g) Initializes the obs_mod lookup table (obs_index), and the number 
!        of observations (od_num_obs).
!     h) Loads the observational data, including station location. By default, 
!        the more accurate RDRSTN routine is used for _all_ observations. 
!        The less accurate STNPOS routine is used only if the observation 
!        is optical and either OPT_USE_EOP is set .true. or RDRSTN fails 
!        due to insufficient data.
!     i) Undeletes all photometry marked "D" if ACCEPT_ALL is set in VARLIST.
!     j) Sets the values of od_t_first and od_t_last.
!
! Inputs:
!   j2000_arg - Reference system flag. True for J2000, false for B1950.
!   time_ref  - Reference epoch (JD)
!   eph_start - Start time of the available ephemerides
!   eph_stop  - Stop time of the available ephemerides
!   solve_verbose - Turns on verbose output.
! Output:
!   radius_body - radius of the body
!   errcod   - Error code. zero indicates success.
! Optional input:
!   file_num - Specifies the internal obs_manager file (default=1)
!
! Notes:
!   - The relevant data for radar obs is stored in the right ascension area.
!
! TODO:
!     - Correctly handle FK5 vs. ICRS differences.
!-------------------------------------------------------------------------------
! Written       Jan 05, 2001       Steve Chesley
! Modified      Mar 02, 2001       Steve Chesley
!      - Added ACCEPT_ALL option.
! Modified      Apr 10, 2001       Steve Chesley
!      - Fixed bug in type 'S' observation handling.
! Modified      May 24, 2001       Steve Chesley
!      - Fixed (bad) bug that caused pre-TPM obs to use geocentric observer.
! Modified      Jul 30, 2001      Steve Chesley
!     - Added special handling of simulated Earth receive time for S/C obs.
! Modified      Sep 21, 2001      Steve Chesley
!     - Added a priori constraints.
! Modified      Oct xx, 2001      Steve Chesley
!     - Times from the obs_mod are now wrt t_ref. Removed de-referencing.
! Modified      Feb 06, 2002      Steve Chesley
!     - Implemented new station position routines based on EOP files.
! Modified      Feb 25, 2002      Steve Chesley
!     - Implemented new propagator (estpar_mod, derivs_mod, etc.)
!     - Disabled SRI file loading.
! Modified      May 21, 2002      Steve Chesley
!     - Print A Priori info only if verbose==.true.
! Modified      Feb 11, 2003       Steve Chesley
!     - Save unit scaling into cos_dec array for radar too.
! Modified      Sep 17, 2003       Steve Chesley
!     - Tweak documentation.
! Modified      Sep 24, 2003       Steve Chesley
!     - Add simulated observations.
! Modified      Oct 13, 2003       Steve Chesley
!     - Changed sim_obs type character from "F" to "+".
! Modified      Jan 23, 2004       Steve Chesley
!     - Change stn_pos error handling for radar.
! Modified      Mar 12, 2004       Steve Chesley
!     - Changed x_rec from a 6-vector to a 3-vector.
!     - Changes associated with od_filter routine
! Modified      Jan 06, 2005       Steve Chesley
!     - Minor bug fixes
! Modified      Jul 01, 2005       Steve Chesley
!     - Only call stn_pos_ini once per invocation.
! Modified      Jun 14, 2007       Steve Chesley
!     - Minor mods for Sim. bs. handling
!     - Allow APR_SRI files to be used
!
! $Id: od_ini.f90,v 1.17 2007/06/14 22:54:06 chesley Exp $
!$==============================================================================

>>>>>>>>>>>>>>>>>>> FILE: ../od/od_epoch.f90 <<<<<<<<<<<<<<<<<<<

subroutine od_epoch(epoch, errcod)

real (kind=wp), intent(out) :: epoch
integer,        intent(out) :: errcod

!===============================================================================
! od_epoch computes the weighted midpoint of the observed arc, rounded to 
!    the nearest day.
!
! Output:
!    epoch  - The epoch of the initial_elements (JED).
!    errcod - error code. Zero indicates success.
! 
! Note:
!    - Radar observations are not included in the computation.
!
!-------------------------------------------------------------------------------
! Written       Steve Chesley Sep. 20, 2000
! Modified      Steve Chesley Mar. 03, 2005
!     - Tweaked documentation. (Output is in JED.)
!
! $Id: od_epoch.f90,v 1.4 2007/06/14 21:42:14 chesley Exp $
!$==============================================================================

>>>>>>>>>>>>>>>>>>> FILE: ../od/od_constraint.f90 <<<<<<<<<<<<<<<<<<<

subroutine od_constraint(n_constr, len_constr, constr, errcod, alt_jac)

integer, intent(in) :: n_constr
integer, intent(in) :: len_constr
real (kind=wp), intent(in) :: constr(len_constr, n_constr)
integer, intent(out) :: errcod
real (kind=wp), optional, intent(in) :: alt_jac(:,:)

!===============================================================================
! od_constraint sets the constraint matrix for use in od_filter. 
!
! Input:
!   constraint - The constraint matrix. The columns are the multiple 
!     constraints, and the number of columns is the number of constraints.
!
! Output:
!   errcod   - error code. Zero indicates success.
!-------------------------------------------------------------------------------
! Written       Mar 12, 2004       Steve Chesley
!
! $Id: od_constraint.f90,v 1.1 2004/12/01 20:49:27 chesley Exp $
!$==============================================================================

>>>>>>>>>>>>>>>>>>> FILE: ../od/od_solve.f90 <<<<<<<<<<<<<<<<<<<

subroutine od_solve(epoch_jd, input_elements, correction, src_matrix, norm_dq, &
   errcod, prt_filter, elem_type)

real (kind=wp),             intent(in)  :: epoch_jd
real (kind=wp),             intent(in)  :: input_elements(:)
real (kind=wp),             intent(out) :: correction(:)
real (kind=wp),             intent(out) :: src_matrix(:)
real (kind=wp),             intent(out) :: norm_dq
integer,                    intent(out) :: errcod
logical,          optional, intent(in)  :: prt_filter
character(len=3), optional, intent(in)  :: elem_type

!===============================================================================
! od_solve computes the differential correction given an orbit and an 
!    observation set.
!
! Input:
!    epoch_jd - the epoch of the initial_elements (JD)
!    input_elements - the initial elements.
! Output:
!    correction - the desired correction to the initial elements, in the same 
!       space as the elements.
!    src_matrix - the vector_stored sqrt covariance matrix
!    norm_dq    - The RMS of the normalized corrections: 
!       sqrt( corr^T . C . corr / n_param), where C is the normal matrix.
!    errcod     - error code. Zero indicates success.
! Optional Inputs:
!    prt_filter - output the sqrt info and sqrt cov matrices to STDOUT.
!       (Default: .false.)
!    elem_type - Type of element set, e.g., CEL or CAR. (Default: CEL)
!
!-------------------------------------------------------------------------------
! Written       Jan 04, 2001      Steve Chesley
! Modified      Mar 05, 2001      Steve Chesley
!     - Optical residuals computed modulo two_pi.
! Modified      Mar 07, 2001      Steve Chesley
!     - Revised logic for forward and backward integration loop.
! Modified      Mar 20, 2001      Steve Chesley
!     - Added constrained correction capability.
! Modified      Mar 26, 2001      Steve Chesley
!     - Improved radar bounce point logic to include warning if undefined.
! Modified      Apr 11, 2001      Steve Chesley
!     - Added calls to objsdt to deal with epoch between emit and obs times.
!     - Changed computation of max_ltime to deal with nonelliptic orbits.
! Modified      Jul 30, 2001      Steve Chesley
!     - Added code to compute Earth receive time delay for S/C obs.
! Modified      Aug 10, 2001      Steve Chesley
!     - Fixed bug dealing with S0 partials
! Modified      Sep 12, 2001      Steve Chesley
!     - Updated matrix multiplication to use vs_mat_mul.
! Modified      Sep 21, 2001      Steve Chesley
!     - Added a priori constraints.
! Modified      Oct 04, 2001      Steve Chesley
!     - Added alt_jac argument
! Modified      Nov 14, 2001      Steve Chesley
!     - Tweaked dimensioning of x_ini and x, and call to cetpv
! Modified      Dec 10, 2001      Steve Chesley
!     - Changed apriori_tmp to allocatable. Added allocate and deallocate 
!       statements.
! Modified      Feb 20, 2002      Steve Chesley
!     - Corrected computation of max_ltime.
! Modified      Feb 25, 2002      Steve Chesley
!     - Substantial reworking IAW upgraded F90 propagator.
!     - Consider parameters not supported yet.
! Modified      May 21, 2002      Steve Chesley
!     - Consider parameters supported.
! Modified      Sep 24, 2003       Steve Chesley
!     - Add simulated observations.
! Modified      Mar 03, 2004       Steve Chesley
!     - Split the routine into two parts: od_pass & od_filter
! Modified      Nov 20, 2008       Steve Chesley
!     - Changed estimation of max_owlt for parabolic/hyperbolic orbits
!     - Abort if state at epoch shows r > 1500 AU or v > 0.2c
!
! $Id: od_solve.f90,v 1.19 2008/11/21 01:49:42 chesley Exp $
!$==============================================================================

>>>>>>>>>>>>>>>>>>> FILE: ../od/od_pass.f90 <<<<<<<<<<<<<<<<<<<

subroutine od_pass(epoch_jd, x_ini, max_owlt, sri_matrix, errcod)

real (kind=wp), intent(in)  :: epoch_jd
real (kind=wp), intent(in)  :: x_ini(:,:)
real (kind=wp), intent(in)  :: max_owlt
real (kind=wp), intent(out) :: sri_matrix(:)
integer,        intent(out) :: errcod

!===============================================================================
! od_pass computes the residuals and forms the SRI Matrix with embedded RHS
!
! Input:
!    epoch_jd - the epoch of the initial_elements (JD)
!    x_ini - Initial cartesian state and partials wrt initial conditions.
!            (The partials could be wrt elements or state.)
! Output:
!    sri_matrix - the vector_stored sqrt information matrix with RHS
!    errcod     - error code. Zero indicates success.
!
!-------------------------------------------------------------------------------
! Written       Mar 03, 2004      Steve Chesley
! Modified      Jan 07, 2005       Steve Chesley
!     - Changed the handling of residuals for simulated observations
! Modified      Nov 20, 2008       Steve Chesley
!     - Changed some variable names to be more clear
!     - Fixed bug related to separation of integrator state and observation
!       state in OBJSTA loop
!     - Properly handle cases where geocentric light time correction 
!       is not accurate enough
!
! $Id: od_pass.f90,v 1.5 2008/11/21 01:54:22 chesley Exp $
!$==============================================================================

>>>>>>>>>>>>>>>>>>> FILE: ../od/od_filter.f90 <<<<<<<<<<<<<<<<<<<

subroutine od_filter(sri_matrix, correction, src_matrix, norm_dq, errcod, &
   prt_filter, elements)

real (kind=wp),           intent(in)  :: sri_matrix(:)
real (kind=wp),           intent(out) :: correction(:)
real (kind=wp),           intent(out) :: src_matrix(:)
real (kind=wp),           intent(out) :: norm_dq
integer,                  intent(out) :: errcod
logical,        optional, intent(in)  :: prt_filter
real (kind=wp), optional, intent(in)  :: elements(:)

!===============================================================================
! od_filter computes the differential correction given an SRI matrix 
!
! Input:
!    sri_matrix - the vector_stored sqrt information matrix with RHS
! Output:
!    correction - the desired correction to the initial elements, in the same 
!       space as the elements.
!    src_matrix - the vector_stored sqrt covariance matrix
!    norm_dq    - The RMS of the normalized corrections: 
!       sqrt( corr^T . C . corr / n_param), where C is the normal matrix.
!    errcod     - error code. Zero indicates success.
! Optional Inputs:
!    prt_filter - output the sqrt info and sqrt cov matrices to STDOUT.
!       (Default: .false.)
!    elements - Input elements. Used only to compute RHS for a priori contraints
!
! Note: If constraints are set then od_filter computes the least squares 
!       correction in the subspace orthogonal to the constraints. If alt_jac
!       is loaded then the constraint must be ALREADY in the alternate space. 
!       (Although the correction returned will be in the working space.)
!
!-------------------------------------------------------------------------------
! Written       Mar 12, 2004      Steve Chesley
!
! $Id: od_filter.f90,v 1.1 2004/12/01 20:49:27 chesley Exp $
!$==============================================================================

>>>>>>>>>>>>>>>>>>> FILE: ../od/od_attrib.f90 <<<<<<<<<<<<<<<<<<<

subroutine od_attrib(epoch, attrib, e_0, x_obs, phi, errcod)

real (kind=wp), intent(out) :: epoch
real (kind=wp), intent(out) :: attrib(4)
real (kind=wp), intent(out) :: e_0(3)
real (kind=wp), intent(out) :: x_obs(6)
real (kind=wp), intent(out) :: phi
integer,        intent(out) :: errcod

!===============================================================================
! od_attrib computes the "attributable" associated with the observations
!
! Output:
!    epoch  - The epoch of the initial_elements after t_ref.
!    errcod - error code. Zero indicates success.
! 
! Note:
!    - Radar observations are not included in the computation.
!
!-------------------------------------------------------------------------------
! Written       Steve Chesley Sep. 20, 2000
!
! $Id: od_attrib.f90,v 1.1 2004/12/01 20:49:27 chesley Exp $
!$==============================================================================

>>>>>>>>>>>>>>>>>>> FILE: ../od/od_to_obs.f90 <<<<<<<<<<<<<<<<<<<

subroutine od_to_obs(errcod)

integer,           intent(out) :: errcod

!===============================================================================
! od_to_obs writes the rejection flags and residuals back to the obs_mod. 
!    This is typically done after a new and final orbit solution is obtained 
!    so that the obs file will reflect the postfit residuals and rejection 
!    flags of the orbit that is in the orb file.
!
! Output:
!   errcod   - error code. Zero indicates success.
!
! Notes:
!    - There must be a call to od_statistics before the call to od_to_obs
!      so that the RSS residual values are current.
!-------------------------------------------------------------------------------
! Written       Sep 18, 2000       Steve Chesley
! Modified      Jul 30, 2001      Steve Chesley
!     - Added code to obs_put Earth receive time delays for S/C obs.
! Modified      Sep 19, 2003      Steve Chesley
!     - Short return if in sim obs mode.
!
! $Id: od_to_obs.f90,v 1.7 2004/12/01 20:49:27 chesley Exp $
!$==============================================================================

>>>>>>>>>>>>>>>>>>> FILE: ../od/od_get.f90 <<<<<<<<<<<<<<<<<<<

subroutine od_get(errcod, od_obs_index, od_reject, od_obs_type, od_obs_code, &
   od_time_ut, od_x_rec, od_ra_sig, od_dec_sig, od_x_obj, od_ra_res, &
   od_dec_res, od_ra, od_dec, od_ra_partial, od_dec_partial, &
   od_cos_dec, od_ra_wgt, od_dec_wgt, od_src_constr)

integer,                                intent(out) :: errcod
integer,                      optional, intent(out) :: od_obs_index(:)
character (len=len_reject),   optional, intent(out) :: od_reject(:)
character (len=len_obs_type), optional, intent(out) :: od_obs_type(:)
character (len=len_obs_code), optional, intent(out) :: od_obs_code(:)
real (kind=wp),               optional, intent(out) :: od_time_ut(:)
real (kind=wp),               optional, intent(out) :: od_x_rec(:,:)
real (kind=wp),               optional, intent(out) :: od_ra_sig(:)
real (kind=wp),               optional, intent(out) :: od_dec_sig(:)
real (kind=wp),               optional, intent(out) :: od_x_obj(:,:)
real (kind=wp),               optional, intent(out) :: od_ra_res(:)
real (kind=wp),               optional, intent(out) :: od_dec_res(:)
real (kind=wp),               optional, intent(out) :: od_ra(:)
real (kind=wp),               optional, intent(out) :: od_dec(:)
real (kind=wp),               optional, intent(out) :: od_ra_partial(:,:)
real (kind=wp),               optional, intent(out) :: od_dec_partial(:,:)
real (kind=wp),               optional, intent(out) :: od_cos_dec(:)
real (kind=wp),               optional, intent(out) :: od_ra_wgt(:)
real (kind=wp),               optional, intent(out) :: od_dec_wgt(:)
real (kind=wp),               optional, intent(out) :: od_src_constr(:)

!===============================================================================
! od_get retrieves selected arrays from the od_mod private module data.
!
! Output:
!   errcod   - Error code. zero indicates success.
! Optional output:
!   All of the arrays in the list above.
! Notes: 
!    - Normally the calling program will get the od_num_obs and then 
!        allocate the target arrays before calling this routine.
!    - The size of the target array is not compared to the size module array.
!    - The od_stat flag is not actually necessary to retrieve od_src_constr, 
!        but, since od_stat_flag implies od_solve_flag, the requirement is 
!        conservative, and simpler to implement.
!-------------------------------------------------------------------------------
! Written       Sep 18, 2000       Steve Chesley
! Modified      Mar 12, 2004       Steve Chesley
!    - Added od_src_constr optional argument.
!
! $Id: od_get.f90,v 1.3 2004/12/01 20:49:27 chesley Exp $
!$==============================================================================

>>>>>>>>>>>>>>>>>>> FILE: ../od/od_put.f90 <<<<<<<<<<<<<<<<<<<

subroutine od_put(od_reject, od_ra, od_dec, od_cos_dec, od_ra_wgt)

character (len=len_reject),   optional, intent(in) :: od_reject(:)
real (kind=wp),               optional, intent(in) :: od_ra(:)
real (kind=wp),               optional, intent(in) :: od_dec(:)
real (kind=wp),               optional, intent(in) :: od_cos_dec(:)
real (kind=wp),               optional, intent(in) :: od_ra_wgt(:)

!===============================================================================
! od_put puts the selected array back into the od_mod private data repository.
!
! Optional input:
!   od_reject  - Rejection character codes
!   od_ra      - Right ascencion
!   od_dec     - Declination
!   od_cos_dec - Cosine(dec) * sec_per_rad
!   od_ra_wgt  - cos_dec / ra_sig
!
! Notes:
!    - Normally the calling program will get the od_num_obs and then 
!      allocate the target arrays before calling this routine.
!    - The size of the target array is not compared to the size module array.
!-------------------------------------------------------------------------------
! Written       Sep 18, 2000       Steve Chesley
!
! $Id: od_put.f90,v 1.3 2004/12/01 20:49:27 chesley Exp $
!$==============================================================================

>>>>>>>>>>>>>>>>>>> FILE: ../od/od_statistics.f90 <<<<<<<<<<<<<<<<<<<

subroutine od_statistics(correction, rms_res, max_res_diff, rd_flag, errcod)

real (kind=wp), intent(in)  :: correction(:)
real (kind=wp), intent(out) :: rms_res
real (kind=wp), intent(out) :: max_res_diff
logical,        intent(out) :: rd_flag
integer,        intent(out) :: errcod

!===============================================================================
! od_statistics computes numerous statistical values and auxiliary data.
!
! Input:
!   correction   - the most recently applied correction.
! Output:
!   rms_res      - The RMS of normailzed residuals.
!   max_res_diff - The maximum computed residual difference.
!   rd_flag      - Flag to indicate whether residual differences were computed.
!   errcod       - error code. Zero indicates success.
!
! Notes:
!   - The od_solve_flag must be set by a call to od_solve before entering 
!     this routine.
!   - The computation of residual differences is governed by the status of 
!     res_pred_flag. If this flag is not set then the predicted residuals 
!     were not computed in the previous pass thus the residual differences 
!     cannot be computed for this pass.
!-------------------------------------------------------------------------------
! Written       Sep 20, 2000       Steve Chesley
! Modified      Apr 05, 2001       Steve Chesley
!   - Change to permit satellite only data arcs.
! Modified      Jul 02, 2002       Steve Chesley
!   - Changed mean_res terms to reflect mean of weighted residuals.
! Modified      Feb 11, 2003       Steve Chesley
!   - Fixed bug in computation of residual differences
! Modified      Sep 19, 2003       Steve Chesley
!   - Tweaked max residual search to ignore sim obs.
! Modified      Nov 19, 2003       Steve Chesley
!     - Add chi_square
! Modified      Jan 06, 2005       Steve Chesley
!     - Removed maximum residual checking as this is done in od_prt_sum
!
! $Id: od_statistics.f90,v 1.8 2005/01/06 20:47:07 chesley Exp $
!$==============================================================================

>>>>>>>>>>>>>>>>>>> FILE: ../od/od_stat_summ.f90 <<<<<<<<<<<<<<<<<<<

subroutine od_stat_summ(n_opt_sat, n_del, n_dop, rms_vec, errcod)

integer,        intent(out) :: n_opt_sat
integer,        intent(out) :: n_del
integer,        intent(out) :: n_dop
real (kind=wp), intent(out) :: rms_vec(*)
integer,        intent(out) :: errcod

!===============================================================================
! od_stat_summ returns summary data from the od_statistics call.
!
! Output:
!   n_opt_sat - Number of optical observations, including satellite obs.
!   n_del     - Number of radar delay observations.
!   n_dop     - Number of radar Doppler observations.
!   rms_vec   - Vector of various RMS values:
!                  1: Unweighted, RA (arcsec)
!                  2: Unweighted, Dec (arcsec)
!                  3: Unweighted, Total optical (arcsec)
!                  4: Weighted, Total optical
!                  5: Normalized, optical only
!                  6: Normalized, delay
!                  7: Normalized, Doppler
!                  8: Normalized, optical+radar
!   errcod    - Zero indicates success.
!
! Notes:
!   - The od_stat_flag must be set by a call to od_statistics before entering 
!     this routine.
!
!-------------------------------------------------------------------------------
! Written       Nov 13, 2000       Steve Chesley
! Modified      Nov 15, 2001       Steve Chesley
!   - If (.not.od_stat_flag) then call od_statistics, rather then returning 
!     with error.
! Modified      Jan 25, 2002       Steve Chesley
!   - Documentation.
! Modified      Mar 06, 2002       Steve Chesley
!   - Changed varaible n_p to n_elems based on new F90 propagator routines.
!
! $Id: od_stat_summ.f90,v 1.5 2004/12/01 20:49:27 chesley Exp $
!$==============================================================================

>>>>>>>>>>>>>>>>>>> FILE: ../od/od_prt_sum.f90 <<<<<<<<<<<<<<<<<<<

subroutine od_prt_sum(errcod, object, soln_num, producer, pleph,&
   sbeph, small_bodies, num_small_bodies)

integer,                     intent(out):: errcod
character (len=*), optional, intent(in) :: object
character (len=*), optional, intent(in) :: soln_num
character (len=*), optional, intent(in) :: producer
character (len=*), optional, intent(in) :: pleph
character (len=*), optional, intent(in) :: sbeph
character (len=*), optional, intent(in) :: small_bodies(*)
integer,           optional, intent(in) :: num_small_bodies

!===============================================================================
! od_prt_res_sum prints to STDOUT a summary of a particular orbit solution.
!
! Output:
!   errcod   - error code. Zero indicates success.
!
! Optional Inputs: (All self-explanatory)
!   object
!   soln_num
!   producer
!   pleph
!   sbeph
!   small_bodies
!   num_small_bodies
!-------------------------------------------------------------------------------
! Written       Sep 20, 2000       Steve Chesley
! Modified      Apr 11, 2001       Steve Chesley
!     - Improved reporting if for objects wthout ground-based optical 
!       observations.
! Modified      Jul 31, 2001       Steve Chesley
!     - Minor bug fix in satellite section.
! Modified      Jul 02, 2002       Steve Chesley
!     - Changed "Mean" to "Mean weighted residual" to reflect change in 
!       output from od_statisitics.
! Modified      Sep 24, 2003       Steve Chesley
!     - Output number if simulated observations
! Modified      Nov 19, 2003       Steve Chesley
!     - Add chi_square
! Modified      Feb 01, 2006       Steve Chesley
!     - Increased width of field for no. sat. obs.
!
! $Id: od_prt_sum.f90,v 1.9 2007/06/14 21:46:13 chesley Exp $
!$==============================================================================

>>>>>>>>>>>>>>>>>>> FILE: ../od/od_prt_res.f90 <<<<<<<<<<<<<<<<<<<

subroutine od_prt_res(unit, errcod, system, wrt_del)

integer,                     intent(in)  :: unit
integer,                     intent(out) :: errcod
character (len=*), optional, intent(in)  :: system
logical,           optional, intent(in)  :: wrt_del

!===============================================================================
! od_prt_res prints the Sqrt Information and Sqrt Covariance matrices 
!     and the conditioning number to STDOUT.
!
! Inputs:
!   unit     - Output file unit. Headers are only printed if unit == 6.
! Output:
!   errcod   - error code. Zero indicates success.
! Optional inputs:
!   system   - Reference system appears in the header. (Normally J2000 or B1950)
!   wrt_del  - Determines whether deleted observations should be printed.
!-------------------------------------------------------------------------------
! Written       Sep 20, 2000       Steve Chesley
! Modified      Mar 06, 2001       Steve Chesley
!    - Made system optional.
!    - Added wrt_del optional argument.
! Modified      Jan 25, 2002       Steve Chesley
!    - Included optional system ('sys') in header.
! Modified      Sep 24, 2003       Steve Chesley
!     - Add simulated observations.
!
! $Id: od_prt_res.f90,v 1.6 2004/12/01 20:49:27 chesley Exp $
!$==============================================================================

>>>>>>>>>>>>>>>>>>> FILE: ../od/od_clear.f90 <<<<<<<<<<<<<<<<<<<

subroutine od_clear(errcod)

integer, intent(out) :: errcod

!===============================================================================
! od_clear clears the od_mod data in anticipation of another call to od_ini. 
!    This routine simply deallocates the allocatable arrays and updates the 
!    od_ini_flag. This call is not needed for the first call to od_ini, 
!    but is needed before every subsequent call. 
!
! Output:
!   errcod   - error code. Zero indicates success.
!-------------------------------------------------------------------------------
! Written       Dec 15, 2000       Steve Chesley
! Modified      Feb 26, 2002       Steve Chesley
!    - Added ert_delay to deallocate statement.
! Modified      Mar 06, 2002       Steve Chesley
!    - Removed p_name from deallocate statement.
! Modified      Jul 02, 2005       Steve Chesley
!    - Added src_matrix2 and rot_constr to deallocate statement.
!
! $Id: od_clear.f90,v 1.6 2007/06/14 21:41:29 chesley Exp $
!$==============================================================================

>>>>>>>>>>>>>>>>>>> FILE: ../od/mag_mod.f90 <<<<<<<<<<<<<<<<<<<

module mag_mod

use constants_mod, only: wp
use obs_mod, only: obs_get, put_resid, obsmsg, &
   len_magnitude, len_mag_reject, len_reject, len_obs_type
use od_mod, only: od_get, od_msg, od_num_obs

implicit none
save
private

! Public procedures
public :: mag_est, mag_est_com

! Public variables
public :: mag_msg

! Public module variables
character (len=256) :: mag_msg

!=============================================================================
! This is the mag_mod module.
!
! Public routines: 
!    mag_est - Estimate the asteroid's magnitude.
!    mag_est_com - Estimate the comet's magnitudes.
! Public variables:
!    mag_msg - Diagnostic message in case of error.
! Notes:
!
! TODO: For outlier handling, add runtime options and reporting. 
!-----------------------------------------------------------------------------
! Written       Oct 22, 2000       Steve Chesley
! Modified      Mar 02, 2001       Steve Chesley
!      - Fixed Misc. bugs.
!      - Added verbose and accept_all arguments.
! Modified      Mar 06, 2002       Steve Chesley
!   - Added "only:" restrictions to use statements.
! Modified      Nov 19, 2003       Steve Chesley
!   - Eliminated "Unknown color" STDERR output.
! Modified      Dec 18, 2003       Steve Chesley
!   - Handle v, J and H bands.
! Modified      Sep 20, 2006       Nick Mastrodemos
!   - Add cometary magnitudes
! Modified      Mar 31, 2005       Steve Chesley
!   - Add "_wp" to inline constants in comp_mag
! Modified      Nov 29, 2007       Steve Chesley
!   - Add 'r' band with V-r = 0.35 (per email from Spahr)
!   - And changed handling of unknown pass bands.
! Modified      Feb 13, 2008       Steve Chesley
!   - Added new pass band and aligned all offsets with the MPC (per Spahr email)
! Modified      May 06, 2008       Steve Chesley
!   - Changed outlier testing to use absolute value.
!   - Changed handling for unknown colors.
!
! $Id: mag_mod.f90,v 1.17 2008/05/06 20:10:27 chesley Exp $
!===============================================================================

contains

subroutine mag_est(g, abs_mag, rms_mag, errcod, update_obs, file_num, &
   verbose, accept_all)

real (kind=wp),    intent(in)  :: g
real (kind=wp),    intent(out) :: abs_mag
real (kind=wp),    intent(out) :: rms_mag
integer,           intent(out) :: errcod 
logical, optional, intent(in)  :: update_obs
integer, optional, intent(in)  :: file_num
logical, optional, intent(in)  :: verbose
logical, optional, intent(in)  :: accept_all

!=============================================================================
! mag_est
!
! Input:
!   g - Slope parameter.
! Output:
!   abs_mag  - The estimated absolute magnitude.
!   errcod   - Error code. Zero indicates success.
! Optional input:
!   update_obs - flag to determine whether or not to dump results to obs_mod.
!   file_num   - Specifies the internal obs_manager file (default=1)
!   verbose    - Highly verbose outlier rejection reporting
!   accept_all - Force initial acceptance of all photometry marked "D"
!
! Notes:
!    - This must be called _after_ the call to od_ini and obs_ini.
!
!$============================================================================

>>>>>>>>>>>>>>>>>>> FILE: ../od/reject_mod.f90 <<<<<<<<<<<<<<<<<<<

module reject_mod

use constants_mod, only: wp, sec_per_rad
use obs_mod, only: len_reject, len_obs_type, len_obs_code
use od_mod, only: od_get, od_put, od_msg, od_num_obs
use matvec_util_mod, only: vsmat_vec_mul
use estpar_mod, only: n_elems

implicit none
save
private

! Public procedures
public :: reject_ini, reject
! Public variables
public :: reject_msg

character (len=128) :: reject_msg

!===============================================================================
! This is the reject_mod module, used for rejection and recovery of 
! outlying observations.
!
! Public routines: 
!    reject_ini - set up variables
!    reject - rejects/recovers outlier observation
!
! Public variables:
!    reject_msg - Diagnostic message in case of error.
!
! NOTE:
!    - There are several arrays that are allocated by reject_ini, 
!      and are never deallocated.
!
! REFERENCE:
!   Carpino, Milani, and Chesley
!   "Error Statistics of Asteroid Optical Astrometric Observations"
!   Icarus, 2001
!
!-------------------------------------------------------------------------------
! Written       Sep 27, 2000       Steve Chesley
! Modified      Jan 4, 2001        Steve Chesley
!    - Implemented orbital error in chi^2 computation
! Modified      May 23, 2001        Steve Chesley
!    - Implemented check to ensure minimum number of obs needed for OD.
!    - Improved documentation.
!    - Implemented check for rejection of last observation in apparition.
!    - Implemented rej_opp argument of reject_ini.
! Modified      May 24, 2001        Steve Chesley
!    - Tweaked obs_avail logic.
! Modified      Feb 04, 2002        Steve Chesley
!    - Fixed carriage control bugs.
! Modified      Feb 15, 2002        Steve Chesley
!    - Add more precision to date reporting.
! Modified      Mar 06, 2002       Steve Chesley
!   - Added "only:" restrictions to use statements.
! Modified      Feb 12, 2003       Steve Chesley
!   - Issue warning when rejecting radar astrometry.
! Modified      Dec 07, 2004       Steve Chesley
!   - Remove obsmsg from 'use obs_mod, only:' list.
!
! $Id: reject_mod.f90,v 1.14 2005/01/06 20:47:07 chesley Exp $
!$==============================================================================

!^########################## REJECT_INI ########################################
subroutine reject_ini(num_param, auto_acc, auto_del, nsigs_acc, nsigs_del, &
   nsigs_fraction, time_ref, verbos, rej_opp, errcod)

integer,        intent(in) :: num_param
logical,        intent(in) :: auto_acc
logical,        intent(in) :: auto_del
real (kind=wp), intent(in) :: nsigs_acc
real (kind=wp), intent(in) :: nsigs_del
real (kind=wp), intent(in) :: nsigs_fraction
real (kind=wp), intent(in) :: time_ref
logical,        intent(in) :: verbos
logical,        intent(in) :: rej_opp
integer,        intent(out) :: errcod

!===============================================================================
! reject_ini loads a few parameters into the module and allocates some of 
!    the arrays that will be used later.
!
! Inputs:
!    num_param - the number of estimated parameters, including the elements. 
!       This is used only to decide whether or not we can afford to reduce 
!       the number of available observations.
!    auto_acc - auto_accept flag
!    auto_del - auto_delete flag
!    nsigs_acc - number of sigmas at which an observation is recovered
!    nsigs_del - number of sigmas at which an observation is rejected
!    nsigs_fraction - the rejection threshold is set to the maximum of 
!       nsigs_del and nsigs_fraction*maximum_included_sigma
!    time_ref - reference epoch (JD)
!    rej_opp - determines whether the last observation in an apparition 
!       may be rejected.
!    verbos - flag to control amount of output
!
! Outputs:
!    errcod - Zero indicates success.
!
!$==============================================================================

!^########################## REJECT ############################################
subroutine reject(src_matrix, change_flag, errcod, od_update)

real (kind=wp),    intent(in)  :: src_matrix(:)
logical,           intent(out) :: change_flag
integer,           intent(out) :: errcod
logical, optional, intent(in)  :: od_update

!===============================================================================
! reject rejects and recovers observational outliers.
!
! Inputs:
!    src_matrix - the sqrt covariance matrix of the current solution
!
! Output:
!    change_flag - Indicates if there was a change in the rejection flags.
!    errcod - Zero indicates success. Negative indicates success, 
!       but with some anomaly.
!
! Optional Input:
!    od_update - Logical flag to indicate whether the od_mod data should be 
!       updated with the results of the rejection pass. This should be set 
!       to false if there will be no further iterations before the obs file 
!       is written back to file.
!
!$==============================================================================

>>>>>>>>>>>>>>>>>>> FILE: ../od/lov_mod.f90 <<<<<<<<<<<<<<<<<<<

module lov_mod

use constants_mod, only: wp, deg_per_rad, sin_eps, cos_eps, gm_sun, two_pi, pi
use od_mod, only: od_solve, od_msg, od_constraint
use od_converge_mod, only: od_converge, converge_msg
use matvec_util_mod, only: mat_inv, vec_unit
use estpar_mod, only: n_param, n_utm, n_utm_rhs

implicit none
save
private ! by default

! Public Procedures:
public :: lov_ini, lov_clear, lov_step, lov_correct, lov_rhs, lov_reset_weakness

! Public module variables:
public :: lov_msg

! Retrievable error message
character (len=180) :: lov_msg

!===============================================================================
! This is the lov_mod. Really.
!
!-------------------------------------------------------------------------------
! Written       Aug 14, 2001         Steve Chesley
! Modified      Oct 12, 2001         Steve Chesley
!    - Major overhaul.
! Modified      Nov 14, 2001         Steve Chesley
!    - Added lov_clear routine
!    - Ensured src_matrix outputs are always dimensioned at n_utm.
! Modified      Dec 18, 2001         Steve Chesley
!    - Tweaked lov_ini to make lov_clear unnecessary, but left lov_clear alone.
! Modified      Dec 27, 2001         Steve Chesley
!    - Tweaked lov_ini to force updating of reference direction to nominal 
!      at every call.
!    - Added lov_reset_weakness to allow updating of reference direction to 
!      nominal without call to lov_ini
! Modified      Mar 06, 2002       Steve Chesley
!   - Added "only:" restrictions to use statements.
!   - Eliminated p_name IAW new F90 propagator routines.
! Modified      May 13, 2002       Steve Chesley
!   - Use estpar_mod for n_param, etc.
! Modified      Feb 19, 2004       Steve Chesley
!   - Tweaked scaling rules in lov_rhs.
! Modified      Oct 18, 2005       Steve Chesley
!   - Allowed WEAK_FRAME=SMA for LOV in semimajor axis direction
!
! $Id: lov_mod.f90,v 1.7 2007/06/14 21:39:31 chesley Exp $
!$==============================================================================

!^########################## LOV_INI ###########################################
subroutine lov_ini(method, weak_frame, rescale, epoch_jed, ref_elements, &
   weakness, errcod, nom_weakness)

character (len=3),        intent(in)  :: method
character (len=3),        intent(in)  :: weak_frame
logical,                  intent(in)  :: rescale
real (kind=wp),           intent(in)  :: epoch_jed
real (kind=wp),           intent(in)  :: ref_elements(:)
real (kind=wp),           intent(out) :: weakness(:)
integer,                  intent(out) :: errcod
real (kind=wp), optional, intent(in)  :: nom_weakness(:)

!===============================================================================
! lov_ini initializes several module variables based on the reference orbit.
!
! Input:
!   method       - LOV integration method, e.g., "EUL" or "RK2"
!   weak_frame   - Reference frame to be used in the LOV integration
!   rescale      - Use scaled variables in the LOV integration frame
!   epoch_jed    - JED epoch of elements
!   ref_elements - Reference element set
!
! Output:
!   weakness     - one-sigma weak direction vector
!   errcod       - error code. Zero indicates success.
!
! Optional Input:
!   nom_weakness - Externally derived nominal weakness, which will be saved
!      as the reference direction instead of the local nominal weakness.
!
!$==============================================================================

!^########################## LOV_RESET_WEAKNESS ################################
subroutine lov_reset_weakness(errcod)
integer,        intent(out) :: errcod

!===============================================================================
! lov_reset_weakness loads the original nominal weak direction back into the 
!   lov_mod reference weak direction.
!
! Output:
!   errcod       - error code. Zero indicates success.
!
!$==============================================================================

!^########################## LOV_CLEAR #########################################
subroutine lov_clear(errcod)

integer,                  intent(out) :: errcod

!===============================================================================
! lov_clear clears the LOV module for another independent analysis.
!
! Output:
!   errcod       - error code. Zero indicates success.
!
!$==============================================================================

!^########################## LOV_CORRECT #######################################
subroutine lov_correct(elements, weakness, src_matrix, rms_resid, errcod)

real (kind=wp), intent(inout) :: elements(:)
real (kind=wp), intent(inout) :: weakness(:)
real (kind=wp), intent(out)   :: src_matrix(:)
real (kind=wp), intent(out)   :: rms_resid
integer,        intent(out)   :: errcod

!===============================================================================
! lov_correct
!
! Input/Output:
!   elements - 
!   weakness   - one-sigma weak eigenvector
! Output:
!   src_matrix - Vector-stored Sqrt Covariance
!   rms_resid  - RMS of residuals
!   errcod     - error code. Zero indicates success.
!$==============================================================================

!^########################## LOV_STEP ##########################################
subroutine lov_step(delta_sig, elements, weakness, src_matrix, errcod)

real (kind=wp), intent(in)    :: delta_sig
real (kind=wp), intent(inout) :: elements(:)
real (kind=wp), intent(inout) :: weakness(:)
real (kind=wp), intent(out)   :: src_matrix(:)
integer,        intent(out)   :: errcod

!===============================================================================
! lov_step
!
! Input:
!   delta_sig  - desired step along LOV
! Input/Output:
!   elements   - starting point and ending point
!   weakness   - corresponding one-sigma weak direction eigenvector
! Output:
!   errcod     - error code. Zero indicates success.
!$==============================================================================

!^########################## LOV_RHS ###########################################
subroutine lov_rhs(elements, weakness, errcod, src_matrix)

real (kind=wp),           intent(in)  :: elements(:)
real (kind=wp),           intent(out) :: weakness(:)
integer,                  intent(out) :: errcod
real (kind=wp), optional, intent(out) :: src_matrix(:)

!===============================================================================
! lov_rhs computes the weak direction times weakness from a given orbit.
!
! Input:
!   elements   - 
! Output:
!   weakness   - one-sigma vector in the weak direction
!   errcod     - Zero indicates success.
! Optional Output:
!   src_matrix - 
!
!$==============================================================================
