[oorb] r283 committed - Last small updates before export to GitHub.

[oo...@googlecode.com]

Revision: 283
Author: mikael.granvik
Date: Mon Jun 15 22:42:10 2015 UTC
Log: Last small updates before export to GitHub.
https://code.google.com/p/oorb/source/detail?r=283

Modified:
/trunk/classes/File_class.f90
/trunk/classes/Observation_class.f90
/trunk/classes/Observations_class.f90
/trunk/classes/PhysicalParameters_class.f90
/trunk/classes/StochasticOrbit_class.f90
/trunk/data/ET-UT.dat
/trunk/data/JPL_ephemeris/ephtester.f90
/trunk/data/updateOBSCODE
/trunk/gnuplot/vov_plot_car.gp
/trunk/main/io.f90
/trunk/main/oorb.conf
/trunk/main/oorb.f90
/trunk/modules/planetary_data.f90
/trunk/modules/statistics.f90

=======================================
--- /trunk/classes/File_class.f90 Wed Feb 15 11:55:33 2012 UTC
+++ /trunk/classes/File_class.f90 Mon Jun 15 22:42:10 2015 UTC
@@ -1,6 +1,6 @@
!====================================================================!
! !
-! Copyright 2002,2003,2004,2005,2006,2007,2008,2009,2010,2011,2012 !
+! Copyright 2002-2014,2015 !
! Mikael Granvik, Jenni Virtanen, Karri Muinonen, Teemu Laakso, !
! Dagmara Oszkiewicz !
! !
@@ -66,7 +66,7 @@
!! </pre>
!!
!! @author MG
-!! @version 2012-01-17
+!! @version 2015-06-15
!!
!! @see Observations_class
!! @see Time_class
@@ -437,7 +437,7 @@

i = 0
DO WHILE (LEN_TRIM(line) /= 0)
- Line_ = " "
+ line_ = " "
IF (line(1:1) /= " ") THEN
i = i + 1
indx = INDEX(line," ")
=======================================
--- /trunk/classes/Observation_class.f90 Fri Aug 22 06:45:25 2014 UTC
+++ /trunk/classes/Observation_class.f90 Mon Jun 15 22:42:10 2015 UTC
@@ -1,6 +1,6 @@
!====================================================================!
! !
-! Copyright 2002-2013,2014 !
+! Copyright 2002-2014,2015 !
! Mikael Granvik, Jenni Virtanen, Karri Muinonen, Teemu Laakso, !
! Dagmara Oszkiewicz !
! !
@@ -29,7 +29,7 @@
!! @see Observations_class
!!
!! @author MG, JV
-!! @version 2014-08-22
+!! @version 2015-06-15
!!
MODULE Observation_cl

@@ -537,12 +537,15 @@
!!
!! Returns error.
!!
- SUBROUTINE addMultinormalDeviate_Obs(this, mean, covariance)
+ SUBROUTINE addMultinormalDeviate_Obs(this, mean, covariance,
combined_covariance)

IMPLICIT NONE
TYPE (Observation), INTENT(inout) :: this
REAL(bp), DIMENSION(6), INTENT(in) :: mean
REAL(bp), DIMENSION(6,6), INTENT(in) :: covariance
+ LOGICAL, OPTIONAL, INTENT(in) :: combined_covariance
+
+ LOGICAL :: combined_covariance_

IF (.NOT. this%is_initialized) THEN
error = .TRUE.
@@ -551,14 +554,28 @@
RETURN
END IF

+ IF (PRESENT(combined_covariance)) THEN
+ combined_covariance_ = combined_covariance
+ ELSE
+ combined_covariance_ = .TRUE.
+ END IF
+
CALL addMultinormalDeviate(this%obs_scoord, mean, covariance)
IF (error) THEN
CALL errorMessage("Observation / addMultinormalDeviate", &
"TRACE BACK", 1)
RETURN
END IF
- ! New noise estimate = old noise estimate + added noise
- this%covariance = this%covariance + covariance
+
+ IF (combined_covariance_) THEN
+ ! New noise covariance = old noise covariance + added noise
covariance
+ this%covariance = this%covariance + covariance
+ ELSE
+ ! New noise covariance = added noise covariance
+ !
+ ! (this is aimed for cases when noise is added to noiseless
synthetic astrometry)
+ this%covariance = covariance
+ END IF

END SUBROUTINE addMultinormalDeviate_Obs

=======================================
--- /trunk/classes/Observations_class.f90 Wed Apr 17 05:07:29 2013 UTC
+++ /trunk/classes/Observations_class.f90 Mon Jun 15 22:42:10 2015 UTC
@@ -1,6 +1,6 @@
!====================================================================!
! !
-! Copyright 2002-2012,2013 !
+! Copyright 2002-2014,2015 !
! Mikael Granvik, Jenni Virtanen, Karri Muinonen, Teemu Laakso, !
! Dagmara Oszkiewicz !
! !
@@ -54,7 +54,7 @@
!! @see StochasticOrbit_class
!!
!! @author MG, JV
-!! @version 2013-04-16
+!! @version 2015-06-15
!!
MODULE Observations_cl

@@ -3104,18 +3104,20 @@
CHARACTER(len=132) :: record
CHARACTER(len=124) :: line1, line2, line_, str, filter, obstype
CHARACTER(len=7) :: nr_str
- CHARACTER(len=4) :: timescale, obsy_code
+ CHARACTER(len=4) :: timescale, obsy_code, strk_class
REAL(bp), DIMENSION(:,:), POINTER :: planeph
REAL(bp), DIMENSION(:,:), ALLOCATABLE :: element_arr
REAL(bp), DIMENSION(6,6) :: covariance
REAL(bp), DIMENSION(6) :: coordinates, stdev_, mean
- REAL(bp), DIMENSION(3) :: position, velocity
+ REAL(bp), DIMENSION(3) :: position, velocity, pos1, pos2
REAL(bp) :: day, sec, arcsec, mag, ra, dec, jd, mjd_utc, dt, &
ecl_lon, ecl_lat, angscan, pos_unc_along, pos_unc_across, &
vel_unc_along, vel_unc_across, rot_angle, correlation, &
- rmin, rarcmin, mag_unc, s2n, mjd_tt
+ rmin, rarcmin, mag_unc, s2n, mjd_tt, mjd_tcb, strk_len, &
+ strk_len_unc, strk_direction, strk_direction_unc
INTEGER :: i, j, err, year, month, hour, min, deg, arcmin, &
- nlines, coord_unit, indx, iobs, irecord, norb, ccd
+ nlines, coord_unit, indx, iobs, irecord, norb, ccd, &
+ border1, border2
LOGICAL, DIMENSION(6) :: obs_mask
LOGICAL :: discovery, converttonewformat

@@ -4257,8 +4259,8 @@
obs_mask = (/ .FALSE., .TRUE., .TRUE., .FALSE., .FALSE., .FALSE. /)

! Origin of observation dates: 1.0 Jan 2010 ( = JD 2455197.5 = MJD
55197.0)
- mjd_utc = 55197.0_bp
- i = 0
+ mjd_tcb = 55197.0_bp
+ i = 0
DO

READ(getUnit(obsf), "(A)", iostat=err) line
@@ -4283,9 +4285,6 @@
END IF
covariance(3,2) = covariance(2,3)

- covariance(2,:) = covariance(2,:) * COS(position(3))
- covariance(:,2) = covariance(:,2) * COS(position(3))
-
i = i + 1
IF (i == 1) THEN
discovery = .TRUE.
@@ -4293,7 +4292,7 @@
discovery = .FALSE.
END IF

- CALL NEW(t, mjd_utc + dt, "utc")
+ CALL NEW(t, mjd_tcb + dt, "TCB")
IF (error) THEN
CALL errorMessage("Observations / readObservationFile", &
"TRACE BACK (115)", 1)
@@ -4360,6 +4359,243 @@

END DO

+ CASE ("geo")
+
+ covariance = 0.0_bp
+ position = 0.0_bp
+ velocity = 0.0_bp
+ obs_mask = (/ .FALSE., .TRUE., .TRUE., .FALSE., .FALSE., .FALSE. /)
+
+ i = 0
+ DO
+
+ READ(getUnit(obsf), "(A)", iostat=err) line
+ IF (err > 0) THEN
+ error = .TRUE.
+ CALL errorMessage("Observations / readObservationFile", &
+ "Error while reading observations from file (1).", 1)
+ RETURN
+ ELSE IF (err < 0) THEN ! end-of-file
+ EXIT
+ ELSE IF (line(1:1) == "#") THEN
+ CYCLE
+ END IF
+ READ(line, *, iostat=err) number, obsy_code, mjd_utc,
position(2), &
+ position(3), covariance(2,2), covariance(3,3),
coordinates(1:6)
+ IF (err /= 0) THEN
+ error = .TRUE.
+ CALL errorMessage("Observations / readObservationFile", &
+ "Error while reading observations from file (2).", 1)
+ RETURN
+ END IF
+
+ i = i + 1
+ IF (i == 1) THEN
+ discovery = .TRUE.
+ ELSE
+ discovery = .FALSE.
+ END IF
+
+ CALL NEW(t, mjd_utc, "UTC")
+ IF (error) THEN
+ CALL errorMessage("Observations / readObservationFile", &
+ "TRACE BACK (115)", 1)
+ RETURN
+ END IF
+ obsy = getObservatory(obsies, "247")
+ IF (error) THEN
+ CALL errorMessage("Observations / readObservationFile", &
+ "TRACE BACK (120)", 1)
+ RETURN
+ END IF
+ ! Transform geocentric spacecraft coordinates to
+ ! heliocentric spacecraft coordinates
+ mjd_tt = getMJD(t, "TT")
+ IF (error) THEN
+ CALL errorMessage("Observations / readObservationFile", &
+ "TRACE BACK (121)", 1)
+ RETURN
+ END IF
+ planeph => JPL_ephemeris(mjd_tt, 3, 11, error)
+ IF (error) THEN
+ CALL errorMessage("Observations / readObservationFile", &
+ "TRACE BACK (122)", 1)
+ RETURN
+ END IF
+ CALL NEW(obsy_ccoord, coordinates + planeph(1,:), "equatorial",
t)
+ IF (error) THEN
+ CALL errorMessage("Observations / readObservationFile", &
+ "TRACE BACK (125)", 1)
+ RETURN
+ END IF
+ DEALLOCATE(planeph, stat=err)
+ CALL NEW(obs_scoord, position, velocity, "equatorial", t)
+ IF (error) THEN
+ CALL errorMessage("Observations / readObservationFile", &
+ "TRACE BACK (130)", 1)
+ RETURN
+ END IF
+ satellite_ccoord = getObservatoryCCoord(obsies, "500", t)
+ CALL rotateToEquatorial(satellite_ccoord)
+ CALL rotateToEquatorial(obsy_ccoord)
+ coordinates = getCoordinates(obsy_ccoord) -
getCoordinates(satellite_ccoord)
+ CALL NULLIFY(satellite_ccoord)
+ CALL NEW(satellite_ccoord, coordinates, "equatorial", t)
+ CALL NULLIFY(this%obs_arr(i))
+ CALL NEW(this%obs_arr(i), number=number, designation=" ", &
+ discovery=discovery, note1=" ", note2="V", &
+ obs_scoord=obs_scoord, covariance=covariance, &
+ obs_mask=obs_mask, mag=99.9_bp, filter=" ", &
+ obsy=obsy, obsy_ccoord=obsy_ccoord, &
+ satellite_ccoord=satellite_ccoord, &
+ coord_unit=2)
+ IF (error) THEN
+ CALL errorMessage("Observations / readObservationFile", &
+ "TRACE BACK (135)", 1)
+ RETURN
+ END IF
+
+ CALL NULLIFY(obs_scoord)
+ CALL NULLIFY(t)
+ CALL NULLIFY(obsy)
+ CALL NULLIFY(obsy_ccoord)
+ CALL NULLIFY(satellite_ccoord)
+
+ END DO
+
+ CASE ("strk")
+
+
+ DEALLOCATE(this%obs_arr, this%obs_note_arr, stat=err)
+ ALLOCATE(this%obs_arr(2*nlines), this%obs_note_arr(2*nlines),
stat=err)
+ IF (err /= 0) THEN
+ error = .TRUE.
+ CALL errorMessage("Observations / readObservationFile", &
+ "Could not allocate memory.", 1)
+ RETURN
+ END IF
+ this%obs_note_arr = " "
+
+ number = ""
+ covariance = 0.0_bp
+ position = 0.0_bp
+ velocity = 0.0_bp
+ obs_mask = (/ .FALSE., .TRUE., .TRUE., .FALSE., .FALSE., .FALSE. /)
+
+ i = 0
+ DO
+
+ READ(getUnit(obsf), "(A)", iostat=err) line
+ IF (err > 0) THEN
+ error = .TRUE.
+ CALL errorMessage("Observations / readObservationFile", &
+ "Error while reading observations from file (1).", 1)
+ RETURN
+ ELSE IF (err < 0) THEN ! end-of-file
+ EXIT
+ ELSE IF (line(1:1) == "#") THEN
+ CYCLE
+ END IF
+ READ(line, *, iostat=err) designation, mjd_utc, dt, &
+ obsy_code, position(2), covariance(2,2), position(3), &
+ covariance(3,3), covariance(2,3), strk_len, &
+ strk_len_unc, strk_direction, strk_direction_unc, &
+ mag, mag_unc, pos1(2), pos1(3), border1, pos2(2), &
+ pos2(3), border2, strk_class
+ IF (err /= 0) THEN
+ error = .TRUE.
+ CALL errorMessage("Observations / readObservationFile", &
+ "Error while reading observations from file (2).", 1)
+ RETURN
+ END IF
+
+ ! Change correlation to covariance
+ covariance(2,3) = covariance(2,3) * covariance(2,2) *
covariance(3,3)
+ covariance(3,2) = covariance(2,3)
+ ! Change standard deviations to variances
+ covariance(2,2) = covariance(2,2)**2
+ covariance(3,3) = covariance(3,3)**2
+ ! Change degrees to radians
+ position = position*rad_deg
+ covariance = covariance*(rad_deg**2)
+ pos1 = pos1*rad_deg
+ pos2 = pos2*rad_deg
+ strk_len = strk_len*rad_deg
+ strk_len_unc = strk_len_unc*rad_deg
+ strk_direction = strk_direction*rad_deg
+ strk_direction_unc = strk_direction_unc*rad_deg
+
+ ! Streak start point
+ i = i + 1
+ discovery = .FALSE.
+ CALL NEW(t, mjd_utc-dt/2.0_bp, "UTC")
+ IF (error) THEN
+ CALL errorMessage("Observations / readObservationFile", &
+ "TRACE BACK (115)", 1)
+ RETURN
+ END IF
+ CALL NEW(obs_scoord, pos1, velocity, "equatorial", t)
+ IF (error) THEN
+ CALL errorMessage("Observations / readObservationFile", &
+ "TRACE BACK (130)", 1)
+ RETURN
+ END IF
+ obsy = getObservatory(obsies, obsy_code)
+ obsy_ccoord = getObservatoryCCoord(obsies, obsy_code, t)
+ CALL NULLIFY(this%obs_arr(i))
+ CALL NEW(this%obs_arr(i), number=number,
designation=designation, &
+ discovery=discovery, note1=" ", note2="C", &
+ obs_scoord=obs_scoord, covariance=covariance, &
+ obs_mask=obs_mask, mag=mag, filter=" ", &
+ obsy=obsy, obsy_ccoord=obsy_ccoord)
+ IF (error) THEN
+ CALL errorMessage("Observations / readObservationFile", &
+ "TRACE BACK (135)", 1)
+ RETURN
+ END IF
+ CALL NULLIFY(obs_scoord)
+ CALL NULLIFY(t)
+ CALL NULLIFY(obsy)
+ CALL NULLIFY(obsy_ccoord)
+
+ ! Streak end point
+ i = i + 1
+ discovery = .FALSE.
+ CALL NEW(t, mjd_utc+dt/2.0_bp, "UTC")
+ IF (error) THEN
+ CALL errorMessage("Observations / readObservationFile", &
+ "TRACE BACK (115)", 1)
+ RETURN
+ END IF
+ CALL NEW(obs_scoord, pos2, velocity, "equatorial", t)
+ IF (error) THEN
+ CALL errorMessage("Observations / readObservationFile", &
+ "TRACE BACK (130)", 1)
+ RETURN
+ END IF
+ obsy = getObservatory(obsies, obsy_code)
+ obsy_ccoord = getObservatoryCCoord(obsies, obsy_code, t)
+ CALL NULLIFY(this%obs_arr(i))
+ CALL NEW(this%obs_arr(i), number=number,
designation=designation, &
+ discovery=discovery, note1=" ", note2="C", &
+ obs_scoord=obs_scoord, covariance=covariance, &
+ obs_mask=obs_mask, mag=mag, filter=" ", &
+ obsy=obsy, obsy_ccoord=obsy_ccoord)
+ IF (error) THEN
+ CALL errorMessage("Observations / readObservationFile", &
+ "TRACE BACK (135)", 1)
+ RETURN
+ END IF
+ CALL NULLIFY(obs_scoord)
+ CALL NULLIFY(t)
+ CALL NULLIFY(obsy)
+ CALL NULLIFY(obsy_ccoord)
+
+ END DO
+
+
+
+
CASE ("sentinel")

covariance = 0.0_bp
@@ -4394,9 +4630,6 @@
END IF
covariance(3,2) = covariance(2,3)

- covariance(2,:) = covariance(2,:) * COS(position(3))
- covariance(:,2) = covariance(:,2) * COS(position(3))
-
i = i + 1
IF (i == 1) THEN
discovery = .TRUE.
=======================================
--- /trunk/classes/PhysicalParameters_class.f90 Fri Aug 22 08:15:05 2014 UTC
+++ /trunk/classes/PhysicalParameters_class.f90 Mon Jun 15 22:42:10 2015 UTC
@@ -1,6 +1,6 @@
!====================================================================!
! !
-! Copyright 2002-2013,2014 !
+! Copyright 2002-2014,2015 !
! Mikael Granvik, Jenni Virtanen, Karri Muinonen, Teemu Laakso, !
! Dagmara Oszkiewicz !
! !
@@ -30,7 +30,7 @@
!! @see StochasticOrbit_class
!!
!! @author MG
-!! @version 2014-08-22
+!! @version 2015-06-15
!!
MODULE PhysicalParameters_cl

@@ -700,6 +700,12 @@
END IF
END DO
obsy_ccoord_arr => reallocate(obsy_ccoord_arr)
+ IF (.NOT.ASSOCIATED(obsy_ccoord_arr)) THEN
+ CALL warningMessage("PhysicalParameters / estimateHAndG", &
+ "0 observations contain brightness information -> " // &
+ "estimation of H and G is impossible.", 1)
+ RETURN
+ END IF
nmag = SIZE(obsy_ccoord_arr)
ALLOCATE(obs_mag_arr_(nmag), filter_arr_(nmag), obs_mag_unc_arr_(nmag))
j = 0
=======================================
--- /trunk/classes/StochasticOrbit_class.f90 Fri Apr 17 07:27:11 2015 UTC
+++ /trunk/classes/StochasticOrbit_class.f90 Mon Jun 15 22:42:10 2015 UTC
File is too large to display a diff.
=======================================
--- /trunk/data/ET-UT.dat Sun Jul 26 07:53:31 2009 UTC
+++ /trunk/data/ET-UT.dat Mon Jun 15 22:42:10 2015 UTC
@@ -197,10 +197,15 @@
1 1 2001 64.09
1 1 2002 64.30
1 1 2003 64.47
-1 1 2004 65.
-1 1 2005 66.
-1 1 2006 66.
-1 1 2007 67.
-1 1 2008 68.
+1 1 2004 64.57
+1 1 2005 64.6876
+1 1 2006 64.8452
+1 1 2007 65.1464
+1 1 2008 65.4574
+1 1 2009 65.7768
+1 1 2010 66.0699
+1 1 2011 66.3246
+1 1 2012 66.6030
+1 1 2013 66.9069
1 1 2050 114.
1 1 3000 250.
=======================================
--- /trunk/data/JPL_ephemeris/ephtester.f90 Fri Jan 22 09:14:43 2010 UTC
+++ /trunk/data/JPL_ephemeris/ephtester.f90 Mon Jun 15 22:42:10 2015 UTC
@@ -1,6 +1,6 @@
!====================================================================!
! !
-! Copyright 2002,2003,2004,2005,2006,2007,2008,2009,2010 !
+! Copyright 2002-2014,2015 !
! Mikael Granvik, Jenni Virtanen, Karri Muinonen, Teemu Laakso, !
! Dagmara Oszkiewicz !
! !
@@ -29,7 +29,7 @@
!! Pluto are included.
!!
!! @author MG
-!! @version 2010-01-21
+!! @version 2015-06-15
!!
PROGRAM ephtester

@@ -69,7 +69,8 @@
END IF
END IF
i = i + 1
- IF (ncenter == 0 .OR. ncenter > 11 .OR. ntarget > 11) THEN
+ ! IF (ncenter == 0 .OR. ncenter > 11 .OR. ntarget > 11) THEN
+ IF (ncenter == 0 .OR. ncenter > 13 .OR. ntarget > 13) THEN
CYCLE
END IF
mjd_tt = jd_tt - 2400000.5_rprec8
@@ -79,11 +80,19 @@
WRITE(0,*) "***** PREVIOUS NOT AN ERROR IN SW BUT IN JPL TEST FILE
*****"
CYCLE
END IF
- IF (ABS(ephemeris(1,ncoord)-correct_value) > 10.0E-13_rprec8) THEN
- WRITE(0,*) "***** COMPARISON ERROR OCCURRED *****"
- STOP
+
+ WRITE(*,"(I0,1X,A,1X,F10.2,3(1X,I2),3(1X,F20.15))") &
+ ieph_type, TRIM(str), jd_tt, ntarget, ncenter, ncoord, &
+ correct_value, ephemeris(1,ncoord), &
+ ephemeris(1,ncoord)-correct_value
+
+ ! IF (ABS(ephemeris(1,ncoord)-correct_value) > 10.0E-13_rprec8)
THEN
+ ! WRITE(0,*) "***** COMPARISON ERROR OCCURRED *****"
+ ! STOP
+ ! END IF
+ IF (ASSOCIATED(ephemeris)) THEN
+ DEALLOCATE(ephemeris,stat=err)
END IF
- DEALLOCATE(ephemeris)
END DO
CLOSE(12)
CALL JPL_ephemeris_nullify()
=======================================
--- /trunk/data/updateOBSCODE Wed Sep 9 02:20:21 2009 UTC
+++ /trunk/data/updateOBSCODE Mon Jun 15 22:42:10 2015 UTC
@@ -2,7 +2,7 @@
#
#====================================================================#
# #
-# Copyright 2002,2003,2004,2005,2006,2007,2008,2009 #
+# Copyright 2002-2014,2015 #
# Mikael Granvik, Jenni Virtanen, Karri Muinonen, Teemu Laakso, #
# Dagmara Oszkiewicz #
# #
@@ -28,9 +28,9 @@
# re-formats them for OpenOrb.
#
# Author: MG
-# Date: 2009-07-27
+# Date: 2015-06-15

-wget http://www.cfa.harvard.edu/iau/lists/ObsCodes.html
+wget http://www.minorplanetcenter.net/iau/lists/ObsCodes.html
sed -e '2d' ObsCodes.html | grep -v "<" > OBSCODE.dat
rm -f ObsCodes.html

=======================================
--- /trunk/gnuplot/vov_plot_car.gp Fri Sep 16 10:27:02 2011 UTC
+++ /trunk/gnuplot/vov_plot_car.gp Mon Jun 15 22:42:10 2015 UTC
@@ -1,6 +1,6 @@
#====================================================================#
# #
-# Copyright 2002,2003,2004,2005,2006,2007,2008,2009,2010,2011 #
+# Copyright 2002-2014,2015 #
# Mikael Granvik, Jenni Virtanen, Karri Muinonen, Teemu Laakso, #
# Dagmara Oszkiewicz #
# #
@@ -26,7 +26,7 @@
# ranges.
#
# Author: MG
-# Version: 2011-09-16
+# Version: 2015-06-15
#
reset
set terminal postscript enhanced portrait linewidth 1.0 8.0
@@ -49,28 +49,28 @@
set ylabel 'z [AU]'
plot 'vov_orbits.out' using 1:3 pt 7, \
'vov_nominal_orbit.out' using 1:3 with points pt 3 ps 3.0, \
-'vov_sampling_grid.out' using 1:3:9:10 with yerrorbars lt 1
+'vov_sampling_grid.out' using 1:3:7:10 with yerrorbars lt 1
set size 0.5,0.33
set origin 0.0,0.33
set xlabel 'x [AU]'
set ylabel 'dx/dt [AU/d]'
plot 'vov_orbits.out' using 1:4 pt 7, \
'vov_nominal_orbit.out' using 1:4 with points pt 3 ps 3.0, \
-'vov_sampling_grid.out' using 1:4:11:12 with yerrorbars lt 1
+'vov_sampling_grid.out' using 1:4:7:12 with yerrorbars lt 1
set size 0.5,0.33
set origin 0.5,0.33
set xlabel 'x [AU]'
set ylabel 'dy/dt [AU/d]'
plot 'vov_orbits.out' using 1:5 pt 7, \
'vov_nominal_orbit.out' using 1:5 with points pt 3 ps 3.0, \
-'vov_sampling_grid.out' using 1:5:13:14 with yerrorbars lt 1
+'vov_sampling_grid.out' using 1:5:7:14 with yerrorbars lt 1
set size 0.5,0.33
set origin 0.0,0.0
set xlabel 'x [AU]'
set ylabel 'dz/dt [AU/d]'
plot 'vov_orbits.out' using 1:6 pt 7, \
'vov_nominal_orbit.out' using 1:6 with points pt 3 ps 3.0, \
-'vov_sampling_grid.out' using 1:6:15:16 with yerrorbars lt 1
+'vov_sampling_grid.out' using 1:6:7:16 with yerrorbars lt 1
set size 0.5,0.33
set origin 0.5,0.0
set xlabel 'x [AU]'
=======================================
--- /trunk/main/io.f90 Wed Feb 25 08:09:08 2015 UTC
+++ /trunk/main/io.f90 Mon Jun 15 22:42:10 2015 UTC
@@ -1,6 +1,6 @@
!====================================================================!
! !
-! Copyright 2002-2013,2014 !
+! Copyright 2002-2014,2015 !
! Mikael Granvik, Jenni Virtanen, Karri Muinonen, Teemu Laakso, !
! Dagmara Oszkiewicz !
! !
@@ -27,7 +27,7 @@
!! called from main programs.
!!
!! @author MG, JV
-!! @version 2014-08-22
+!! @version 2015-06-15
!!
MODULE io

@@ -1278,7 +1278,7 @@
END IF
CASE ("sor.ran.obs")
IF (PRESENT(sor_random_obs)) THEN
-! sor_random_obs = .TRUE.
+ ! sor_random_obs = .TRUE.
READ(par_val, *, iostat=err) sor_random_obs
IF (err /= 0) THEN
error = .TRUE.
=======================================
--- /trunk/main/oorb.conf Thu Oct 13 09:04:18 2011 UTC
+++ /trunk/main/oorb.conf Mon Jun 15 22:42:10 2015 UTC
@@ -1,6 +1,6 @@
#====================================================================#
# #
-# Copyright 2002,2003,2004,2005,2006,2007,2008,2009,2010,2011 #
+# Copyright 2002-2014,2015 #
# Mikael Granvik, Jenni Virtanen, Karri Muinonen, Teemu Laakso, #
# Dagmara Oszkiewicz #
# #
@@ -27,7 +27,7 @@
# lines starting with "#" are ignored.
#
# Author: MG
-# Version: 2011-10-13
+# Version: 2015-06-15


# GENERAL SETTINGS
@@ -39,9 +39,10 @@
verbose.error: 1

# Name of binary file containing JPL's planetary ephemerides,
-# typically deXXX.dat, where XXX is either 405 or 406. Note that the
-# name used *must* include either 405 or 406 since there is no other
-# way to differentiate between the contents.
+# typically deXXX.dat, where XXX is either 405 or 406, or IMCCE's
+# INPOP planetary ephemerides that can be read with code adhering to
+# the testeph.f algorithm such as the little-endian TDB version of
+# INPOP10b.
planetary_ephemeris_fname: de405.dat


@@ -64,7 +65,7 @@

# Element type to be used in the resulting orbital element file
# [ keplerian | cartesian | cometary ]
-element_type_out: keplerian
+element_type_out: cartesian

# Automatic plot generation (T=yes/F=no)
plot.results: T
@@ -95,7 +96,7 @@
#epoch.mjd: 53800.0

# Toggle outlier rejection
-#outlier_rejection:
+outlier_rejection: T

# Outlier criterion (sigma multiplier)
outlier.multiplier: 4.0
@@ -136,7 +137,7 @@
integrator: bulirsch-stoer

# Integrator step length (in days)
-integration_step: 2.0
+integration_step: 1.0

# Relativistic corrections
relativity: T
@@ -151,7 +152,7 @@
integrator_init: bulirsch-stoer

# Integrator step length of the initial orbit (in days)
-integration_step_init: 2.0
+integration_step_init: 1.0


# STATISTICAL PARAMETERS
@@ -203,7 +204,7 @@
#apriori.rho.max :


-# STATISTICAL ORBITAL RANGING
+# MC / MCMC STATISTICAL ORBITAL RANGING

# Method for solving the two-point boundary-value problem
# [ continued fraction | p-iteration | n-body amoeba ]
@@ -216,7 +217,7 @@
sor.norb: 2000

# Maximum number of trial orbits
-sor.ntrial: 10000000
+sor.ntrial: 1000000

# Method for solving the two-point boundary-value problem during
# the preliminary steps of stepwise Ranging
@@ -416,7 +417,10 @@

# PHYSICAL PARAMETERS

-# Toggle estimation of H(alpha=0) magnitude
+# Toggle estimation of H(alpha=0) magnitude. Turning on this option
+# requires that the astrometry data file contain useful magnitude
+# information. This generally means it has to have at least three
+# values.
# [ T | F ]
pp.H_estimation: T

=======================================
--- /trunk/main/oorb.f90 Fri Apr 17 07:48:32 2015 UTC
+++ /trunk/main/oorb.f90 Mon Jun 15 22:42:10 2015 UTC
@@ -26,7 +26,7 @@
!! Main program for various tasks that include orbit computation.
!!
!! @author MG
-!! @version 2015-03-11
+!! @version 2015-06-15
!!
PROGRAM oorb

@@ -877,7 +877,6 @@

frame = get_cl_option("--frame=", "ecliptic")

-
! Initialize random number generator using current working
! directory and system clock
IF (.NOT.get_cl_option("--fixran",.FALSE.)) THEN
@@ -1966,7 +1965,7 @@
str = TRIM(id) // "_"// TRIM(str)
! WRITE RANGING OUTPUT FILE WITH
CALL NEW(out_file, TRIM(id) // ".sor")
-! CALL NEW(out_file, TRIM(str) // ".sor")
+ ! CALL NEW(out_file, TRIM(str) // ".sor")
CALL OPEN(out_file)
IF (error) THEN
CALL errorMessage("oorb / ranging", &
@@ -2026,7 +2025,7 @@
END IF
IF (separately) THEN
CALL NEW(out_file, TRIM(id) // ".orb")
-! CALL NEW(out_file, TRIM(str) // ".orb")
+ ! CALL NEW(out_file, TRIM(str) // ".orb")
CALL OPEN(out_file)
IF (error) THEN
CALL errorMessage("oorb / ranging", &
@@ -2154,7 +2153,6 @@
STOP
END IF
IF (compress) THEN
-
CALL system("gzip -f " // TRIM(str)
// "_ranging_residual_stamps.eps")
END IF
IF (plot_open) THEN
@@ -2167,7 +2165,7 @@
STOP
END IF

- ! Topocentric ranges
+ ! Topocentric ranges
CALL NEW(tmp_file, "sor_histo.out")
CALL OPEN(tmp_file)
IF (error) THEN
@@ -2204,7 +2202,6 @@
"TRACE BACK (176)", 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))
@@ -2262,7 +2259,7 @@
pdf=pdf_arr_cmp, std_dev=stdev, errstr=errstr)
CASE ("mcmc")
CALL moments(elements_arr(:,j), &
- std_dev=stdev, errstr=errstr)
+ std_dev=stdev, errstr=errstr)
END SELECT
WRITE(getUnit(tmp_file), "(F22.15,1X)", &
advance="no") stdev
@@ -2303,8 +2300,8 @@
TRIM(element_type_comp_prm) // &
"_results.eps* &")
END IF
-! CALL system("cp " // TRIM(str) // &
-! "_ranging_ranges.out sor_ranges.out")
+ ! CALL system("cp " // TRIM(str) // &
+ ! "_ranging_ranges.out sor_ranges.out")
CALL system("gnuplot " // TRIM(gnuplot_scripts_dir)
// "/sor_plot_range.gp")
CALL system("cp sor_ranges.eps " // TRIM(str) // &
"_ranging_ranges.eps")
@@ -2326,7 +2323,6 @@
WRITE(stdout,"(3(1X,A))") "Object", &
TRIM(id), "successfully processed."
END IF
-
END IF
CALL NULLIFY(storb)
CALL NULLIFY(orb)
@@ -2339,7 +2335,6 @@



-
CASE ("simplex")

! Orbit optimization using the downhill simplex method.
@@ -7368,7 +7363,6 @@
mag = mags(j) - 5.0_bp*LOG10(SQRT(heliocentric_r2*ephemeris_r2))

! Filter transformations to V:
- !
IF (obsy_code_arr(j) == "F51") THEN

! PS1SC internal analysis by A. Fitzsimmons 08-11-2011.
=======================================
--- /trunk/modules/planetary_data.f90 Fri Oct 26 06:45:48 2012 UTC
+++ /trunk/modules/planetary_data.f90 Mon Jun 15 22:42:10 2015 UTC
@@ -1,6 +1,6 @@
!====================================================================!
! !
-! Copyright 2002-2011,2012 !
+! Copyright 2002-2014,2015 !
! Mikael Granvik, Jenni Virtanen, Karri Muinonen, Teemu Laakso, !
! Dagmara Oszkiewicz !
! !
@@ -49,7 +49,7 @@
!!</pre>
!!
!! @author MG, TL
-!! @version 2012-10-26
+!! @version 2015-06-15
!!
MODULE planetary_data

=======================================
--- /trunk/modules/statistics.f90 Wed Nov 30 07:53:19 2011 UTC
+++ /trunk/modules/statistics.f90 Mon Jun 15 22:42:10 2015 UTC
@@ -1,6 +1,6 @@
!====================================================================!
! !
-! Copyright 2002,2003,2004,2005,2006,2007,2008,2009,2010,2011 !
+! Copyright 2002-2014,2015 !
! Mikael Granvik, Jenni Virtanen, Karri Muinonen, Teemu Laakso, !
! Dagmara Oszkiewicz !
! !
@@ -26,7 +26,7 @@
!! Tools for statistics.
!!
!! @author MG
-!! @version 2011-11-30
+!! @version 2015-06-15
!!
MODULE statistics

@@ -203,7 +203,7 @@
REAL(rprec8), DIMENSION(:), INTENT(in), OPTIONAL :: pdf
LOGICAL, DIMENSION(:), OPTIONAL, INTENT(in) :: mask
REAL(rprec8), OPTIONAL, INTENT(out) :: mean, std_dev,
skew, kurt
- CHARACTER(len=*), INTENT(inout) :: errstr
+ CHARACTER(len=*), INTENT(inout) :: errstr

REAL(rprec8), DIMENSION(:), ALLOCATABLE :: pdf_
REAL(rprec8) :: mean_, std_dev_