[oorb] r275 committed - Tentative versions of autoMCMCRanging4 and MCMCRanging4 (random-walk r...
2 views
Skip to first unread message
oo...@googlecode.com
Mar 10, 2015, 3:58:45 AM3/10/15
to oo...@googlegroups.com
Revision: 275
Author: jxvir...@gmail.com
Date: Tue Mar 10 07:58:16 2015 UTC
Log: Tentative versions of autoMCMCRanging4 and MCMCRanging4
(random-walk ranging)
https://code.google.com/p/oorb/source/detail?r=275
Modified:
/trunk/main/oorb.f90
=======================================
--- /trunk/main/oorb.f90 Thu Feb 26 16:16:35 2015 UTC
+++ /trunk/main/oorb.f90 Tue Mar 10 07:58:16 2015 UTC
@@ -277,7 +277,7 @@
solar_alt, solar_alt_max, solelon_max, solelon_min,
generat_multiplier, stdev, &
step, sun_moon_r2, sunlat, sunlon, &
timespan, tlat, tlon, toclat, toclon, tsclat, tsclon, &
- ta_s, ta_c, fak, ecc_anom, true_anom
+ ta_s, ta_c, fak, ecc_anom, true_anom, chi2, chi2_min
INTEGER, DIMENSION(:), POINTER :: &
repetition_arr_cmp, &
repetition_arr_in
@@ -2278,18 +2278,27 @@
STOP
END IF
IF (info_verb >= 2) THEN
- WRITE(stdout,"(4X,A,1X,F5.2,1X,A)") &
+ WRITE(stdout,"(4X,A,1X,F6.2,1X,A)") &
"Observational arc: ", dt, "days"
END IF
- IF (sor_rho_init(3) > HUGE(sor_rho_init(3))/2) THEN
- ! Initialize the rho2-rho1 range based on the observational
timespan.
- IF (dt > 10) THEN
- sor_rho_init(3) = -0.5_bp
- ELSE
- sor_rho_init(3) = -0.05_bp*dt
- END IF
- sor_rho_init(4) = -sor_rho_init(3)
+!Below not valid for mcmc?
+!!$ IF (sor_rho_init(3) > HUGE(sor_rho_init(3))/2) THEN
+!!$ ! Initialize the rho2-rho1 range based on the observational
timespan.
+!!$ IF (dt > 10) THEN
+!!$ sor_rho_init(3) = -0.5_bp
+!!$ ELSE
+!!$ sor_rho_init(3) = -0.05_bp*dt
+!!$ END IF
+!!$ sor_rho_init(4) = -sor_rho_init(3)
+!!$ END IF
+
+ CALL NEW(storb, obss_sep(i))
+ IF (error) THEN
+ CALL errorMessage("oorb / mcmcranging", &
+ "TRACE BACK (65)", 1)
+ STOP
END IF
+
IF (ASSOCIATED(id_arr_storb_in) .AND. ALLOCATED(storb_arr_in)) THEN
DO j=1,SIZE(id_arr_storb_in)
IF (TRIM(id) == TRIM(id_arr_storb_in(j))) THEN
@@ -2328,6 +2337,10 @@
"TRACE BACK (33)", 1)
STOP
END IF
+ !ML orbit
+ ref_orb=getNominalOrbit(storb_arr_in(j), ml_orbit=.true.)
+ call setParameters(storb, orb_ml=ref_orb)
+ sor_iterate_bounds(:) = .false.
END IF
END IF
IF (.NOT.exist(epoch)) THEN
@@ -2366,12 +2379,6 @@
ELSE
CALL NULLIFY(t)
t = copy(epoch)
- END IF
- CALL NEW(storb, obss_sep(i))
- IF (error) THEN
- CALL errorMessage("oorb / mcmcranging", &
- "TRACE BACK (65)", 1)
- STOP
END IF
CALL setParameters(storb, &
dyn_model=dyn_model, &
@@ -2410,7 +2417,7 @@
END IF
SELECT CASE (sor_type_prm)
CASE (1)
- CALL MCMCRanging2(storb)
+ CALL MCMCRanging3(storb) ! MCMCRanging3 = "basic" MCMC
CASE (2)
!!$ CALL setParameters(storb, &
!!$ dyn_model_init=dyn_model_init, &
@@ -2420,7 +2427,8 @@
!!$ accept_multiplier=accwin_multiplier, &
!!$ apriori_rho_min=apriori_rho_min, &
!!$ set_acceptance_window=.FALSE.)
- CALL autoMCMCRanging3(storb)
+! CALL autoMCMCRanging3(storb)
+ CALL autoMCMCRanging4(storb)
CASE default
CALL errorMessage("oorb / mcmcranging", &
"Unknown type of ranging:", 1)
@@ -2661,7 +2669,8 @@
END IF
DO j=1,6
CALL moments(elements_arr(:,j), &
- pdf=pdf_arr_cmp, std_dev=stdev, errstr=errstr)
+! pdf=pdf_arr_cmp, std_dev=stdev, errstr=errstr)
+ pdf=repetition_arr_cmp, std_dev=stdev, errstr=errstr)
WRITE(getUnit(tmp_file), "(F22.15,1X)", &
advance="no") stdev
stdev_arr(j)=stdev
@@ -2731,6 +2740,605 @@
END IF
END DO
+CASE ("random-walk_ranging")
+
+ ! Orbital inversion using random-walk ranging (Muinonen et al. 2015),
+ ! that is, unlike in MCMCRanging, the acceptance value is based not
+ ! on the quotient of proposal p.d.f.:s, but on the comparison of
+ ! difference of proposal and previous chi^2 values to the predefined
+ ! boundary value, usually 20.1 (corresponding to 3*sigma).
+
+ ! To be done:
+ ! 1. Decide where to put the ran_walk flag (currently in
StOrb/MCMCRanging4)
+ ! and put it there.
+ ! 2. Implement the "take-over" automatisation. If user-given
initial parameters
+ ! fail to produce a decent orbit, assume own initial guesses for
range values.
+ ! These values are 2.5 +/- 0.2 au (MBOs) or 2.8 +/- 2.8 au (also
to include NEOs
+ ! and Jupiter Trojans) and 43 +/- 8 au (TNOs).
+
+ CALL NULLIFY(epoch)
+ dyn_model = " "
+ integrator = " "
+ integration_step = -1.0_bp
+ apriori_a_max = -1.0_bp
+ apriori_a_min = -1.0_bp
+ apriori_periapsis_max = -1.0_bp
+ apriori_periapsis_min = -1.0_bp
+ apriori_apoapsis_max = -1.0_bp
+ apriori_apoapsis_min = -1.0_bp
+ apriori_rho_min = -1.0_bp
+ sor_type_prm = -1
+ sor_2point_method = " "
+ sor_2point_method_sw = " "
+ sor_norb = -1
+ sor_norb_sw = -1
+ sor_ntrial = -1
+ sor_ntrial_sw = -1
+ sor_niter = -1
+ sor_rho_init = HUGE(sor_rho_init)
+ generat_multiplier = -1.0_bp
+ sor_genwin_offset = -1.0_bp
+ sor_iterate_bounds = .TRUE.
+ accwin_multiplier = -1.0_bp
+ gaussian_rho = .FALSE.
+ regularized = .FALSE.
+ write_residuals = .FALSE.
+ CALL readConfigurationFile(conf_file, &
+ t0=epoch, &
+ dyn_model=dyn_model, &
+ perturbers=perturbers, &
+ integrator=integrator, &
+ integration_step=integration_step, &
+ dyn_model_init=dyn_model_init, &
+ integrator_init=integrator_init, &
+ integration_step_init=integration_step_init, &
+ apriori_a_max=apriori_a_max, &
+ apriori_a_min=apriori_a_min, &
+ apriori_periapsis_max=apriori_periapsis_max, &
+ apriori_periapsis_min=apriori_periapsis_min, &
+ apriori_apoapsis_max=apriori_apoapsis_max, &
+ apriori_apoapsis_min=apriori_apoapsis_min, &
+ apriori_rho_min=apriori_rho_min, &
+ sor_type=sor_type_prm, sor_2point_method=sor_2point_method, &
+ sor_2point_method_sw=sor_2point_method_sw, &
+ sor_norb=sor_norb, sor_norb_sw=sor_norb_sw, &
+ sor_ntrial=sor_ntrial, sor_ntrial_sw=sor_ntrial_sw, &
+ sor_niter=sor_niter, sor_rho_init=sor_rho_init, &
+ generat_multiplier=generat_multiplier, &
+ sor_genwin_offset=sor_genwin_offset, &
+ sor_iterate_bounds=sor_iterate_bounds, &
+ accwin_multiplier=accwin_multiplier, &
+ regularized_pdf=regularized, &
+ sor_random_obs=random_obs, sor_rho_gauss=gaussian_rho, &
+ write_residuals=write_residuals, &
+ ls_correction_factor=ls_correction_factor, &
+ ls_element_mask=ls_element_mask, &
+ ls_niter_major_max=ls_niter_major_max, &
+ ls_niter_major_min=ls_niter_major_min, &
+ ls_niter_minor=ls_niter_minor)
+ IF (error) THEN
+ CALL errorMessage("oorb / random-walk_ranging", &
+ "TRACE BACK (5)", 1)
+ STOP
+ END IF
+ IF (get_cl_option("--epoch-mjd-tt=", .FALSE.)) THEN
+ CALL NULLIFY(epoch)
+ ! New epoch given as MJD TT
+ mjd_tt = get_cl_option("--epoch-mjd-tt=", 0.0_bp)
+ CALL NEW(epoch, mjd_tt, "TT")
+ IF (error) THEN
+ CALL errorMessage("oorb / random-walk_ranging", &
+ "TRACE BACK (10)", 1)
+ STOP
+ END IF
+ END IF
+
+ ! Print header before printing first orbit:
+ first = .TRUE.
+ DO i=1,SIZE(obss_sep,dim=1)
+ id = getID(obss_sep(i))
+ IF (error) THEN
+ CALL errorMessage("oorb / random-walk_ranging", &
+ "TRACE BACK (15)", 1)
+ STOP
+ END IF
+ IF (info_verb >= 2) THEN
+ WRITE(stdout,"(1X,I0,3(A),1X,I0,A,I0)") i, &
+ ". observation set (", TRIM(id), ")."
+ END IF
+ nobs = getNrOfObservations(obss_sep(i))
+ IF (error) THEN
+ CALL errorMessage("oorb / random-walk_ranging", &
+ "TRACE BACK (20)", 1)
+ STOP
+ END IF
+ dt = getObservationalTimespan(obss_sep(i))
+ IF (error) THEN
+ CALL errorMessage("oorb / random-walk_ranging", &
+ "TRACE BACK (40)", 1)
+ STOP
+ END IF
+ IF (info_verb >= 2) THEN
+ WRITE(stdout,"(4X,A,1X,F5.2,1X,A)") &
+ "Observational arc: ", dt, "days"
+ END IF
+ IF (sor_rho_init(3) > HUGE(sor_rho_init(3))/2) THEN
+ ! Initialize the rho2-rho1 range based on the observational
timespan.
+ IF (dt > 10) THEN
+ sor_rho_init(3) = -0.5_bp
+ ELSE
+ sor_rho_init(3) = -0.05_bp*dt
+ END IF
+ sor_rho_init(4) = -sor_rho_init(3)
+ END IF
+ !chi2_min assuming the error model is perfect
+ obs_masks => getObservationMasks(obss_sep(i))
+ chi2_min_init=REAL(COUNT(obs_masks), bp)
+ print*,"Chi2min_init ", chi2_min_init
+ IF (ASSOCIATED(id_arr_storb_in) .AND. ALLOCATED(storb_arr_in)) THEN
+ DO j=1,SIZE(id_arr_storb_in)
+ IF (TRIM(id) == TRIM(id_arr_storb_in(j))) THEN
+ EXIT
+ END IF
+ END DO
+ IF (j <= SIZE(id_arr_storb_in)) THEN
+ CALL setParameters(storb_arr_in(j), &
+ dyn_model=dyn_model_init, &
+ perturbers=perturbers, &
+ integrator=integrator_init, &
+ integration_step=integration_step_init, &
+ accept_multiplier=accwin_multiplier, &
+ apriori_rho_min=apriori_rho_min, &
+ set_acceptance_window=.FALSE.)
+ IF (error) THEN
+ CALL errorMessage("oorb / random-walk_ranging", &
+ "TRACE BACK (30)", 1)
+ STOP
+ END IF
+ CALL constrainRangeDistributions(storb_arr_in(j),
obss_sep(i))
+ IF (error) THEN
+ CALL errorMessage("oorb / random-walk_ranging", &
+ "TRACE BACK (31)", 1)
+ STOP
+ END IF
+ ! Allow user to decide whether rho bounds are updated from
orbit file
+ ! (but always update chi2min if orbit file is present)
+ IF (ALL(sor_iterate_bounds)) then
+ CALL setRangeBounds(storb_arr_in(j))
+ IF (error) THEN
+ CALL errorMessage("oorb / random-walk_ranging", &
+ "TRACE BACK (32)", 1)
+ STOP
+ END IF
+
+ sor_rho_init(1:4) = getRangeBounds(storb_arr_in(j))
+ IF (error) THEN
+ CALL errorMessage("oorb / random-walk_ranging", &
+ "TRACE BACK (33)", 1)
+ STOP
+ END IF
+ end if
+ ! Initialize chi2min
+ orb_arr_cmp => getSampleOrbits(storb_arr_in(j))
+ IF (error) THEN
+ CALL errorMessage("oorb / random-walk_ranging", &
+ "TRACE BACK (34)", 1)
+ STOP
+ END IF
+ chi2_min=HUGE(chi2_min)
+ DO k=1,SIZE(orb_arr_cmp,dim=1)
+ chi2=getChi2(storb_arr_in(j), orb_arr_cmp(k))
+ if (chi2 < chi2_min_init) chi2_min_init=chi2
+ END DO
+ END IF
+ END IF
+! print*,chi2_min_init
+ IF (.NOT.exist(epoch)) THEN
+ CALL NULLIFY(t)
+ obs = getObservation(obss_sep(i),1)
+ IF (error) THEN
+ CALL errorMessage("oorb / random-walk_ranging", &
+ "TRACE BACK (45)", 1)
+ STOP
+ END IF
+ t = getTime(obs)
+ IF (error) THEN
+ CALL errorMessage("oorb / random-walk_ranging", &
+ "TRACE BACK (50)", 1)
+ STOP
+ END IF
+ CALL NULLIFY(obs)
+ mjd = getMJD(t, "tt")
+ IF (error) THEN
+ CALL errorMessage("oorb / random-walk_ranging", &
+ "TRACE BACK (55)", 1)
+ STOP
+ END IF
+ CALL NULLIFY(t)
+ IF (get_cl_option("--exact-mid-epoch", .FALSE.)) THEN
+ mjd = mjd+dt/2.0_bp
+ ELSE
+ mjd = REAL(NINT(mjd+dt/2.0_bp),bp)
+ END IF
+ CALL NEW(t, mjd, "tt")
+ IF (error) THEN
+ CALL errorMessage("oorb / random-walk_ranging", &
+ "TRACE BACK (60)", 1)
+ STOP
+ END IF
+ ELSE
+ CALL NULLIFY(t)
+ t = copy(epoch)
+ END IF
+ CALL NEW(storb, obss_sep(i))
+ IF (error) THEN
+ CALL errorMessage("oorb / random-walk_ranging", &
+ "TRACE BACK (65)", 1)
+ STOP
+ END IF
+ CALL setParameters(storb, &
+ dyn_model=dyn_model, &
+ perturbers=perturbers, &
+ integrator=integrator, &
+ integration_step=integration_step, &
+ outlier_rejection=outlier_rejection_prm, &
+ outlier_multiplier=outlier_multiplier_prm, &
+ t_inv=t, &
+ element_type=element_type_comp_prm, &
+ dchi2_rejection = dchi2_rejection, &
+ dchi2_max = dchi2_max, &
+ regularized_pdf = regularized, &
+ accept_multiplier=accwin_multiplier, &
+ apriori_a_max=apriori_a_max, apriori_a_min=apriori_a_min, &
+ apriori_periapsis_max=apriori_periapsis_max, &
+ apriori_periapsis_min=apriori_periapsis_min, &
+ apriori_apoapsis_max=apriori_apoapsis_max, &
+ apriori_apoapsis_min=apriori_apoapsis_min, &
+ apriori_rho_min=apriori_rho_min, &
+ sor_2point_method=sor_2point_method, &
+ sor_2point_method_sw=sor_2point_method_sw, &
+ sor_norb_sw=sor_norb_sw, sor_ntrial_sw=sor_ntrial_sw, &
+ sor_norb=sor_norb, sor_ntrial=sor_ntrial, &
+ sor_rho1_l=sor_rho_init(1), sor_rho1_u=sor_rho_init(2), &
+ sor_rho2_l=sor_rho_init(3), sor_rho2_u=sor_rho_init(4), &
+ sor_niter=sor_niter, sor_iterate_bounds=sor_iterate_bounds, &
+ sor_random_obs_selection=.FALSE., &
+ gaussian_pdf=gaussian_rho, &
+ generat_multiplier=generat_multiplier, &
+ sor_generat_offset=sor_genwin_offset, &
+ chi2_min = chi2_min_init)
+ IF (error) THEN
+ CALL errorMessage("oorb / random-walk_ranging", &
+ "TRACE BACK (70)", 1)
+ STOP
+ END IF
+ SELECT CASE (sor_type_prm)
+ CASE (1)
+ CALL MCMCRanging4(storb) ! MCMCRanging3 = "basic" MCMC
+ CASE (2)
+!!$ CALL setParameters(storb, &
+!!$ dyn_model_init=dyn_model_init, &
+!!$ perturbers=perturbers, &
+!!$ integrator_init=integrator_init, &
+!!$ integration_step_init=integration_step_init, &
+!!$ accept_multiplier=accwin_multiplier, &
+!!$ apriori_rho_min=apriori_rho_min, &
+!!$ set_acceptance_window=.FALSE.)
+ !CALL autoMCMCRanging3(storb)
+ CALL errorMessage("oorb / random-walk_ranging", &
+ "Automatic ranging not implemented for
random-walk_ranging:", 1)
+ CASE default
+ CALL errorMessage("oorb / random-walk_ranging", &
+ "Unknown type of ranging:", 1)
+ IF (err_verb >= 1) THEN
+ WRITE(stderr,"(2X,I0)") sor_type_prm
+ END IF
+ STOP
+ END SELECT
+ IF (error) THEN
+ CALL errorMessage("oorb / random-walk_ranging", &
+ "MCMC failed:", 1)
+ IF (err_verb >= 1) THEN
+ WRITE(stderr,"(3A,1X,I0)") "ID: ", TRIM(id), &
+ " and number of observations: ", nobs
+ END IF
+ error = .FALSE.
+ CALL NEW(out_file, "problematic_observation_sets" // "." // &
+ TRIM(observation_format_out))
+ CALL setPositionAppend(out_file)
+ CALL OPEN(out_file)
+ IF (error) THEN
+ CALL errorMessage("oorb / random-walk_ranging", &
+ "TRACE BACK (90)", 1)
+ STOP
+ END IF
+ CALL writeObservationFile(obss_sep(i), getUnit(out_file), &
+ TRIM(observation_format_out))
+ IF (error) THEN
+ CALL errorMessage("oorb / random-walk_ranging", &
+ "TRACE BACK (95)", 1)
+ STOP
+ END IF
+ CALL NULLIFY(out_file)
+ IF (error) THEN
+ CALL errorMessage("oorb / random-walk_ranging", &
+ "TRACE BACK (100)", 1)
+ STOP
+ END IF
+ ELSE
+ IF (info_verb >= 2) THEN
+ WRITE(stdout,"(3(1X,A))") "Mcmc inversion for object", &
+ TRIM(id), "is ready."
+ END IF
+ CALL NEW(out_file, TRIM(id) // ".sor")
+ CALL OPEN(out_file)
+ IF (error) THEN
+ CALL errorMessage("oorb / random-walk_ranging", &
+ "TRACE BACK (105)", 1)
+ STOP
+ END IF
+ ! WRITE OBSERVATIONS:
+ CALL writeObservationFile(obss_sep(i), getUnit(out_file), &
+ TRIM(observation_format_out))
+ IF (error) THEN
+ CALL errorMessage("oorb / random-walk_ranging", &
+ "TRACE BACK (110)", 1)
+ STOP
+ END IF
+ WRITE(getUnit(out_file),"(A)") "#"
+ ! WRITE MCMC PARAMETERS
+ CALL writeSORResults(storb, obss_sep(i), getUnit(out_file))
+ IF (error) THEN
+ CALL errorMessage("oorb / mrandom-walk_ranging", &
+ "TRACE BACK (115)", 1)
+ STOP
+ END IF
+ CALL NULLIFY(out_file)
+ ! WRITE ORBITAL-ELEMENT PDF
+ orb_arr_cmp => getSampleOrbits(storb)
+ IF (error) THEN
+ CALL errorMessage("oorb / random-walk_ranging", &
+ "TRACE BACK (125)", 1)
+ STOP
+ END IF
+ pdf_arr_cmp => getPDFValues(storb, element_type_comp_prm)
+ IF (error) THEN
+ CALL errorMessage("oorb / random-walk_ranging", &
+ "TRACE BACK (130)", 1)
+ STOP
+ END IF
+ rchi2_arr_cmp => getReducedChi2Distribution(storb)
+ IF (error) THEN
+ CALL errorMessage("oorb / random-walk_ranging", &
+ "TRACE BACK (135)", 1)
+ STOP
+ END IF
+ CALL getResults(storb, repetition_arr_cmp=repetition_arr_cmp)
+!!$ CALL getResults(storb, &
+!!$ reg_apr_arr=reg_apr_arr_cmp, &
+!!$ jac_arr=jac_arr_cmp)
+!!$ IF (error) THEN
+!!$ CALL errorMessage("oorb / mcmcranging ", &
+!!$ "TRACE BACK (140)", 1)
+!!$ STOP
+!!$ END IF
+ IF (separately) THEN
+ CALL NEW(out_file, TRIM(id) // ".orb")
+ CALL OPEN(out_file)
+ IF (error) THEN
+ CALL errorMessage("oorb / random-walk_ranging", &
+ "TRACE BACK (145)", 1)
+ STOP
+ END IF
+ lu_orb_out = getUnit(out_file)
+ END IF
+ DO j=1,SIZE(orb_arr_cmp,dim=1)
+ IF (orbit_format_out == "orb") THEN
+ CALL writeOpenOrbOrbitFile(lu_orb_out, &
+ print_header=j==1, &
+ element_type_out=element_type_out_prm, &
+ id=id, &
+ orb=orb_arr_cmp(j), &
+ element_type_pdf=element_type_comp_prm, &
+ pdf=pdf_arr_cmp(j), &
+ rchi2=rchi2_arr_cmp(j), &
+ repetitions=repetition_arr_cmp(j), &
+ mjd=mjd_epoch)
+!!$ reg_apr=reg_apr_arr_cmp(j), &
+!!$ jac_sph_inv=jac_arr_cmp(j,1), &
+!!$ jac_car_kep=jac_arr_cmp(j,2), &
+!!$ jac_equ_kep=jac_arr_cmp(j,3))
+ ELSE IF (orbit_format_out == "des") THEN
+ CALL errorMessage("oorb / random-walk_rangingg", &
+ "DES format not yet supported for MCMC output.", 1)
+ STOP
+ END IF
+ IF (error) THEN
+ CALL errorMessage("oorb / random-walk_ranging", &
+ "TRACE BACK (150)", 1)
+ STOP
+ END IF
+ END DO
+ IF (separately) THEN
+ CALL NULLIFY(out_file)
+ IF (compress) THEN
+ CALL system("gzip -f " // TRIM(id) // ".orb")
+ END IF
+ END IF
+ ! WRITE RESIDUALS
+ IF (write_residuals) THEN
+ CALL NEW(out_file, TRIM(out_fname) // ".res")
+ CALL setPositionAppend(out_file)
+ CALL OPEN(out_file)
+ IF (error) THEN
+ CALL errorMessage("oorb / random-walk_ranging", &
+ "TRACE BACK (155)", 1)
+ STOP
+ END IF
+ CALL writeResiduals(storb, obss_sep(i), getUnit(out_file))
+ IF (error) THEN
+ CALL errorMessage("oorb / random-walk_ranging", &
+ "TRACE BACK (160)", 1)
+ STOP
+ END IF
+ CALL NULLIFY(out_file)
+ END IF
+ IF (plot_results) THEN
+ CALL toString(dt, str, error, frmt="(F10.2)")
+ IF (error) THEN
+ CALL errorMessage("oorb / random-walk_ranging", &
+ "TRACE BACK (165)", 1)
+ STOP
+ END IF
+ str = TRIM(id) // "_"// TRIM(str)
+ CALL makeResidualStamps(storb, obss_sep(i), TRIM(str) // &
+ "_random-walk_ranging_residual_stamps.eps")
+ IF (error) THEN
+ CALL errorMessage("oorb / random-walk_ranging", &
+ "TRACE BACK (170)", 1)
+ STOP
+ END IF
+ IF (compress) THEN
+ CALL system("gzip -f " // TRIM(str)
// "_random-walk_ranging_residual_stamps.eps")
+ END IF
+ IF (plot_open) THEN
+ CALL system("gv " // TRIM(str)
// "_random-walk_ranging_residual_stamps.eps* &")
+ END IF
+ ALLOCATE(elements_arr(SIZE(orb_arr_cmp,dim=1),7), stat=err)
+ IF (err /= 0) THEN
+ CALL errorMessage("oorb / random-walk_ranging", &
+ "Could not allocate memory (3).", 1)
+ STOP
+ END IF
+ CALL NEW(tmp_file,
TRIM(str)// "_random-walk_ranging_orbits.out")
+ CALL OPEN(tmp_file)
+ IF (error) THEN
+ CALL errorMessage("oorb / random-walk_ranging ", &
+ "TRACE BACK (175)", 1)
+ STOP
+ END IF
+ DO j=1,SIZE(orb_arr_cmp,dim=1)
+ IF (element_type_comp_prm == "cartesian") THEN
+ CALL rotateToEcliptic(orb_arr_cmp(j))
+ END IF
+ elements_arr(j,1:6) = getElements(orb_arr_cmp(j),
element_type_comp_prm)
+ IF (error) THEN
+ CALL errorMessage("oorb / random-walk_ranging", &
+ "TRACE BACK (180)", 1)
+ STOP
+ END IF
+ IF (element_type_comp_prm == "keplerian") THEN
+ elements_arr(j,3:6) = elements_arr(j,3:6)/rad_deg
+ END IF
+ elements_arr(j,7) = pdf_arr_cmp(j)
+ t = getTime(orb_arr_cmp(j))
+ IF (error) THEN
+ CALL errorMessage("oorb / random-walk_ranging ", &
+ "TRACE BACK (185)", 1)
+ STOP
+ END IF
+ WRITE(getUnit(tmp_file),*) &
+ elements_arr(j,1:6), &
+ pdf_arr_cmp(j), &
+ getCalendarDateString(t,"tdt")
+ IF (error) THEN
+ CALL errorMessage("oorb / random-walk_ranging ", &
+ "TRACE BACK (190)", 1)
+ STOP
+ END IF
+ CALL NULLIFY(t)
+ END DO
+ CALL NULLIFY(tmp_file)
+ CALL NEW(tmp_file, TRIM(str) // &
+ "_random-walk_ranging_sample_standard_deviations.out")
+ CALL OPEN(tmp_file)
+ IF (error) THEN
+ CALL errorMessage("oorb / random-walk_ranging ", &
+ "TRACE BACK (195)", 1)
+ STOP
+ END IF
+ WRITE(getUnit(tmp_file), "(F22.15,1X)", &
+ advance="no") &
+ getObservationalTimespan(obss_sep(i))
+ IF (error) THEN
+ CALL errorMessage("oorb / random-walk_ranging ", &
+ "TRACE BACK (200)", 1)
+ STOP
+ END IF
+ DO j=1,6
+ CALL moments(elements_arr(:,j), &
+ pdf=pdf_arr_cmp, std_dev=stdev, errstr=errstr)
+ WRITE(getUnit(tmp_file), "(F22.15,1X)", &
+ advance="no") stdev
+ stdev_arr(j)=stdev
+ END DO
+ WRITE(getUnit(tmp_file),*)
+ CALL NULLIFY(tmp_file)
+ DEALLOCATE(elements_arr, stat=err)
+ IF (err /= 0) THEN
+ CALL errorMessage("oorb / random-walk_ranging", &
+ "Could not deallocate memory (5).", 1)
+ STOP
+ END IF
+ ! Make plot using gnuplot:
+ CALL system("cp " // TRIM(str) // &
+ "_random-walk_ranging_orbits.out sor_orbits.out")
+ IF (element_type_comp_prm == "cartesian") THEN
+ CALL system("gnuplot " // TRIM(gnuplot_scripts_dir)
// "/sor_plot_car.gp")
+ ELSE
+ if (stdev_arr(1) < 1.0_bp) then
+ CALL system("gnuplot " // TRIM(gnuplot_scripts_dir)
// "/sor_plot_kep_nolog.gp")
+ else
+ CALL system("gnuplot " // TRIM(gnuplot_scripts_dir)
// "/sor_plot_kep.gp")
+ end if
+ END IF
+ CALL system("cp sor_results.eps " // TRIM(str) // &
+ "_random-walk_ranging_" // TRIM(element_type_comp_prm)
// &
+ "_results.eps")
+ CALL system("rm -f sor_orbits.out sor_results.eps " // &
+ TRIM(str) // "_random-walk_ranging_orbits.out ")
+!!$ // TRIM(str) // &
+!!$ "_random-walk_ranging_sample_standard_deviations.out")
+ IF (compress) THEN
+ CALL system("gzip -f " // TRIM(str)
// "_random-walk_ranging_" // &
+ TRIM(element_type_comp_prm) // "_results.eps")
+ END IF
+ IF (plot_open) THEN
+ CALL system("gv " // TRIM(str) // "_random-walk_ranging_"
// &
+ TRIM(element_type_comp_prm) // &
+ "_results.eps* &")
+ END IF
+ END IF
+ DO j=1,SIZE(orb_arr_cmp)
+ CALL NULLIFY(orb_arr_cmp(j))
+ END DO
+ DEALLOCATE(orb_arr_cmp, stat=err)
+ DEALLOCATE(pdf_arr_cmp, stat=err)
+ DEALLOCATE(rchi2_arr_cmp, stat=err)
+ DEALLOCATE(repetition_arr_cmp, stat=err)
+!!$ DEALLOCATE(reg_apr_arr_cmp, stat=err)
+!!$ DEALLOCATE(jac_arr_cmp, stat=err)
+ IF (err /= 0) THEN
+ CALL errorMessage("oorb / random-walk_ranging", &
+ "Could not deallocate memory (10).", 1)
+ STOP
+ END IF
+ IF (info_verb >= 2) THEN
+ WRITE(stdout,"(3(1X,A))") "Object", &
+ TRIM(id), "successfully processed."
+ END IF
+
+ END IF
+ CALL NULLIFY(storb)
+ CALL NULLIFY(orb)
+ IF (info_verb > 2) THEN
+ WRITE(stdout,*)
+ WRITE(stdout,*)
+ END IF
+ END DO
+
+
+!*********************
CASE ("simplex")
@@ -6499,25 +7107,25 @@
IF (info_verb >= 2 .AND. dyn_model /= "2-body") THEN
CALL getCalendarDate(epoch0, "TT", year0, month0, day0)
CALL getCalendarDate(epoch, "TT", year1, month1, day1)
- WRITE(stderr,'(A)') ""
- WRITE(stderr,'(A,I0,A,I0,A,F0.5,A,I0,A,I0,A,F0.5,A)') &
- "Planetary encounters by object " //
TRIM(id_arr_storb_in(i)) // &
- " between ", year0, "-", month0, "-", day0, " and ",
&
- year1, "-", month1, "-", day1, ":"
+!~ WRITE(stderr,'(A)') ""
+!~ WRITE(stderr,'(A,I0,A,I0,A,F0.5,A,I0,A,I0,A,F0.5,A)') &
+!~ "Planetary encounters by object " //
TRIM(id_arr_storb_in(i)) // &
+!~ " between ", year0, "-", month0, "-", day0, "
and ", &
+!~ year1, "-", month1, "-", day1, ":"
DO j=1,SIZE(encounters,dim=1)
- WRITE(stderr,'(A,I0,A)') "Orbit #", j, ":"
+!~ WRITE(stderr,'(A,I0,A)') "Orbit #", j, ":"
DO k=1,SIZE(encounters,dim=2)
IF (encounters(j,k,2) < 1.1_bp) THEN
WRITE(stderr,'(A22,1X,A7,1X,A6,1X,F15.6,1X,A,1X,F0.6,1X,A)') &
"!! I-M-P-A-C-T !! with",
planetary_locations(k), &
"at MJD", encounters(j,k,1), "TT given a
stepsize of", &
encounters(j,k,4), "days."
- ELSE
-
WRITE(stderr,'(A22,1X,A7,1X,A6,1X,F15.6,1X,A,1X,F11.8,1X,A,1X,F0.6,1X,A)') &
- "Closest encounter with",
planetary_locations(k), &
- "at MJD", encounters(j,k,1), &
- "TT at a distance of", encounters(j,k,3), &
- "AU given a stepsize of",
encounters(j,k,4), "days."
+!~ ELSE
+!~
WRITE(stderr,'(A22,1X,A7,1X,A6,1X,F15.6,1X,A,1X,F11.8,1X,A,1X,F0.6,1X,A)') &
+!~ "Closest encounter with",
planetary_locations(k), &
+!~ "at MJD", encounters(j,k,1), &
+!~ "TT at a distance of",
encounters(j,k,3), &
+!~ "AU given a stepsize of",
encounters(j,k,4), "days."
END IF
END DO
END DO
@@ -6799,13 +7407,13 @@
STOP
END IF
IF (info_verb >= 2 .AND. dyn_model /= "2-body") THEN
- CALL getCalendarDate(epoch0, "TT", year0, month0, day0)
- CALL getCalendarDate(epoch, "TT", year1, month1, day1)
- WRITE(stderr,'(A)') ""
-
WRITE(stderr,'(A,I0,A,I0,A,I0,A,F0.5,A,I0,A,I0,A,F0.5,A)') &
- "Planetary encounters by object " //
TRIM(id_arr_in(i)) // &
- " (orbit #", i, ") between ", year0, "-", month0, &
- "-", day0, " and ", year1, "-", month1, "-",
day1, ":"
+!~ CALL getCalendarDate(epoch0, "TT", year0, month0, day0)
+!~ CALL getCalendarDate(epoch, "TT", year1, month1, day1)
+!~ WRITE(stderr,'(A)') ""
+!~
WRITE(stderr,'(A,I0,A,I0,A,I0,A,F0.5,A,I0,A,I0,A,F0.5,A)') &
+!~ "Planetary encounters by object " //
TRIM(id_arr_in(i)) // &
+!~ " (orbit #", i, ") between ", year0, "-", month0,
&
+!~ "-", day0, " and ", year1, "-", month1, "-",
day1, ":"
DO j=1,SIZE(encounters,dim=2)
IF (encounters(1,j,2) < 1.1_bp) THEN
! encounters(1,j,2) == 1 == impact
@@ -6813,12 +7421,11 @@
"!! I-M-P-A-C-T !! with",
planetary_locations(j), &
"at MJD", encounters(1,j,1), "TT given a
stepsize of", &
encounters(1,j,4), "days."
- ELSE
- ! encounters(1,j,2) == 2 == non-impact
-
WRITE(stderr,'(A22,1X,A7,1X,A6,1X,F15.6,1X,A,1X,F11.8,1X,A,1X,F0.6,1X,A)') &
- "Closest encounter with",
planetary_locations(j), &
- "at MJD", encounters(1,j,1), "TT at a distance
of", &
- encounters(1,j,3), "AU given a stepsize of",
encounters(1,j,4), "days."
+!~ ELSE
+!~
WRITE(stderr,'(A22,1X,A7,1X,A6,1X,F15.6,1X,A,1X,F11.8,1X,A,1X,F0.6,1X,A)') &
+!~ "Closest encounter with",
planetary_locations(j), &
+!~ "at MJD", encounters(1,j,1), "TT at a
distance of", &
+!~ encounters(1,j,3), "AU given a stepsize
of", encounters(1,j,4), "days."
END IF
END DO
END IF
Author: jxvir...@gmail.com
Date: Tue Mar 10 07:58:16 2015 UTC
Log: Tentative versions of autoMCMCRanging4 and MCMCRanging4
(random-walk ranging)
https://code.google.com/p/oorb/source/detail?r=275
Modified:
/trunk/main/oorb.f90
=======================================
--- /trunk/main/oorb.f90 Thu Feb 26 16:16:35 2015 UTC
+++ /trunk/main/oorb.f90 Tue Mar 10 07:58:16 2015 UTC
@@ -277,7 +277,7 @@
solar_alt, solar_alt_max, solelon_max, solelon_min,
generat_multiplier, stdev, &
step, sun_moon_r2, sunlat, sunlon, &
timespan, tlat, tlon, toclat, toclon, tsclat, tsclon, &
- ta_s, ta_c, fak, ecc_anom, true_anom
+ ta_s, ta_c, fak, ecc_anom, true_anom, chi2, chi2_min
INTEGER, DIMENSION(:), POINTER :: &
repetition_arr_cmp, &
repetition_arr_in
@@ -2278,18 +2278,27 @@
STOP
END IF
IF (info_verb >= 2) THEN
- WRITE(stdout,"(4X,A,1X,F5.2,1X,A)") &
+ WRITE(stdout,"(4X,A,1X,F6.2,1X,A)") &
"Observational arc: ", dt, "days"
END IF
- IF (sor_rho_init(3) > HUGE(sor_rho_init(3))/2) THEN
- ! Initialize the rho2-rho1 range based on the observational
timespan.
- IF (dt > 10) THEN
- sor_rho_init(3) = -0.5_bp
- ELSE
- sor_rho_init(3) = -0.05_bp*dt
- END IF
- sor_rho_init(4) = -sor_rho_init(3)
+!Below not valid for mcmc?
+!!$ IF (sor_rho_init(3) > HUGE(sor_rho_init(3))/2) THEN
+!!$ ! Initialize the rho2-rho1 range based on the observational
timespan.
+!!$ IF (dt > 10) THEN
+!!$ sor_rho_init(3) = -0.5_bp
+!!$ ELSE
+!!$ sor_rho_init(3) = -0.05_bp*dt
+!!$ END IF
+!!$ sor_rho_init(4) = -sor_rho_init(3)
+!!$ END IF
+
+ CALL NEW(storb, obss_sep(i))
+ IF (error) THEN
+ CALL errorMessage("oorb / mcmcranging", &
+ "TRACE BACK (65)", 1)
+ STOP
END IF
+
IF (ASSOCIATED(id_arr_storb_in) .AND. ALLOCATED(storb_arr_in)) THEN
DO j=1,SIZE(id_arr_storb_in)
IF (TRIM(id) == TRIM(id_arr_storb_in(j))) THEN
@@ -2328,6 +2337,10 @@
"TRACE BACK (33)", 1)
STOP
END IF
+ !ML orbit
+ ref_orb=getNominalOrbit(storb_arr_in(j), ml_orbit=.true.)
+ call setParameters(storb, orb_ml=ref_orb)
+ sor_iterate_bounds(:) = .false.
END IF
END IF
IF (.NOT.exist(epoch)) THEN
@@ -2366,12 +2379,6 @@
ELSE
CALL NULLIFY(t)
t = copy(epoch)
- END IF
- CALL NEW(storb, obss_sep(i))
- IF (error) THEN
- CALL errorMessage("oorb / mcmcranging", &
- "TRACE BACK (65)", 1)
- STOP
END IF
CALL setParameters(storb, &
dyn_model=dyn_model, &
@@ -2410,7 +2417,7 @@
END IF
SELECT CASE (sor_type_prm)
CASE (1)
- CALL MCMCRanging2(storb)
+ CALL MCMCRanging3(storb) ! MCMCRanging3 = "basic" MCMC
CASE (2)
!!$ CALL setParameters(storb, &
!!$ dyn_model_init=dyn_model_init, &
@@ -2420,7 +2427,8 @@
!!$ accept_multiplier=accwin_multiplier, &
!!$ apriori_rho_min=apriori_rho_min, &
!!$ set_acceptance_window=.FALSE.)
- CALL autoMCMCRanging3(storb)
+! CALL autoMCMCRanging3(storb)
+ CALL autoMCMCRanging4(storb)
CASE default
CALL errorMessage("oorb / mcmcranging", &
"Unknown type of ranging:", 1)
@@ -2661,7 +2669,8 @@
END IF
DO j=1,6
CALL moments(elements_arr(:,j), &
- pdf=pdf_arr_cmp, std_dev=stdev, errstr=errstr)
+! pdf=pdf_arr_cmp, std_dev=stdev, errstr=errstr)
+ pdf=repetition_arr_cmp, std_dev=stdev, errstr=errstr)
WRITE(getUnit(tmp_file), "(F22.15,1X)", &
advance="no") stdev
stdev_arr(j)=stdev
@@ -2731,6 +2740,605 @@
END IF
END DO
+CASE ("random-walk_ranging")
+
+ ! Orbital inversion using random-walk ranging (Muinonen et al. 2015),
+ ! that is, unlike in MCMCRanging, the acceptance value is based not
+ ! on the quotient of proposal p.d.f.:s, but on the comparison of
+ ! difference of proposal and previous chi^2 values to the predefined
+ ! boundary value, usually 20.1 (corresponding to 3*sigma).
+
+ ! To be done:
+ ! 1. Decide where to put the ran_walk flag (currently in
StOrb/MCMCRanging4)
+ ! and put it there.
+ ! 2. Implement the "take-over" automatisation. If user-given
initial parameters
+ ! fail to produce a decent orbit, assume own initial guesses for
range values.
+ ! These values are 2.5 +/- 0.2 au (MBOs) or 2.8 +/- 2.8 au (also
to include NEOs
+ ! and Jupiter Trojans) and 43 +/- 8 au (TNOs).
+
+ CALL NULLIFY(epoch)
+ dyn_model = " "
+ integrator = " "
+ integration_step = -1.0_bp
+ apriori_a_max = -1.0_bp
+ apriori_a_min = -1.0_bp
+ apriori_periapsis_max = -1.0_bp
+ apriori_periapsis_min = -1.0_bp
+ apriori_apoapsis_max = -1.0_bp
+ apriori_apoapsis_min = -1.0_bp
+ apriori_rho_min = -1.0_bp
+ sor_type_prm = -1
+ sor_2point_method = " "
+ sor_2point_method_sw = " "
+ sor_norb = -1
+ sor_norb_sw = -1
+ sor_ntrial = -1
+ sor_ntrial_sw = -1
+ sor_niter = -1
+ sor_rho_init = HUGE(sor_rho_init)
+ generat_multiplier = -1.0_bp
+ sor_genwin_offset = -1.0_bp
+ sor_iterate_bounds = .TRUE.
+ accwin_multiplier = -1.0_bp
+ gaussian_rho = .FALSE.
+ regularized = .FALSE.
+ write_residuals = .FALSE.
+ CALL readConfigurationFile(conf_file, &
+ t0=epoch, &
+ dyn_model=dyn_model, &
+ perturbers=perturbers, &
+ integrator=integrator, &
+ integration_step=integration_step, &
+ dyn_model_init=dyn_model_init, &
+ integrator_init=integrator_init, &
+ integration_step_init=integration_step_init, &
+ apriori_a_max=apriori_a_max, &
+ apriori_a_min=apriori_a_min, &
+ apriori_periapsis_max=apriori_periapsis_max, &
+ apriori_periapsis_min=apriori_periapsis_min, &
+ apriori_apoapsis_max=apriori_apoapsis_max, &
+ apriori_apoapsis_min=apriori_apoapsis_min, &
+ apriori_rho_min=apriori_rho_min, &
+ sor_type=sor_type_prm, sor_2point_method=sor_2point_method, &
+ sor_2point_method_sw=sor_2point_method_sw, &
+ sor_norb=sor_norb, sor_norb_sw=sor_norb_sw, &
+ sor_ntrial=sor_ntrial, sor_ntrial_sw=sor_ntrial_sw, &
+ sor_niter=sor_niter, sor_rho_init=sor_rho_init, &
+ generat_multiplier=generat_multiplier, &
+ sor_genwin_offset=sor_genwin_offset, &
+ sor_iterate_bounds=sor_iterate_bounds, &
+ accwin_multiplier=accwin_multiplier, &
+ regularized_pdf=regularized, &
+ sor_random_obs=random_obs, sor_rho_gauss=gaussian_rho, &
+ write_residuals=write_residuals, &
+ ls_correction_factor=ls_correction_factor, &
+ ls_element_mask=ls_element_mask, &
+ ls_niter_major_max=ls_niter_major_max, &
+ ls_niter_major_min=ls_niter_major_min, &
+ ls_niter_minor=ls_niter_minor)
+ IF (error) THEN
+ CALL errorMessage("oorb / random-walk_ranging", &
+ "TRACE BACK (5)", 1)
+ STOP
+ END IF
+ IF (get_cl_option("--epoch-mjd-tt=", .FALSE.)) THEN
+ CALL NULLIFY(epoch)
+ ! New epoch given as MJD TT
+ mjd_tt = get_cl_option("--epoch-mjd-tt=", 0.0_bp)
+ CALL NEW(epoch, mjd_tt, "TT")
+ IF (error) THEN
+ CALL errorMessage("oorb / random-walk_ranging", &
+ "TRACE BACK (10)", 1)
+ STOP
+ END IF
+ END IF
+
+ ! Print header before printing first orbit:
+ first = .TRUE.
+ DO i=1,SIZE(obss_sep,dim=1)
+ id = getID(obss_sep(i))
+ IF (error) THEN
+ CALL errorMessage("oorb / random-walk_ranging", &
+ "TRACE BACK (15)", 1)
+ STOP
+ END IF
+ IF (info_verb >= 2) THEN
+ WRITE(stdout,"(1X,I0,3(A),1X,I0,A,I0)") i, &
+ ". observation set (", TRIM(id), ")."
+ END IF
+ nobs = getNrOfObservations(obss_sep(i))
+ IF (error) THEN
+ CALL errorMessage("oorb / random-walk_ranging", &
+ "TRACE BACK (20)", 1)
+ STOP
+ END IF
+ dt = getObservationalTimespan(obss_sep(i))
+ IF (error) THEN
+ CALL errorMessage("oorb / random-walk_ranging", &
+ "TRACE BACK (40)", 1)
+ STOP
+ END IF
+ IF (info_verb >= 2) THEN
+ WRITE(stdout,"(4X,A,1X,F5.2,1X,A)") &
+ "Observational arc: ", dt, "days"
+ END IF
+ IF (sor_rho_init(3) > HUGE(sor_rho_init(3))/2) THEN
+ ! Initialize the rho2-rho1 range based on the observational
timespan.
+ IF (dt > 10) THEN
+ sor_rho_init(3) = -0.5_bp
+ ELSE
+ sor_rho_init(3) = -0.05_bp*dt
+ END IF
+ sor_rho_init(4) = -sor_rho_init(3)
+ END IF
+ !chi2_min assuming the error model is perfect
+ obs_masks => getObservationMasks(obss_sep(i))
+ chi2_min_init=REAL(COUNT(obs_masks), bp)
+ print*,"Chi2min_init ", chi2_min_init
+ IF (ASSOCIATED(id_arr_storb_in) .AND. ALLOCATED(storb_arr_in)) THEN
+ DO j=1,SIZE(id_arr_storb_in)
+ IF (TRIM(id) == TRIM(id_arr_storb_in(j))) THEN
+ EXIT
+ END IF
+ END DO
+ IF (j <= SIZE(id_arr_storb_in)) THEN
+ CALL setParameters(storb_arr_in(j), &
+ dyn_model=dyn_model_init, &
+ perturbers=perturbers, &
+ integrator=integrator_init, &
+ integration_step=integration_step_init, &
+ accept_multiplier=accwin_multiplier, &
+ apriori_rho_min=apriori_rho_min, &
+ set_acceptance_window=.FALSE.)
+ IF (error) THEN
+ CALL errorMessage("oorb / random-walk_ranging", &
+ "TRACE BACK (30)", 1)
+ STOP
+ END IF
+ CALL constrainRangeDistributions(storb_arr_in(j),
obss_sep(i))
+ IF (error) THEN
+ CALL errorMessage("oorb / random-walk_ranging", &
+ "TRACE BACK (31)", 1)
+ STOP
+ END IF
+ ! Allow user to decide whether rho bounds are updated from
orbit file
+ ! (but always update chi2min if orbit file is present)
+ IF (ALL(sor_iterate_bounds)) then
+ CALL setRangeBounds(storb_arr_in(j))
+ IF (error) THEN
+ CALL errorMessage("oorb / random-walk_ranging", &
+ "TRACE BACK (32)", 1)
+ STOP
+ END IF
+
+ sor_rho_init(1:4) = getRangeBounds(storb_arr_in(j))
+ IF (error) THEN
+ CALL errorMessage("oorb / random-walk_ranging", &
+ "TRACE BACK (33)", 1)
+ STOP
+ END IF
+ end if
+ ! Initialize chi2min
+ orb_arr_cmp => getSampleOrbits(storb_arr_in(j))
+ IF (error) THEN
+ CALL errorMessage("oorb / random-walk_ranging", &
+ "TRACE BACK (34)", 1)
+ STOP
+ END IF
+ chi2_min=HUGE(chi2_min)
+ DO k=1,SIZE(orb_arr_cmp,dim=1)
+ chi2=getChi2(storb_arr_in(j), orb_arr_cmp(k))
+ if (chi2 < chi2_min_init) chi2_min_init=chi2
+ END DO
+ END IF
+ END IF
+! print*,chi2_min_init
+ IF (.NOT.exist(epoch)) THEN
+ CALL NULLIFY(t)
+ obs = getObservation(obss_sep(i),1)
+ IF (error) THEN
+ CALL errorMessage("oorb / random-walk_ranging", &
+ "TRACE BACK (45)", 1)
+ STOP
+ END IF
+ t = getTime(obs)
+ IF (error) THEN
+ CALL errorMessage("oorb / random-walk_ranging", &
+ "TRACE BACK (50)", 1)
+ STOP
+ END IF
+ CALL NULLIFY(obs)
+ mjd = getMJD(t, "tt")
+ IF (error) THEN
+ CALL errorMessage("oorb / random-walk_ranging", &
+ "TRACE BACK (55)", 1)
+ STOP
+ END IF
+ CALL NULLIFY(t)
+ IF (get_cl_option("--exact-mid-epoch", .FALSE.)) THEN
+ mjd = mjd+dt/2.0_bp
+ ELSE
+ mjd = REAL(NINT(mjd+dt/2.0_bp),bp)
+ END IF
+ CALL NEW(t, mjd, "tt")
+ IF (error) THEN
+ CALL errorMessage("oorb / random-walk_ranging", &
+ "TRACE BACK (60)", 1)
+ STOP
+ END IF
+ ELSE
+ CALL NULLIFY(t)
+ t = copy(epoch)
+ END IF
+ CALL NEW(storb, obss_sep(i))
+ IF (error) THEN
+ CALL errorMessage("oorb / random-walk_ranging", &
+ "TRACE BACK (65)", 1)
+ STOP
+ END IF
+ CALL setParameters(storb, &
+ dyn_model=dyn_model, &
+ perturbers=perturbers, &
+ integrator=integrator, &
+ integration_step=integration_step, &
+ outlier_rejection=outlier_rejection_prm, &
+ outlier_multiplier=outlier_multiplier_prm, &
+ t_inv=t, &
+ element_type=element_type_comp_prm, &
+ dchi2_rejection = dchi2_rejection, &
+ dchi2_max = dchi2_max, &
+ regularized_pdf = regularized, &
+ accept_multiplier=accwin_multiplier, &
+ apriori_a_max=apriori_a_max, apriori_a_min=apriori_a_min, &
+ apriori_periapsis_max=apriori_periapsis_max, &
+ apriori_periapsis_min=apriori_periapsis_min, &
+ apriori_apoapsis_max=apriori_apoapsis_max, &
+ apriori_apoapsis_min=apriori_apoapsis_min, &
+ apriori_rho_min=apriori_rho_min, &
+ sor_2point_method=sor_2point_method, &
+ sor_2point_method_sw=sor_2point_method_sw, &
+ sor_norb_sw=sor_norb_sw, sor_ntrial_sw=sor_ntrial_sw, &
+ sor_norb=sor_norb, sor_ntrial=sor_ntrial, &
+ sor_rho1_l=sor_rho_init(1), sor_rho1_u=sor_rho_init(2), &
+ sor_rho2_l=sor_rho_init(3), sor_rho2_u=sor_rho_init(4), &
+ sor_niter=sor_niter, sor_iterate_bounds=sor_iterate_bounds, &
+ sor_random_obs_selection=.FALSE., &
+ gaussian_pdf=gaussian_rho, &
+ generat_multiplier=generat_multiplier, &
+ sor_generat_offset=sor_genwin_offset, &
+ chi2_min = chi2_min_init)
+ IF (error) THEN
+ CALL errorMessage("oorb / random-walk_ranging", &
+ "TRACE BACK (70)", 1)
+ STOP
+ END IF
+ SELECT CASE (sor_type_prm)
+ CASE (1)
+ CALL MCMCRanging4(storb) ! MCMCRanging3 = "basic" MCMC
+ CASE (2)
+!!$ CALL setParameters(storb, &
+!!$ dyn_model_init=dyn_model_init, &
+!!$ perturbers=perturbers, &
+!!$ integrator_init=integrator_init, &
+!!$ integration_step_init=integration_step_init, &
+!!$ accept_multiplier=accwin_multiplier, &
+!!$ apriori_rho_min=apriori_rho_min, &
+!!$ set_acceptance_window=.FALSE.)
+ !CALL autoMCMCRanging3(storb)
+ CALL errorMessage("oorb / random-walk_ranging", &
+ "Automatic ranging not implemented for
random-walk_ranging:", 1)
+ CASE default
+ CALL errorMessage("oorb / random-walk_ranging", &
+ "Unknown type of ranging:", 1)
+ IF (err_verb >= 1) THEN
+ WRITE(stderr,"(2X,I0)") sor_type_prm
+ END IF
+ STOP
+ END SELECT
+ IF (error) THEN
+ CALL errorMessage("oorb / random-walk_ranging", &
+ "MCMC failed:", 1)
+ IF (err_verb >= 1) THEN
+ WRITE(stderr,"(3A,1X,I0)") "ID: ", TRIM(id), &
+ " and number of observations: ", nobs
+ END IF
+ error = .FALSE.
+ CALL NEW(out_file, "problematic_observation_sets" // "." // &
+ TRIM(observation_format_out))
+ CALL setPositionAppend(out_file)
+ CALL OPEN(out_file)
+ IF (error) THEN
+ CALL errorMessage("oorb / random-walk_ranging", &
+ "TRACE BACK (90)", 1)
+ STOP
+ END IF
+ CALL writeObservationFile(obss_sep(i), getUnit(out_file), &
+ TRIM(observation_format_out))
+ IF (error) THEN
+ CALL errorMessage("oorb / random-walk_ranging", &
+ "TRACE BACK (95)", 1)
+ STOP
+ END IF
+ CALL NULLIFY(out_file)
+ IF (error) THEN
+ CALL errorMessage("oorb / random-walk_ranging", &
+ "TRACE BACK (100)", 1)
+ STOP
+ END IF
+ ELSE
+ IF (info_verb >= 2) THEN
+ WRITE(stdout,"(3(1X,A))") "Mcmc inversion for object", &
+ TRIM(id), "is ready."
+ END IF
+ CALL NEW(out_file, TRIM(id) // ".sor")
+ CALL OPEN(out_file)
+ IF (error) THEN
+ CALL errorMessage("oorb / random-walk_ranging", &
+ "TRACE BACK (105)", 1)
+ STOP
+ END IF
+ ! WRITE OBSERVATIONS:
+ CALL writeObservationFile(obss_sep(i), getUnit(out_file), &
+ TRIM(observation_format_out))
+ IF (error) THEN
+ CALL errorMessage("oorb / random-walk_ranging", &
+ "TRACE BACK (110)", 1)
+ STOP
+ END IF
+ WRITE(getUnit(out_file),"(A)") "#"
+ ! WRITE MCMC PARAMETERS
+ CALL writeSORResults(storb, obss_sep(i), getUnit(out_file))
+ IF (error) THEN
+ CALL errorMessage("oorb / mrandom-walk_ranging", &
+ "TRACE BACK (115)", 1)
+ STOP
+ END IF
+ CALL NULLIFY(out_file)
+ ! WRITE ORBITAL-ELEMENT PDF
+ orb_arr_cmp => getSampleOrbits(storb)
+ IF (error) THEN
+ CALL errorMessage("oorb / random-walk_ranging", &
+ "TRACE BACK (125)", 1)
+ STOP
+ END IF
+ pdf_arr_cmp => getPDFValues(storb, element_type_comp_prm)
+ IF (error) THEN
+ CALL errorMessage("oorb / random-walk_ranging", &
+ "TRACE BACK (130)", 1)
+ STOP
+ END IF
+ rchi2_arr_cmp => getReducedChi2Distribution(storb)
+ IF (error) THEN
+ CALL errorMessage("oorb / random-walk_ranging", &
+ "TRACE BACK (135)", 1)
+ STOP
+ END IF
+ CALL getResults(storb, repetition_arr_cmp=repetition_arr_cmp)
+!!$ CALL getResults(storb, &
+!!$ reg_apr_arr=reg_apr_arr_cmp, &
+!!$ jac_arr=jac_arr_cmp)
+!!$ IF (error) THEN
+!!$ CALL errorMessage("oorb / mcmcranging ", &
+!!$ "TRACE BACK (140)", 1)
+!!$ STOP
+!!$ END IF
+ IF (separately) THEN
+ CALL NEW(out_file, TRIM(id) // ".orb")
+ CALL OPEN(out_file)
+ IF (error) THEN
+ CALL errorMessage("oorb / random-walk_ranging", &
+ "TRACE BACK (145)", 1)
+ STOP
+ END IF
+ lu_orb_out = getUnit(out_file)
+ END IF
+ DO j=1,SIZE(orb_arr_cmp,dim=1)
+ IF (orbit_format_out == "orb") THEN
+ CALL writeOpenOrbOrbitFile(lu_orb_out, &
+ print_header=j==1, &
+ element_type_out=element_type_out_prm, &
+ id=id, &
+ orb=orb_arr_cmp(j), &
+ element_type_pdf=element_type_comp_prm, &
+ pdf=pdf_arr_cmp(j), &
+ rchi2=rchi2_arr_cmp(j), &
+ repetitions=repetition_arr_cmp(j), &
+ mjd=mjd_epoch)
+!!$ reg_apr=reg_apr_arr_cmp(j), &
+!!$ jac_sph_inv=jac_arr_cmp(j,1), &
+!!$ jac_car_kep=jac_arr_cmp(j,2), &
+!!$ jac_equ_kep=jac_arr_cmp(j,3))
+ ELSE IF (orbit_format_out == "des") THEN
+ CALL errorMessage("oorb / random-walk_rangingg", &
+ "DES format not yet supported for MCMC output.", 1)
+ STOP
+ END IF
+ IF (error) THEN
+ CALL errorMessage("oorb / random-walk_ranging", &
+ "TRACE BACK (150)", 1)
+ STOP
+ END IF
+ END DO
+ IF (separately) THEN
+ CALL NULLIFY(out_file)
+ IF (compress) THEN
+ CALL system("gzip -f " // TRIM(id) // ".orb")
+ END IF
+ END IF
+ ! WRITE RESIDUALS
+ IF (write_residuals) THEN
+ CALL NEW(out_file, TRIM(out_fname) // ".res")
+ CALL setPositionAppend(out_file)
+ CALL OPEN(out_file)
+ IF (error) THEN
+ CALL errorMessage("oorb / random-walk_ranging", &
+ "TRACE BACK (155)", 1)
+ STOP
+ END IF
+ CALL writeResiduals(storb, obss_sep(i), getUnit(out_file))
+ IF (error) THEN
+ CALL errorMessage("oorb / random-walk_ranging", &
+ "TRACE BACK (160)", 1)
+ STOP
+ END IF
+ CALL NULLIFY(out_file)
+ END IF
+ IF (plot_results) THEN
+ CALL toString(dt, str, error, frmt="(F10.2)")
+ IF (error) THEN
+ CALL errorMessage("oorb / random-walk_ranging", &
+ "TRACE BACK (165)", 1)
+ STOP
+ END IF
+ str = TRIM(id) // "_"// TRIM(str)
+ CALL makeResidualStamps(storb, obss_sep(i), TRIM(str) // &
+ "_random-walk_ranging_residual_stamps.eps")
+ IF (error) THEN
+ CALL errorMessage("oorb / random-walk_ranging", &
+ "TRACE BACK (170)", 1)
+ STOP
+ END IF
+ IF (compress) THEN
+ CALL system("gzip -f " // TRIM(str)
// "_random-walk_ranging_residual_stamps.eps")
+ END IF
+ IF (plot_open) THEN
+ CALL system("gv " // TRIM(str)
// "_random-walk_ranging_residual_stamps.eps* &")
+ END IF
+ ALLOCATE(elements_arr(SIZE(orb_arr_cmp,dim=1),7), stat=err)
+ IF (err /= 0) THEN
+ CALL errorMessage("oorb / random-walk_ranging", &
+ "Could not allocate memory (3).", 1)
+ STOP
+ END IF
+ CALL NEW(tmp_file,
TRIM(str)// "_random-walk_ranging_orbits.out")
+ CALL OPEN(tmp_file)
+ IF (error) THEN
+ CALL errorMessage("oorb / random-walk_ranging ", &
+ "TRACE BACK (175)", 1)
+ STOP
+ END IF
+ DO j=1,SIZE(orb_arr_cmp,dim=1)
+ IF (element_type_comp_prm == "cartesian") THEN
+ CALL rotateToEcliptic(orb_arr_cmp(j))
+ END IF
+ elements_arr(j,1:6) = getElements(orb_arr_cmp(j),
element_type_comp_prm)
+ IF (error) THEN
+ CALL errorMessage("oorb / random-walk_ranging", &
+ "TRACE BACK (180)", 1)
+ STOP
+ END IF
+ IF (element_type_comp_prm == "keplerian") THEN
+ elements_arr(j,3:6) = elements_arr(j,3:6)/rad_deg
+ END IF
+ elements_arr(j,7) = pdf_arr_cmp(j)
+ t = getTime(orb_arr_cmp(j))
+ IF (error) THEN
+ CALL errorMessage("oorb / random-walk_ranging ", &
+ "TRACE BACK (185)", 1)
+ STOP
+ END IF
+ WRITE(getUnit(tmp_file),*) &
+ elements_arr(j,1:6), &
+ pdf_arr_cmp(j), &
+ getCalendarDateString(t,"tdt")
+ IF (error) THEN
+ CALL errorMessage("oorb / random-walk_ranging ", &
+ "TRACE BACK (190)", 1)
+ STOP
+ END IF
+ CALL NULLIFY(t)
+ END DO
+ CALL NULLIFY(tmp_file)
+ CALL NEW(tmp_file, TRIM(str) // &
+ "_random-walk_ranging_sample_standard_deviations.out")
+ CALL OPEN(tmp_file)
+ IF (error) THEN
+ CALL errorMessage("oorb / random-walk_ranging ", &
+ "TRACE BACK (195)", 1)
+ STOP
+ END IF
+ WRITE(getUnit(tmp_file), "(F22.15,1X)", &
+ advance="no") &
+ getObservationalTimespan(obss_sep(i))
+ IF (error) THEN
+ CALL errorMessage("oorb / random-walk_ranging ", &
+ "TRACE BACK (200)", 1)
+ STOP
+ END IF
+ DO j=1,6
+ CALL moments(elements_arr(:,j), &
+ pdf=pdf_arr_cmp, std_dev=stdev, errstr=errstr)
+ WRITE(getUnit(tmp_file), "(F22.15,1X)", &
+ advance="no") stdev
+ stdev_arr(j)=stdev
+ END DO
+ WRITE(getUnit(tmp_file),*)
+ CALL NULLIFY(tmp_file)
+ DEALLOCATE(elements_arr, stat=err)
+ IF (err /= 0) THEN
+ CALL errorMessage("oorb / random-walk_ranging", &
+ "Could not deallocate memory (5).", 1)
+ STOP
+ END IF
+ ! Make plot using gnuplot:
+ CALL system("cp " // TRIM(str) // &
+ "_random-walk_ranging_orbits.out sor_orbits.out")
+ IF (element_type_comp_prm == "cartesian") THEN
+ CALL system("gnuplot " // TRIM(gnuplot_scripts_dir)
// "/sor_plot_car.gp")
+ ELSE
+ if (stdev_arr(1) < 1.0_bp) then
+ CALL system("gnuplot " // TRIM(gnuplot_scripts_dir)
// "/sor_plot_kep_nolog.gp")
+ else
+ CALL system("gnuplot " // TRIM(gnuplot_scripts_dir)
// "/sor_plot_kep.gp")
+ end if
+ END IF
+ CALL system("cp sor_results.eps " // TRIM(str) // &
+ "_random-walk_ranging_" // TRIM(element_type_comp_prm)
// &
+ "_results.eps")
+ CALL system("rm -f sor_orbits.out sor_results.eps " // &
+ TRIM(str) // "_random-walk_ranging_orbits.out ")
+!!$ // TRIM(str) // &
+!!$ "_random-walk_ranging_sample_standard_deviations.out")
+ IF (compress) THEN
+ CALL system("gzip -f " // TRIM(str)
// "_random-walk_ranging_" // &
+ TRIM(element_type_comp_prm) // "_results.eps")
+ END IF
+ IF (plot_open) THEN
+ CALL system("gv " // TRIM(str) // "_random-walk_ranging_"
// &
+ TRIM(element_type_comp_prm) // &
+ "_results.eps* &")
+ END IF
+ END IF
+ DO j=1,SIZE(orb_arr_cmp)
+ CALL NULLIFY(orb_arr_cmp(j))
+ END DO
+ DEALLOCATE(orb_arr_cmp, stat=err)
+ DEALLOCATE(pdf_arr_cmp, stat=err)
+ DEALLOCATE(rchi2_arr_cmp, stat=err)
+ DEALLOCATE(repetition_arr_cmp, stat=err)
+!!$ DEALLOCATE(reg_apr_arr_cmp, stat=err)
+!!$ DEALLOCATE(jac_arr_cmp, stat=err)
+ IF (err /= 0) THEN
+ CALL errorMessage("oorb / random-walk_ranging", &
+ "Could not deallocate memory (10).", 1)
+ STOP
+ END IF
+ IF (info_verb >= 2) THEN
+ WRITE(stdout,"(3(1X,A))") "Object", &
+ TRIM(id), "successfully processed."
+ END IF
+
+ END IF
+ CALL NULLIFY(storb)
+ CALL NULLIFY(orb)
+ IF (info_verb > 2) THEN
+ WRITE(stdout,*)
+ WRITE(stdout,*)
+ END IF
+ END DO
+
+
+!*********************
CASE ("simplex")
@@ -6499,25 +7107,25 @@
IF (info_verb >= 2 .AND. dyn_model /= "2-body") THEN
CALL getCalendarDate(epoch0, "TT", year0, month0, day0)
CALL getCalendarDate(epoch, "TT", year1, month1, day1)
- WRITE(stderr,'(A)') ""
- WRITE(stderr,'(A,I0,A,I0,A,F0.5,A,I0,A,I0,A,F0.5,A)') &
- "Planetary encounters by object " //
TRIM(id_arr_storb_in(i)) // &
- " between ", year0, "-", month0, "-", day0, " and ",
&
- year1, "-", month1, "-", day1, ":"
+!~ WRITE(stderr,'(A)') ""
+!~ WRITE(stderr,'(A,I0,A,I0,A,F0.5,A,I0,A,I0,A,F0.5,A)') &
+!~ "Planetary encounters by object " //
TRIM(id_arr_storb_in(i)) // &
+!~ " between ", year0, "-", month0, "-", day0, "
and ", &
+!~ year1, "-", month1, "-", day1, ":"
DO j=1,SIZE(encounters,dim=1)
- WRITE(stderr,'(A,I0,A)') "Orbit #", j, ":"
+!~ WRITE(stderr,'(A,I0,A)') "Orbit #", j, ":"
DO k=1,SIZE(encounters,dim=2)
IF (encounters(j,k,2) < 1.1_bp) THEN
WRITE(stderr,'(A22,1X,A7,1X,A6,1X,F15.6,1X,A,1X,F0.6,1X,A)') &
"!! I-M-P-A-C-T !! with",
planetary_locations(k), &
"at MJD", encounters(j,k,1), "TT given a
stepsize of", &
encounters(j,k,4), "days."
- ELSE
-
WRITE(stderr,'(A22,1X,A7,1X,A6,1X,F15.6,1X,A,1X,F11.8,1X,A,1X,F0.6,1X,A)') &
- "Closest encounter with",
planetary_locations(k), &
- "at MJD", encounters(j,k,1), &
- "TT at a distance of", encounters(j,k,3), &
- "AU given a stepsize of",
encounters(j,k,4), "days."
+!~ ELSE
+!~
WRITE(stderr,'(A22,1X,A7,1X,A6,1X,F15.6,1X,A,1X,F11.8,1X,A,1X,F0.6,1X,A)') &
+!~ "Closest encounter with",
planetary_locations(k), &
+!~ "at MJD", encounters(j,k,1), &
+!~ "TT at a distance of",
encounters(j,k,3), &
+!~ "AU given a stepsize of",
encounters(j,k,4), "days."
END IF
END DO
END DO
@@ -6799,13 +7407,13 @@
STOP
END IF
IF (info_verb >= 2 .AND. dyn_model /= "2-body") THEN
- CALL getCalendarDate(epoch0, "TT", year0, month0, day0)
- CALL getCalendarDate(epoch, "TT", year1, month1, day1)
- WRITE(stderr,'(A)') ""
-
WRITE(stderr,'(A,I0,A,I0,A,I0,A,F0.5,A,I0,A,I0,A,F0.5,A)') &
- "Planetary encounters by object " //
TRIM(id_arr_in(i)) // &
- " (orbit #", i, ") between ", year0, "-", month0, &
- "-", day0, " and ", year1, "-", month1, "-",
day1, ":"
+!~ CALL getCalendarDate(epoch0, "TT", year0, month0, day0)
+!~ CALL getCalendarDate(epoch, "TT", year1, month1, day1)
+!~ WRITE(stderr,'(A)') ""
+!~
WRITE(stderr,'(A,I0,A,I0,A,I0,A,F0.5,A,I0,A,I0,A,F0.5,A)') &
+!~ "Planetary encounters by object " //
TRIM(id_arr_in(i)) // &
+!~ " (orbit #", i, ") between ", year0, "-", month0,
&
+!~ "-", day0, " and ", year1, "-", month1, "-",
day1, ":"
DO j=1,SIZE(encounters,dim=2)
IF (encounters(1,j,2) < 1.1_bp) THEN
! encounters(1,j,2) == 1 == impact
@@ -6813,12 +7421,11 @@
"!! I-M-P-A-C-T !! with",
planetary_locations(j), &
"at MJD", encounters(1,j,1), "TT given a
stepsize of", &
encounters(1,j,4), "days."
- ELSE
- ! encounters(1,j,2) == 2 == non-impact
-
WRITE(stderr,'(A22,1X,A7,1X,A6,1X,F15.6,1X,A,1X,F11.8,1X,A,1X,F0.6,1X,A)') &
- "Closest encounter with",
planetary_locations(j), &
- "at MJD", encounters(1,j,1), "TT at a distance
of", &
- encounters(1,j,3), "AU given a stepsize of",
encounters(1,j,4), "days."
+!~ ELSE
+!~
WRITE(stderr,'(A22,1X,A7,1X,A6,1X,F15.6,1X,A,1X,F11.8,1X,A,1X,F0.6,1X,A)') &
+!~ "Closest encounter with",
planetary_locations(j), &
+!~ "at MJD", encounters(1,j,1), "TT at a
distance of", &
+!~ encounters(1,j,3), "AU given a stepsize
of", encounters(1,j,4), "days."
END IF
END DO
END IF
Reply all
Reply to author
Forward
0 new messages