diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..9e36487 --- /dev/null +++ b/.gitignore @@ -0,0 +1,18 @@ +.idea +__pycache__ +build +ModelUsed.dat +OSSOS-model.out +*.log +*.ipynb +*.so +*.o +*.mod +Classical/Classical +Commons/Driver +Commons/GiMeObj.f95 +Detached/Detached +Inner/Inner +Plutinos/Plutinos +Scattering/Scattering +Twotinos/Twotinos \ No newline at end of file diff --git a/Commons/Driver.f95 b/Commons/Driver.f95 new file mode 100644 index 0000000..495c15d --- /dev/null +++ b/Commons/Driver.f95 @@ -0,0 +1,181 @@ +program Driver +!-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-* +! +! File Driver.f95 +! +! J.-M. Petit Observatoire de Besac +! Version 1: May 2013 +! +! The purpose of this program is to compare models of the Kuiper Belt to +! reality by actually comparing what the surveys (say CFEPS or OSSOS) +! would have found if the model was a good representation of the real +! world, to the objects that were actually found. +! +! The workflow of the survey simulator driver is as follows: +! +! Loop (some condition): +! call GiMeObj(arg_list_1) +! Check model ended: +! set exit condition +! call Detos1(arg_list_2) +! Check detection and tracking: +! store results +! +! The GiMeObj routine is in charge of providing a single new object at +! each call. The Detos1 routine determines is the proposed object would +! have been detected by the survey. +! +! The survey simulator expects orbital elements with respect to ecliptic +! reference frame. +! +! Logical unit numbers 7 to 19 are reserved to use by Driver.f and +! SurveySubs.f and should not be used by GiMeObj or any other routine +! that you may add to the driver. Please use logical unit numbers +! starting from 20. +! +! As currently written, the driver includes a file 'GiMeObj.f' +! containing the definition of the model. The correct way to use this +! feature is to have one's GiMeObj routine in a file and +! create a symobolic link: +! +! ln -s GiMeObj.f +! +! For more information, please refer to README files in the source directory. +! +!-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-* + + use gimeobjut + + implicit none + + integer, parameter :: screen = 6, keybd = 5, verbose = 9 + integer :: lun_h + type(t_orb_m) :: o_m +! color array NEEDS to be length 10 or more! + real (kind=8) :: h, epoch, hmax, rn_iter + integer :: ierr, seed, n_iter, nmax, nchar, values(8), seed_mod + character(80) :: det_outfile, line + character(100) :: comments + character(10) :: time + character(8) :: date + character(5) :: zone + logical :: keep_going, rec + + lun_h = 10 + keep_going = .true. + +! Get arguments +! Seed for random number generator + read (keybd, *, err=9999) seed + seed_mod = seed +! -maximum number of trials (<0) or calibrated number of objects + read (keybd, *, err=9999) nmax +! Maximum value of H to draw + read (keybd, *, err=9999) hmax +! Do we want to record the objects inside GiMeObj ? + read (keybd, *, err=9999) rec +! Name for the model objects outfile + read (keybd, '(a)', err=9999) det_outfile + det_outfile = strip_comment(det_outfile) + +! print *, n_track_max + +! Open output files and write header + open (unit=lun_h, file=det_outfile, status='new', err=9500) + + call date_and_time(date, time, zone, values) + write (lun_h, '(a17,a23,2x,a5)') '# Creation time: ', & + date(1:4)//'-'//date(5:6)//'-'//date(7:8)//'T' & + //time(1:2)//':'//time(3:4)//':'//time(5:10), zone + write (lun_h, '(''#'')') + write (lun_h, '(''# Seed: '', i10)') seed + write (lun_h, '(''#'')') + write (lun_h, '(a,a)') & + '# a e i Omega omega M', & + ' H epoch comment ' + +! Initialize counters + n_iter = 0 + rn_iter = 0.d0 + +! Open and read in object distribution +100 continue + +! Main loop: loop on objects + if (keep_going) then + +! Select object. + nchar = 0 + call GiMeObj (seed_mod, nmax, hmax, rec, & + o_m, epoch, h, comments, nchar, ierr) + + if (ierr .eq. -10) then ! Something wrong with this object, go to next one + goto 100 + else if (ierr .eq. -20) then ! Something very wrong happend, stop + write (screen, '(a)') 'GiMeObj returned -20, stopping.' + goto 2010 + else if (ierr .eq. 100) then ! Reached end of model, prepare to stop + keep_going = .false. + end if + +! Count number of iterations; may exceed 2**31, so use a real*8 helper + n_iter = n_iter + 1 + if (n_iter .gt. 2000000000) then + rn_iter = rn_iter + dble(n_iter) + n_iter = 0 + if (rn_iter .ge. 9.d9) goto 2000 + end if + + write (lun_h, 9000) o_m%a, o_m%e, o_m%inc/drad, o_m%node/drad, & + o_m%peri/drad, o_m%m/drad, h, epoch, comments(1:nchar) + + goto 100 +! end of the if( keep_going ) loop + end if + +2000 continue + write (lun_h, '(''#'')') + write (lun_h, '(''# Total number of objects: '', f11.0)') & + rn_iter + dble(n_iter) + close (lun_h) + + call exit (0) + +9000 format (f9.4,1x,5(f8.4,1x),f6.2,1x,f13.5,1x,a) + +9500 continue + write (screen, *) 'File --- "', det_outfile, '" already exists. ' + goto 9502 + +9502 continue + write (screen, *) 'Make sure "', det_outfile, & + '" do not exist and restart ModelDriver.' + call exit(-1) + +9999 continue + write (screen, *) 'Usage: ModelDriver < input' + write (screen, *) + write (screen, *) 'This will read the following from the keyboard:' + write (screen, *) '' + write (screen, *) '' + write (screen, *) '' + write (screen, *) 'rec' + write (screen, *) '' + write (screen, *) 'where:' + write (screen, *) & + ': integer used as seed for the random number generator' + write (screen, *) ': -maximum number of trials if < 0' + write (screen, *) ' Calibrated number of object if >= 0' + write (screen, *) ': maximum value of H up to which to draw' + write (screen, *) & + ': do we want to record the model objects inside GiMeObj' + write (screen, *) ': output file for model objects' + + call exit(0) + + 2010 continue + write (lun_h,*) "# Failed while generating model object." + call exit(-20) + + +end program Driver diff --git a/Commons/Makefile b/Commons/Makefile new file mode 100644 index 0000000..a0722ca --- /dev/null +++ b/Commons/Makefile @@ -0,0 +1,70 @@ +U = datadec ioutils elemutils rot +OBJU = $(addsuffix .o,${U}) +SRCU = $(addsuffix .f95,${U}) +FPPU = $(addsuffix .fpp,${U}) +MODU = $(addsuffix .mod,${U}) + +M = modelutils GiMeObj +OBJM = $(addsuffix .o,${M}) +SRCM = $(addsuffix .f95,${M}) +FPPM = $(addsuffix .fpp,${M}) +MODM = $(addsuffix .mod,${M}) + + +GOBJ = GiMeObjF95 + +FC = gfortran +FPP = gfortran -E -x f95-cpp-input -fPIC -DPYTHON + +FFLAGS = -O3 -fPIC -x f95-cpp-input +FFLAGP = -O3 +CFLAGS = +CPPFLAGS = + +TOPTARGETS := clean test +GIMEOBJ ?= OSSOS-classical-model + +all: Driver + +%.fpp : %.f95 + $(FPP) $< -o $@ + +%.o : %.f95 + $(FC) $(FFLAGS) -c $< + +.PHONY: clean + +clean: rmtmp + rm -f _$(GOBJ)*.so Driver GiMeObj.f95 + +rmtmp: + rm -f *~ *.o *.mod *.fpp f90wrap*f90 + rm -f $(GOBJ).py + rm -f *.log LOG .f2py_f2cmap + +.PHONY: link + +GiMeObj.f95: link + +link: + \rm -f GiMeObj.f95 + \rm -f GiMeObj.o + ln -s $(GIMEOBJ).f95 GiMeObj.f95 + +Driver: link Driver.f95 $(OBJU) ${SRCM} Makefile + $(FC) $(FFLAGP) -o Driver $(OBJU) ${SRCM} Driver.f95 + +Modules: _$(GOBJ).so Makefile + echo "Modules have been built" + +_$(GOBJ).so: link $(FPPU) $(OBJU) $(GIMEOBJ).f95 Makefile + \rm -f _$(GOBJ).so + f90wrap -k kind_map -m $(GOBJ) $(FPPU) 1>f90wrap_GiMeObj.log 2>&1 + f2py-f90wrap --fcompiler=$(FC) -c -m _$(GOBJ) $(OBJU) f90wrap_*.f90 1>f2py_GiMeObj.log 2>&1 + \rm f90wrap_*.f90 + \rm *.fpp *.mod *.o + sed -e 's/import _$(GOBJ)/from . import _$(GOBJ)/' -i.bck $(GOBJ).py + rm $(GOBJ).py.bck + rm .f2py_f2cmap + mv _$(GOBJ)*.so _$(GOBJ).so + diff --git a/Commons/datadec.f95 b/Commons/datadec.f95 new file mode 100644 index 0000000..a8c0c32 --- /dev/null +++ b/Commons/datadec.f95 @@ -0,0 +1,36 @@ +module datadec + + ! define length of array parameters + integer, parameter :: nw_max = 20, nparmax = 20 + + ! define some useful constants + real (kind=8), parameter :: Pi = 3.141592653589793238d0, drad = Pi/180.0D0, & + TwoHours = 2.d0/24.d0, TwoPi = 2.0d0*Pi, eps = 1.d-14, & + Flatten = 0.003352813d0, & ! flattening of earth, 1/298.257 + Equat_Rad = 6378137.0d0 ! equatorial radius of earth, meters + + real (kind=8), parameter :: km2AU = 149597870.700d0 + + + real (kind=8), parameter :: gmb = 1.d0+1.d0/6023600.0d0+1.d0/408523.71d0 & + +1.d0/328900.56d0+1.d0/3098708.0d0+1.d0/1047.3486d0+1.d0/3497.898d0 & + +1.d0/22902.98d0+1.d0/19412.24d0+1.d0/1.35d8 + + ! Internal variables + real (kind=8) :: om_lim_low, om_lim_high + common /om_lim_com/ om_lim_low, om_lim_high + + ! define data type to represent survey efficiency and pointings, and objects + type t_orb_m + real (kind=8) :: a, e, inc, node, peri, m + end type t_orb_m + + type t_orb_p + real (kind=8) :: a, e, inc, node, peri, tperi + end type t_orb_p + + type t_v3d + real (kind=8) :: x, y, z + end type t_v3d + +end module datadec diff --git a/Commons/elemutils.f95 b/Commons/elemutils.f95 new file mode 100644 index 0000000..1cbb0f7 --- /dev/null +++ b/Commons/elemutils.f95 @@ -0,0 +1,437 @@ +module elemutils + + use datadec + +contains +! \subroutine{coord\_cart} + + subroutine coord_cart (mu, o_m, p, v) + +! This routine transforms delaunay variables into cartisian +! variables. +! +! ANGLES ARE GIVEN IN RADIAN !!!! +! +! \subsection{Arguments} +! \subsubsection{Definitions} +! \begin{verse} +! \verb|a, e, inc, node, peri, m| = osculating elements \\ +! \verb|x, y, z, vx, vy, vz| = cartesian variables: $X, Y, Z, Px, Py, Pz$ +! \end{verse} +! +! Angles are in [radian] +! +! \subsubsection{Declarations} +! +!f2py intent(in) mu +!f2py intent(in) o_m +!f2py intent(out) p +!f2py intent(out) v + implicit none + + type(t_orb_m), intent(in) :: o_m + type(t_v3d), intent(out) :: p, v + real (kind=8), intent(in) :: mu + +! \subsection{Variables} +! \subsubsection{Definitions} +! \begin{verse} +! \verb|cos_e, sin_e, cos_i, sin_i| = sinus and cosines of $E$ and +! $i$ \\ +! \verb|delau| = Delaunay variables: $l, \cos(g), \sin(g), \cos(h), +! \sin(h), L, G, H$ \\ +! \verb|e| = eccentric anomaly \\ +! \verb|mat| = rotation matrix \\ +! \verb|q_vec, qp| = $q$ and $\dot q$ \\ +! \verb|tmp| = temporary variable +! \end{verse} +! +! \subsubsection{Declarations} + integer :: i + real (kind=8) :: delau(8), cos_e, cos_i, e, mat(3,3), q_vec(2), qp(2), & + sin_e, sin_i, tmp, signe, f, de, fp, fpp, fppp + +! Computation of sinus and cosines of angles. + + signe = 1.d0 + if (o_m%a .lt. 0.d0) signe = -1.d0 + cos_i = dcos(o_m%inc) + sin_i = dsqrt(1.d0 - cos_i**2) + delau(2) = dcos(o_m%peri) + delau(3) = dsin(o_m%peri) + delau(4) = dcos(o_m%node) + delau(5) = dsin(o_m%node) + delau(1) = o_m%m - int(o_m%m/TwoPi)*TwoPi + delau(6) = signe*dsqrt(mu*o_m%a*signe) + delau(7) = abs(delau(6))*dsqrt((1.d0 - o_m%e**2)*signe) + +! Rotation matrix. +! The rotation matrix is the composition of 3 matrices: $R_{xq} = +! R_3(-h) \cdot R_1(-i) \cdot R_3(-g)$: +! \begin{displaymath} +! R_{xq} = \left(\matrix{ +! \cos(h)\cos(g)-\frac{H}{G}\sin(h)\sin(g)& +! -\cos(h)\sin(g)-\frac{H}{G}\sin(h)\cos(g)& +! \sqrt{1-\frac{H^2}{G^2}}\sin(h) \cr +! \sin(h)\cos(g)+\frac{H}{G}\cos(h)\sin(g)& +! -\sin(h)\sin(g)+\frac{H}{G}\cos(h)\cos(g)& +! -\sqrt{1-\frac{H^2}{G^2}}\cos(h) \cr +! \sqrt{1-\frac{H^2}{G^2}}\sin(g)& +! \sqrt{1-\frac{H^2}{G^2}}\cos(g)& +! \frac{H}{G} \cr}\right); +! \end{displaymath} + + mat(1,1) = delau(4)*delau(2) - cos_i*delau(5)*delau(3) + mat(1,2) = -delau(4)*delau(3) - cos_i*delau(5)*delau(2) + mat(2,1) = delau(5)*delau(2) + cos_i*delau(4)*delau(3) + mat(2,2) = -delau(5)*delau(3) + cos_i*delau(4)*delau(2) + mat(3,1) = sin_i*delau(3) + mat(3,2) = sin_i*delau(2) + +! Eccentric anomaly. +! We solve iteratively the equation: +! \begin{displaymath} +! E - e \sin(E) = l +! \end{displaymath} +! using the accelerated Newton's method (see Danby). + + e = delau(1) + sign(.85d0, dsin(delau(1)))*o_m%e + i = 0 +1000 continue + sin_e = o_m%e*dsin(e) + f = e - sin_e - delau(1) + if (dabs(f) .gt. 1.d-14) then + cos_e = o_m%e*dcos(e) + fp = 1.d0 - cos_e + fpp = sin_e + fppp = cos_e + de = -f/fp + de = -f/(fp + de*fpp/2.d0) + de = -f/(fp + de*fpp/2.d0 + de*de*fppp/6.d0) + e = e + de + i = i + 1 + if (i .lt. 20) goto 1000 + write (6, *) 'COORD_CART: No convergence: ', i, f + write (6, *) mu, e, de + write (6, *) o_m%a, o_m%e, o_m%inc + write (6, *) o_m%node, o_m%peri, o_m%m + stop + end if +1100 continue + +! Coordinates relative to the orbit. +! The cartisian coordinate are given by $\vec X = R_{xq} \vec q$ +! and $\vec P = R_{xq} \dot{\vec q}$, where: +! \begin{eqnarray*} +! \vec q & = & \left(\frac{L^2}{\mu}(\cos(E) - e), +! \frac{GL}{\mu}\sin(E), 0\right), \\ +! \dot{\vec q} & = & \frac{\mu}{L(1 - e\cos(E))} +! \left(-\sin(E), \frac{G}{L}\cos(E), 0\right) +! \end{eqnarray*} + + cos_e = dcos(e) + sin_e = dsin(e) + q_vec(1) = delau(6)**2*(cos_e - o_m%e)/mu + q_vec(2) = delau(7)*delau(6)*sin_e/mu + tmp = mu/(delau(6)*(1.d0 - o_m%e*cos_e)) + qp(1) = -sin_e*tmp + qp(2) = delau(7)*cos_e*tmp/delau(6) + +! Cartisian coordinates + + p%x = mat(1,1)*q_vec(1) + mat(1,2)*q_vec(2) + p%y = mat(2,1)*q_vec(1) + mat(2,2)*q_vec(2) + p%z = mat(3,1)*q_vec(1) + mat(3,2)*q_vec(2) + v%x = mat(1,1)*qp(1) + mat(1,2)*qp(2) + v%y = mat(2,1)*qp(1) + mat(2,2)*qp(2) + v%z = mat(3,1)*qp(1) + mat(3,2)*qp(2) + + return + end subroutine coord_cart + +! \subroutine{pos\_cart} + + subroutine pos_cart (o_m, p) + +! This routine transforms delaunay variables into cartisian +! variables, positions only. +! +! \subsection{Arguments} +! \subsubsection{Definitions} +! \begin{verse} +! \verb|a, e, inc, node, peri, m| = osculating elements \\ +! \verb|x, y, z| = cartesian variables: $X, Y, Z$ +! \end{verse} +! +! Angles are in [radian] +! +! \subsubsection{Declarations} +! +!f2py intent(in) o_m +!f2py intent(out) p + implicit none + + type(t_orb_m), intent(in) :: o_m + type(t_v3d), intent(out) :: p + +! \subsection{Variables} +! \subsubsection{Definitions} +! \begin{verse} +! \verb|cos_e, sin_e, cos_i, sin_i| = sinus and cosines of $E$ and +! $i$ \\ +! \verb|delau| = Delaunay variables: $l, \cos(g), \sin(g), \cos(h), +! \sin(h), L, G, H$ \\ +! \verb|e| = eccentric anomaly \\ +! \verb|mat| = rotation matrix \\ +! \verb|q_vec| = $q$ \\ +! \end{verse} +! +! \subsubsection{Declarations} + integer :: i + real (kind=8) :: delau(8), cos_e, cos_i, e, mat(3,3), q_vec(2), & + sin_e, sin_i, signe, f, de, fp, fpp, fppp + +! Computation of sinus and cosines of angles. + + signe = 1.d0 + if (o_m%a .lt. 0.d0) signe = -1.d0 + cos_i = dcos(o_m%inc) + sin_i = dsqrt(1.d0 - cos_i**2) + delau(2) = dcos(o_m%peri) + delau(3) = dsin(o_m%peri) + delau(4) = dcos(o_m%node) + delau(5) = dsin(o_m%node) + delau(1) = o_m%m - int(o_m%m/TwoPi)*TwoPi + delau(6) = signe*dsqrt(o_m%a*signe) + delau(7) = abs(delau(6))*dsqrt((1.d0 - o_m%e**2)*signe) + +! Rotation matrix. +! The rotation matrix is the composition of 3 matrices: $R_{xq} = +! R_3(-h) \cdot R_1(-i) \cdot R_3(-g)$: +! \begin{displaymath} +! R_{xq} = \left(\matrix{ +! \cos(h)\cos(g)-\frac{H}{G}\sin(h)\sin(g)& +! -\cos(h)\sin(g)-\frac{H}{G}\sin(h)\cos(g)& +! \sqrt{1-\frac{H^2}{G^2}}\sin(h) \cr +! \sin(h)\cos(g)+\frac{H}{G}\cos(h)\sin(g)& +! -\sin(h)\sin(g)+\frac{H}{G}\cos(h)\cos(g)& +! -\sqrt{1-\frac{H^2}{G^2}}\cos(h) \cr +! \sqrt{1-\frac{H^2}{G^2}}\sin(g)& +! \sqrt{1-\frac{H^2}{G^2}}\cos(g)& +! \frac{H}{G} \cr}\right); +! \end{displaymath} + + mat(1,1) = delau(4)*delau(2) - cos_i*delau(5)*delau(3) + mat(1,2) = -delau(4)*delau(3) - cos_i*delau(5)*delau(2) + mat(2,1) = delau(5)*delau(2) + cos_i*delau(4)*delau(3) + mat(2,2) = -delau(5)*delau(3) + cos_i*delau(4)*delau(2) + mat(3,1) = sin_i*delau(3) + mat(3,2) = sin_i*delau(2) + +! Eccentric anomaly. +! We solve iteratively the equation: +! \begin{displaymath} +! E - e \sin(E) = l +! \end{displaymath} +! using the accelerated Newton's method (see Danby). + + e = delau(1) + sign(.85d0, dsin(delau(1)))*o_m%e + i = 0 +1000 continue + sin_e = o_m%e*dsin(e) + f = e - sin_e - delau(1) + if (dabs(f) .gt. 1.d-14) then + cos_e = o_m%e*dcos(e) + fp = 1.d0 - cos_e + fpp = sin_e + fppp = cos_e + de = -f/fp + de = -f/(fp + de*fpp/2.d0) + de = -f/(fp + de*fpp/2.d0 + de*de*fppp/6.d0) + e = e + de + i = i + 1 + if (i .lt. 20) goto 1000 + write (6, *) 'POS_CART: No convergence: ', i, f + write (6, *) e, de + write (6, *) o_m%a, o_m%e, o_m%inc + write (6, *) o_m%node, o_m%peri, o_m%m + stop + end if +1100 continue + +! Coordinates relative to the orbit. +! The cartisian coordinate are given by $\vec X = R_{xq} \vec q$ +! and $\vec P = R_{xq} \dot{\vec q}$, where: +! \begin{eqnarray*} +! \vec q & = & \left(\frac{L^2}{\mu}(\cos(E) - e), +! \frac{GL}{\mu}\sin(E), 0\right), \\ +! \dot{\vec q} & = & \frac{\mu}{L(1 - e\cos(E))} +! \left(-\sin(E), \frac{G}{L}\cos(E), 0\right) +! \end{eqnarray*} + + cos_e = dcos(e) + sin_e = dsin(e) + q_vec(1) = delau(6)**2*(cos_e - o_m%e) + q_vec(2) = delau(7)*delau(6)*sin_e + +! Cartisian coordinates + + p%x = mat(1,1)*q_vec(1) + mat(1,2)*q_vec(2) + p%y = mat(2,1)*q_vec(1) + mat(2,2)*q_vec(2) + p%z = mat(3,1)*q_vec(1) + mat(3,2)*q_vec(2) + + return + end subroutine pos_cart + +! \subroutine{osc\_el} + + subroutine osc_el (mu, p, v, o_m) + +! This routine transforms cartisian variables into delaunay +! variables. +! +! \subsection{Arguments} +! \subsubsection{Definitions} +! \begin{verse} +! \verb|a, e, inc, node, peri, m| = osculating elements \\ +! \verb|cart| = cartesian variables: $X, Y, Z, Px, Py, Pz$ +! \end{verse} +! +! \subsubsection{Declarations} +!f2py intent(in) mu +!f2py intent(in) p +!f2py intent(in) v +!f2py intent(out) o_m + implicit none + + type(t_v3d), intent(in) :: p, v + type(t_orb_m), intent(out) :: o_m + real (kind=8), intent(in) :: mu + +! \subsection{Variables} +! \subsubsection{Definitions} +! \begin{verse} +! \verb|cos_i, sin_i| = sinus and cosines of $i$ \\ +! \verb|delau| = Delaunay variables: $l, \cos(g), \sin(g), \cos(h), +! \sin(h), L, G, H$ \\ +! \verb|e| = eccentric anomaly \\ +! \verb|f| = $f$ true anomaly \\ +! \verb|g| = $g$ argument of pericenter \\ +! \verb|h_vec| = $\vec h = \vec X \times \vec P$ \\ +! \verb|p_vec| = $\vec p = -\mu \frac{\vec X}{r} +! - \vec h \times \vec P$ \\ +! \verb|r| = radial distance \\ +! \verb|tmp1, tmp2| = temporary variables \\ +! \verb|v2| = velocity squared +! \end{verse} +! +! \subsubsection{Declarations} + real (kind=8) :: delau(8), e, f, h_vec(3), p_vec(3), & + r, tmp1, tmp2, v2, signe + +! Computation of angular momentum and eccentricity vector. +! \begin{eqnarray*} +! \vec h & = & \vec X \times \vec P, \\ +! \vec p & = & -\mu \frac{\vec X}{|\vec X|} +! - \vec h \times \vec P +! \end{eqnarray*} + + h_vec(1) = p%y*v%z - p%z*v%y + h_vec(2) = p%z*v%x - p%x*v%z + h_vec(3) = p%x*v%y - p%y*v%x + r = 1.d0/dsqrt(p%x**2 + p%y**2 + p%z**2) + p_vec(1) = -mu*p%x*r - h_vec(2)*v%z + h_vec(3)*v%y + p_vec(2) = -mu*p%y*r - h_vec(3)*v%x + h_vec(1)*v%z + p_vec(3) = -mu*p%z*r - h_vec(1)*v%y + h_vec(2)*v%x + +! Computation of momenta. +! \begin{eqnarray*} +! L & = & \mu\sqrt{\frac{1} +! {\frac{2\mu}{|\vec X|}-|\vec P|^2}}, \\ +! G & = & |\vec h|, \\ +! H & = & h_z +! \end{eqnarray*} + + v2 = v%x**2 + v%y**2 + v%z**2 + tmp1 = 2.d0*mu*r - v2 + signe = 1.d0 + if (tmp1 .lt. 0.d0) signe = -1.d0 + delau(6) = signe*mu/dsqrt(signe*tmp1) + delau(7) = h_vec(1)**2 + h_vec(2)**2 + h_vec(3)**2 + delau(8) = h_vec(3) + o_m%e = dsqrt(dmax1(mu**2 - delau(7)*tmp1, 0.d0))/mu + delau(7) = dsqrt(delau(7)) + + if ((p%z .eq. 0.d0) .and. (v%z .eq. 0.d0)) then + delau(4) = 1.d0 + delau(5) = 0.d0 + tmp1 = 1.d0/dsqrt(p_vec(1)**2 + p_vec(2)**2) + delau(2) = p_vec(1)*tmp1 + delau(3) = p_vec(2)*tmp1 + else + +! Longitude of node. +! \begin{eqnarray*} +! \cos(h) & = & -\frac{h_y}{\sqrt{h_x^2 + h_y^2}}, \\ +! \sin(h) & = & \frac{h_x}{\sqrt{h_x^2 + h_y^2}} +! \end{eqnarray*} + + tmp1 = 1.d0/dsqrt(h_vec(1)**2 + h_vec(2)**2) + delau(4) = -h_vec(2)*tmp1 + delau(5) = h_vec(1)*tmp1 + +! Argument of pericenter. +! Let us call $\vec N$ the vector derived from $\vec h$, pointing +! to the ascending node: +! \begin{displaymath} +! \vec N = \left(-h_y, h_x, 0\right) +! \end{displaymath} +! From this, we get: +! \begin{eqnarray*} +! \cos(g) & = & \frac{\vec N \cdot \vec p}{|\vec N||\vec p|}, \\ +! \sin(g) & = & \frac{\vec N \times \vec p}{|\vec N||\vec p|} +! \cdot \frac{\vec h}{|\vec h|} +! \end{eqnarray*} + + tmp2 = 1.d0/dsqrt(p_vec(1)**2 + p_vec(2)**2 + p_vec(3)**2) + delau(2) = (h_vec(1)*p_vec(2) - h_vec(2)*p_vec(1))*tmp1*tmp2 + delau(3) = ((h_vec(1)**2 + h_vec(2)**2)*p_vec(3) & + - h_vec(3)*(h_vec(1)*p_vec(1) + h_vec(2)*p_vec(2))) & + *tmp1*tmp2/delau(7) + end if + +! Mean anomaly +! We define $\vec X_{orb} = R_1(i) \cdot R_3(h) \vec X$. It turns +! out that $\vec X_{orb} = (r \cos(g+f), r \sin(g+f), 0)$. Hence: +! \begin{eqnarray*} +! \cos(g+f) & = & \cos(h) X + \sin(h) Y, \\ +! \sin(g+f) & = & \cos(i) \left(\cos(h) Y - \sin(h) X\right) +! + \sin(i) Z +! \end{eqnarray*} +! Furthermore, we have the relation: +! \begin{displaymath} +! \tan(\frac{E}{2}) = \sqrt{\frac{1 - e}{1 + e}} +! \tan(\frac{f}{2}) +! \end{displaymath} +! and finally: +! \begin{displaymath} +! l = E - e \sin(E) +! \end{displaymath} + + tmp1 = (p%x*p_vec(1) + p%y*p_vec(2) + p%z*p_vec(3)) + tmp2 = ((p%z*p_vec(2) - p%y*p_vec(3))*h_vec(1) & + + (p%x*p_vec(3) - p%z*p_vec(1))*h_vec(2) & + + (p%y*p_vec(1) - p%x*p_vec(2))*h_vec(3)) & + /delau(7) + f = datan2(tmp2, tmp1) + e = 2.d0*datan(dsqrt(signe*(1.d0 - o_m%e)/(1.d0 + o_m%e))*dtan(f/2.d0)) + o_m%m = e - o_m%e*dsin(e) + o_m%node = datan2(delau(5), delau(4)) + o_m%peri = datan2(delau(3), delau(2)) + o_m%a = signe*delau(6)**2/mu + o_m%inc = dacos(dmax1(dmin1(delau(8)/delau(7),1.d0),-1.d0)) + + return + end subroutine osc_el + +end module elemutils diff --git a/Commons/ioutils.f95 b/Commons/ioutils.f95 new file mode 100644 index 0000000..1d568e5 --- /dev/null +++ b/Commons/ioutils.f95 @@ -0,0 +1,141 @@ +module ioutils + + use datadec + +contains + + subroutine trim (base_name, i1, i2, finished, len) +!-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- +! This routine trims the input character string 'base_name' and returns +! the +! +! Author: Jean-Marc Petit +! version 1: December 2019 +!-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- +! INPUT +! base_name: Input string with file name in it (CH) +! len : Length of input string (I4) +! +! OUPUT +! i1 : Index of first character of file name (I4) +! i2 : Index of last character of file name (I4) +! finished: Whether it reached the end of the string when searching +! for first character +!-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- +!f2py intent(in) base_name +!f2py intent(out) i1 +!f2py intent(out) i2 +!f2py intent(out) finished +!f2py intent(in) len + implicit none + + integer, intent(in) :: len + integer, intent(out) :: i1, i2 + character(*), intent(in) :: base_name + logical, intent(out) :: finished + + finished = .false. + i1 = 1 +100 continue + if ((base_name(i1:i1) .eq. char(0)) & + .or. (base_name(i1:i1) .eq. char(9)) & + .or. (base_name(i1:i1) .eq. ' ')) then + i1 = i1 + 1 + if (i1 .eq. len) then + finished = .true. + return + end if + goto 100 + end if +101 continue + + i2 = len +110 continue + if ((base_name(i2:i2) .eq. char(0)) & + .or. (base_name(i2:i2) .eq. char(9)) & + .or. (base_name(i2:i2) .eq. ' ')) then + if (i2 .eq. i1) goto 111 + i2 = i2 - 1 + goto 110 + end if +111 continue + + return + end subroutine trim + + subroutine read_file_name (base_name, i1, i2, finished, len) +!-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- +! This routine trims the input character string 'base_name' and returns +! the +! +! Author: Jean-Marc Petit +! version 1: September 2017 +!-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- +! INPUT +! base_name: Input string with file name in it (CH) +! len : Length of input string (I4) +! +! OUPUT +! i1 : Index of first character of file name (I4) +! i2 : Index of last character of file name (I4) +! finished: Whether it reached the end of the string when searching +! for first character +!-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- +!f2py intent(in) base_name +!f2py intent(out) i1 +!f2py intent(out) i2 +!f2py intent(out) finished +!f2py intent(in) len + implicit none + + integer, intent(in) :: len + integer, intent(out) :: i1, i2 + character(*), intent(in) :: base_name + logical,intent(out) :: finished + + finished = .false. + i1 = 1 +100 continue + if ((base_name(i1:i1) .eq. char(0)) & + .or. (base_name(i1:i1) .eq. char(9)) & + .or. (base_name(i1:i1) .eq. ' ')) then + i1 = i1 + 1 + if (i1 .eq. len) then + finished = .true. + return + end if + goto 100 + end if +101 continue + + i2 = i1 + 1 +110 continue + if ((base_name(i2:i2) .ne. char(0)) & + .and. (base_name(i2:i2) .ne. char(9)) & + .and. (base_name(i2:i2) .ne. ' ')) then + if (i2 .eq. len) goto 111 + i2 = i2 + 1 + goto 110 + end if + i2 = i2 - 1 +111 continue + + return + end subroutine read_file_name + + character(100) function strip_comment(str) + character (*), intent(in) :: str + integer c_idx, i1, i2 + logical finished + strip_comment = str + c_idx = index(strip_comment, '!') + if (c_idx /= 0) then + strip_comment = strip_comment(1:c_idx-1) + end if + call trim(strip_comment, i1, i2, finished, len(strip_comment)) + strip_comment = strip_comment(i1:i2) + + end function strip_comment + +end module ioutils + diff --git a/Commons/kind_map b/Commons/kind_map new file mode 100755 index 0000000..1724b5b --- /dev/null +++ b/Commons/kind_map @@ -0,0 +1,18 @@ +{ + 'real': { '' : 'float', + '4' : 'float', + 'isp' : 'float', + '8' : 'double', + 'dp' : 'double', + 'idp' : 'double', + 'sng_ad' : 'float', + 'dbl_ad' : 'double'}, + 'complex' : { '' : 'complex_float', + '8' : 'complex_double', + '16' : 'complex_long_double', + 'dp' : 'complex_double'}, + 'integer' : { '2' : 'int', + '4' : 'int', + '8' : 'long_long', + 'dp' : 'long_long' }, +} diff --git a/Commons/modelutils.f95 b/Commons/modelutils.f95 new file mode 100644 index 0000000..2c756a7 --- /dev/null +++ b/Commons/modelutils.f95 @@ -0,0 +1,654 @@ +module modelutils + + use datadec + + interface + function func(nparam, param, inc) + real (kind=8) :: func + integer (kind=4), intent(in) :: nparam + real (kind=8), intent(in) :: param(*), inc + end function func + end interface + + type func_holder + procedure(func), pointer, nopass :: f_ptr => null() + end type func_holder + +contains + subroutine incdism (seed, nparam, param, incmin, incmax, inc, & + dist, ierr, func) +!-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- +! This routine draws randomly a variable according to probability +! density \verb|func| with parameters \verb|param|. Same as previous +! routine, but can remember up to 10 different distributions at a time. +!-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- +! +! J-M. Petit Observatoire de Besancon +! Version 1 : October 2006 +! +!-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- +! INPUT +! seed : Random number generator seed (I4) +! nparam: Number of parameters (I4) +! param : Parameters (n*R8) +! incmin: Minimum inclination (R8) +! incmax: Maximum inclination (R8) +! dist : index of the selected distribution (I4) +! func : probability density function +! +! OUTPUT +! inc : Inclination (R8) +! ierr : Error code +! 0 : nominal run +! 10 : wrong input data +!-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- +!f2py intent(in,out) seed +!f2py intent(in) nparam +!f2py intent(in), depend(nparam) :: param +!f2py intent(in) incmin +!f2py intent(in) incmax +!f2py intent(in) dist +!f2py intent(out) inc +!f2py intent(out) ierr + implicit none + + integer (kind=4), intent(in) :: nparam, dist + integer (kind=4), intent(inout) :: seed + integer (kind=4), intent(out) :: ierr + real (kind=8), intent(in) :: param(nparmax), incmin, incmax + real (kind=8), intent(out) :: inc + type(func_holder), intent(in) :: func + integer :: i, ilo, ihi, di + real (kind=8) :: random + integer (kind=4), parameter :: np = 16384, nd = 10 + real (kind=8), save :: proba(0:np,nd), inctab(0:np,nd) + logical, save :: first(nd) + + data first /.true.,.true.,.true.,.true.,.true., & + .true.,.true.,.true.,.true.,.true./ + + ierr = 0 + di = min(nd, dist) + if (first(di)) then + inctab(0,di) = incmin + proba(0,di) = 0.d0 + do i = 1, np + inctab(i,di) = incmin + dfloat(i)*(incmax-incmin)/dfloat(np) + proba(i,di) = func%f_ptr(nparam, param(1:nparam), inctab(i,di)) + proba(i-1,di) + end do + do i = 1, np + proba(i,di) = proba(i,di)/proba(np,di) + end do + first(di) = .false. + end if + + random = ran_3(seed) + inc = interp(proba(0,di), inctab(0,di), random, np+1) + + return + end subroutine incdism + + real (kind=8) function Variably_tapered(h, params) +!-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- +! This routine returns the number of objects brighter or equal to H +! following an exponentially tapered exponential with parameters in +! params. +!-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- +! +! J-M. Petit Observatoire de Besancon +! Version 1 : October 2021 +! +!-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- +! INPUT +! seed : Random number generator seed (I4) +! params: parameters for the distribution (4*R8) +! +! OUTPUT +! Variably_tapered : Random value of H (R8) +!-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- +!f2py intent(in) h +!f2py intent(in) params + implicit none + + real (kind=8), intent(in) :: params(4), h + + Variably_tapered = 10.d0**(params(3)*3.d0*(h-params(1))/5.d0) & + *exp(-10.d0**(-params(4)*3.d0*(h-params(2))/5.d0)) + + return + end function Variably_tapered + + real (kind=8) function Variably_tapered_diff(h, params) +!-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- +! This routine returns the differential number of objects at mag H per unit +! mag following an exponentially tapered exponential with parameters in +! params. +!-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- +! +! J-M. Petit Observatoire de Besancon +! Version 1 : August 2023 +! +!-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- +! INPUT +! seed : Random number generator seed (I4) +! params: parameters for the distribution (4*R8) +! +! OUTPUT +! Variably_tapered_diff : Random value of H (R8) +!-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- +!f2py intent(in) h +!f2py intent(in) params + implicit none + + real (kind=8), intent(in) :: params(4), h + + Variably_tapered_diff = (log(10.0d0)*3.0d0/5.0d0) & + *10.0d0**(params(3)*3.0d0*(h-params(1))/5.0d0) & + *(params(3) + params(4)*10.0d0**(-params(4)*3.0d0*(h-params(2))/5.0d0)) & + *exp(-10.0d0**(-params(4)*3.0d0*(h-params(2))/5.0d0)) + + return + end function Variably_tapered_diff + + real (kind=8) function offgau (nparam, param, inc) +!-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- +! This routine returns the unnormalized inclination "probability" +! density as a non-zero centered gaussian. +!-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- +! +! J-M. Petit Observatoire de Besancon +! Version 1 : August 2014 +! +!-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- +! INPUT +! nparam: Number of parameters (I4) +! param : Parameters (n*R8) +! inc : Inclination [rad] (R8) +! +! OUPUT +! offgau: Value of the probability (R8) +!-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- +!f2py intent(in) nparam +!f2py intent(in), depend(nparam) :: param +!f2py intent(in) inc + implicit none + + integer (kind=4), intent(in) :: nparam + real (kind=8), intent(in) :: param(*), inc + real (kind=8), parameter :: Pi = 3.141592653589793238d0, TwoPi = 2.0d0*Pi + real (kind=8) :: fe, s1, angle, t0, t1, t3 + + if (nparam .ne. 2) stop + t0 = param(1) + s1 = param(2) + t1 = 2.*s1**2 + angle = mod(inc, TwoPi) + if (angle .gt. Pi) angle = angle - TwoPi + t3 = -(t0-angle)**2 + if (t3 .lt. -300.d0*t1) then + fe = 0.d0 + else + fe = exp(t3/t1) + end if + offgau = dsin(angle)*fe + + return + end function offgau + + real (kind=8) function cold_inc_sq (nparam, param, inc) +!-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- +! This routine returns the unnormalized inclination "probability" +! density for cold objects. +!-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- +! +! J-M. Petit Observatoire de Besancon +! Version 1 : April 2023, 7th +! +!-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- +! INPUT +! nparam: Number of parameters (I4) +! param : Parameters (n*R8) +! inc : Inclination [rad] (R8) +! +! OUPUT +! cold_inc_sq: Value of the probability (R8) +!-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- +!f2py intent(in) nparam +!f2py intent(in), depend(nparam) :: param +!f2py intent(in) inc + implicit none + + integer (kind=4), intent(in) :: nparam + real (kind=8), intent(in) :: param(*), inc + real (kind=8), parameter :: Pi = 3.141592653589793238d0, & + TwoPi = 2.0d0*Pi, drad = Pi/180.0d0 + real (kind=8), save :: angle, a0, p + logical, save :: first + + data a0 /5.0d0/, p /1.0d0/ + data first /.true./ + + if (first) then + a0 = a0*drad + if (nparam .ge. 1) then + a0 = param(1) + end if + if (nparam .ge. 2) then + p = param(2) + end if + first = .false. + end if + + angle = mod(inc, TwoPi) + if ((angle .gt. 0.0d0) .and. (angle .lt. a0)) then + cold_inc_sq = angle**p*(a0**p - angle**p) + else + cold_inc_sq = 0.0d0 + end if + + return + end function cold_inc_sq + + real (kind=8) function warm_inc_sq (nparam, param, inc) +!-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- +! This routine returns the unnormalized inclination "probability" +! density for warm objects. +!-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- +! +! J-M. Petit Observatoire de Besancon +! Version 1 : April 2023, 7th +! +!-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- +! INPUT +! nparam: Number of parameters (I4) +! param : Parameters (n*R8) +! inc : Inclination [rad] (R8) +! +! OUPUT +! warm_inc_sq: Value of the probability (R8) +!-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- +!f2py intent(in) nparam +!f2py intent(in), depend(nparam) :: param +!f2py intent(in) inc + implicit none + + integer (kind=4), intent(in) :: nparam + real (kind=8), intent(in) :: param(*), inc + real (kind=8), parameter :: Pi = 3.141592653589793238d0, & + TwoPi = 2.0d0*Pi, drad = Pi/180.0d0 + real (kind=8), save :: angle, a0, p + logical, save :: first + + data a0 /12.0d0/, p /1.0d0/ + data first /.true./ + + if (first) then + a0 = a0*drad + if (nparam .ge. 1) then + a0 = param(1) + end if + if (nparam .ge. 2) then + p = param(2) + end if + first = .false. + end if + + angle = mod(inc, TwoPi) + if ((angle .gt. 0.0d0) .and. (angle .lt. a0)) then + warm_inc_sq = angle**p*(a0**p - angle**p) + else + warm_inc_sq = 0.0d0 + end if + + return + end function warm_inc_sq + + real (kind=8) function H_dist_cold_2(seed, nparam, hparam) +!-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- +! This routine draws randomly a number according to the cold belt H_r +! distribution, represented by an exponentially tapered exponential, +! with parameters I've fitted on the OSSOS cold belt data. +!-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- +! +! J-M. Petit Observatoire de Besancon +! Version 1 : October 2021 +! Version 2 : December 2021- updated parameter values. +! Version 3 : January 2022 - forcing slope at small sizes. +! +!-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- +! INPUT +! seed : Random number generator seed (I4) +! nparam: Number of parameters (I4) +! hparam: Parameters for asymptotic slope(s) (n*R8) +! hparam(1): start of asymptote +! hparam(2): contrast at start of asymptote +! hparam(3): slope of asymptote +! hparam(4): end of asymptote +! +! OUTPUT +! H_dist_cold : Random value of H (R8) +!-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- +!f2py intent(in,out) seed + implicit none + + integer (kind=4), parameter :: np = 16384 + integer (kind=4), intent(inout) :: seed + integer (kind=4), intent(in) :: nparam + real (kind=8), intent(in) :: hparam(4) + integer (kind=4) :: i + real (kind=8) :: params(4), random, h_min, h_max, n, c + real (kind=8), save :: proba(0:np), htab(0:np) + logical, save :: first + + data first /.true./ + data params /-2.6d0, 8.1d0, 0.666d0, 0.42d0/ + data h_min /4.6d0/ + + if (first) then + h_max = hparam(nparam) + htab(0) = h_min + proba(0) = 1.d-10 + n = Variably_tapered(hparam(1), params) + c = hparam(2) + do i = 1, np + htab(i) = h_min + dfloat(i)*(h_max-h_min)/dfloat(np) + if (htab(i) .lt. hparam(1)) then + proba(i) = Variably_tapered(htab(i), params) + else + proba(i) = n & + + n*c*(10.0d0**(hparam(3)*(htab(i)-hparam(1))) - 1.0d0) + end if + end do + do i = 0, np + proba(i) = proba(i)/proba(np) + end do + first = .false. + end if + + random = ran_3(seed) + H_dist_cold_2 = interp(proba, htab, random, np+1) + + return + end function H_dist_cold_2 + + subroutine H_diff_hot_5(nparam, hparam, np, hs, dist) +!-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- +! This routine computes the differential distribution of the hot belt H_r, +! represented by an exponentially tapered exponential, with parameters +! fitted on the OSSOS cold belt data, then scaled to hot. +! This version provides a continuous differential function, except for the divot. +! +! This version has modified parameters to obtain a cumulative distribution +! that looks like, and has the same normalisation as H_dist_hot_3. +! +!-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- +! +! J-M. Petit Observatoire de Besancon +! Version 1 : June 2024 - From H_dist_hot_4 +! +!-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- +! INPUT +! nparam: Number of parameters (I4) +! hparam: Parameters for asymptotic slope(s) (n*R8) +! hparam(1): start of asymptote +! hparam(2): contrast at start of asymptote +! hparam(3): slope of asymptote +! hparam(4): end of asymptote +! np : Number of points to return +! hs : Values of H at which we want the distribution +! +! OUTPUT +! hs(0) : H_min +! dist : Values of the cumulative distribution +!-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- +!f2py intent(in) nparam +!f2py intent(in), depend(nparam) :: hparam +!f2py intent(in) np +!f2py intent(in, out) hs +!f2py intent(out) dist + implicit none + + integer (kind=4), intent(in) :: nparam, np + real (kind=8), intent(in) :: hparam(*) + real (kind=8), intent(inout) :: hs(0:np) + real (kind=8), intent(out) :: dist(0:np) + integer (kind=4) :: i + real (kind=8), save :: params(4), h_min, h_max + real (kind=8), save :: a1, a2, a3, sl1, sl2, sl3, h1, h2, h3, scale + + data params /-2.6d0, 8.1d0, 0.666d0, 0.42d0/ + data h_min /-1.d0/ + data sl1 /0.13d0/, sl2 /0.7d0/ + data h1 /2.5d0/, h2 /5.8d0/ + data scale /2.2d0/ + + h_max = hparam(nparam) + h3 = hparam(1) + a2 = scale*Variably_tapered_diff(h2, params) + a3 = hparam(2)*scale*Variably_tapered_diff(h3, params) + sl3 = hparam(3) + a1 = a2*10.0d0**(sl2*(h1-h2)) +! The normalisation is done with the exponentially tapered exponential, +! as fitted on the OSSOS cold component, then scaled by 2. This +! determines the normalisation of the exponential between H = h2 +! and H = h1 +! N( from reference frame to ecliptic ; +! -1 ==> from ecliptic to reference frame +! v_in : input vector v_in(3) in the initial frame +! eps : Epsilon, inclination of ref frame [rad] +! om : Omega, longitude of node [rad] +! +! OUTPUT +! v_out : output vector in the final frame v_out(3) +! ierr : ierr = 0 :: normal exit +! ierr = 100 :: error ieqec neither 1 or -1 . +! +!*************************************************************** +!f2py intent(in) ieqec +!f2py intent(in) v_in +!f2py intent(in) eps +!f2py intent(in) om +!f2py intent(out) v_out +!f2py intent(out) ierr + implicit none + + type(t_v3d), intent(in) :: v_in + type(t_v3d), intent(out) :: v_out + integer, intent(in) :: ieqec + integer, intent(out) :: ierr + real (kind=8), intent(in) :: eps, om +! obliquity at J20000 in arcsec + real (kind=8), parameter :: Pi = 3.141592653589793238d0, & + secrad = Pi/180.d0/3600.d0 + real (kind=8) :: ww(3), coseps, sineps, cosom, sinom + + coseps = dcos(eps) + sineps = dsin(eps) + cosom = dcos(om) + sinom = dsin(om) + +! regular exit + ierr = 0 + +! to allow a call like :: call ref_ecl(ieqec,vv,vv,ierr) + ww(1) = v_in%x + ww(2) = v_in%y + ww(3) = v_in%z + + if (ieqec .eq. 1) then + v_out%x = cosom*ww(1) - sinom*(coseps*ww(2) - sineps*ww(3)) + v_out%y = sinom*ww(1) + cosom*(coseps*ww(2) - sineps*ww(3)) + v_out%z = sineps*ww(2) + coseps*ww(3) + else if (ieqec .eq. -1) then + v_out%x = cosom*ww(1) + sinom*ww(2) + v_out%y = coseps*(-sinom*ww(1) + cosom*ww(2)) + sineps*ww(3) + v_out%z = - sineps*(-sinom*ww(1) + cosom*ww(2)) + coseps*ww(3) + else +! anomalous exit ieqec not allowed + ierr = 100 + end if + + return + end subroutine ref_ecl + + subroutine ref_ecl_osc(ieqec, o_mi, o_mo, eps, om, ierr) +!-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- +! This routine converts osculating elements back and forth between +! given reference frame and ecliptic plane. +! Uses ref_ecl to do the work. +! Reference frame is given by Omega (om), the node longitude (rotation +! around Z axis of ecliptic) and Epsilon (eps), the inclination of the +! reference frame (rotation around node axis of ecliptic). +! +! ANGLES ARE GIVEN IN RADIAN !!!! +! +!-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- +!f2py intent(in) ieqec +!f2py intent(in) o_mi +!f2py intent(in) eps +!f2py intent(in) om +!f2py intent(out) o_mo +!f2py intent(out) ierr + implicit none + + type(t_orb_m), intent(in) :: o_mi + type(t_orb_m), intent(out) :: o_mo + integer, intent(in) :: ieqec + integer, intent(out) :: ierr + real (kind=8), intent(in) :: eps, om + real (kind=8), parameter :: Pi = 3.141592653589793238d0, TwoPi = 2.d0*Pi, & + drad = Pi/180.d0, mu = TwoPi**2 + type(t_orb_m) :: o_md + type(t_v3d) :: posi, poso, veli, velo + + o_md%a = o_mi%a + o_md%e = o_mi%e + o_md%inc = o_mi%inc + o_md%node = o_mi%node + o_md%peri = o_mi%peri + o_md%m = o_mi%m + call coord_cart (mu, o_md, posi, veli) + call ref_ecl (ieqec, posi, poso, eps, om, ierr) + call ref_ecl (ieqec, veli, velo, eps, om, ierr) + call osc_el (mu, poso, velo, o_mo) + call ztopi (o_mo%inc) + if (o_mo%inc .gt. Pi) o_mo%inc = o_mo%inc - Pi + call ztopi (o_mo%node) + call ztopi (o_mo%peri) + call ztopi (o_mo%m) + + return + end subroutine ref_ecl_osc + + subroutine forced_plane_damp(a, inc, ifd, Omfd) +!-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-* +! This routine returns the inclination iforced and node Omforced of the +! forced plane as a function of semimajor-axis in the region of the +! main classical Kuiper belt, according to second order theory with all +! 8 planets. +! +! This version damps the forced plane to the invariable plane according +! to inclination. +! +! The fit is valid only between 35 au and 47.74 au (2:1 MMR). +! Returns 0 for a < 35 au, and the invariable plane for a > 47.74 au +! +! Returns degrees. +! +!-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-* +! +! J-M. Petit Observatoire de Besancon +! Version 1 : November 2021 +! +!-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-* +! INPUT +! a : semimajor-axis [au] (R8) +! inc : Inclination [deg] (R8) +! +! OUTPUT +! ifd : Inclination of forced plane [deg] (R8) +! Omfd : Node of forced plane [deg] (R8) +! +!-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-* +! +! Set of F2PY directives to create a Python module +! +!f2py intent(in) a +!f2py intent(in) inc +!f2py intent(out) ifd +!f2py intent(out) Omfd + implicit none + +! Calling arguments + real (kind=8), intent(in) :: a, inc + real (kind=8), intent(out) :: ifd, Omfd + +! Internal variables + real (kind=8) :: ci0(4), ci1(3), ci2(3), ci3(3), co0(3), co1(3), co2(5), & + co3(5), alpha, beta, gamma, delta, om_lim_low, om_lim_high, damp, y + real (kind=8), parameter :: epsilon = 5713.86d0/3600.d0, & + omega = 387390.8d0/3600.d0, ce_damp = 16.d0, wi_damp = 6.d0 + + data & + ci0 /-0.115657583d0,34.8097343d0,7.79198557d-02,-1.06252408d0/, & + ci1 /-0.392018467d0, 40.5282974d0, 0.137285933/, & + ci2 /-0.391531527d0, 40.6837921d0, 0.190758109/, & + ci3 /0.391087621d0, 40.2941360d0, 0.146276891/, & + co0 /971.35006612732775d0, -50.061061665365997d0, & + 0.74715152013989472d0/, & + co1 /4694.6567730091410d0, -249.62252879980477d0, & + 3.4214148236506192d0/, & + co2 /11.7058372d0, -294.139587d0, 30.8862801d0, & + -1049.41772d0, 15.2022839d0/, & + co3 /30.4611244d0, -1203.70593d0, 2.47961617d0, & + -16.6130295d0, 302.805359d0/, & + alpha /1.d0/, & + om_lim_low /200.d0/, & + om_lim_high /20.d0/ + + common /om_lim_com/ om_lim_low, om_lim_high + + damp(y) = 1.d0*(1.d0-tanh((y-ce_damp)/wi_damp))/2.d0 + + if (a .lt. 35.d0) then + ifd = 0.d0 + Omfd = 0.d0 + else if (a .lt. 37.d0) then + ifd = (ci0(1)/(a - ci0(2)) + ci0(3)*a + ci0(4)) + Omfd = co0(1) + co0(2)*a + co0(3)*a*a + else if (a .lt. 39.0d0) then + ifd = 10.d0**(ci1(1)/(a - ci1(2)) + ci1(3)) + Omfd = co1(1) + co1(2)*a + co1(3)*a*a + else if (a .lt. 40.0d0) then + ifd = 2.477d0 + Omfd = (om_lim_low-163.187d0)/(40.0d0-39.d0)*(a-40.0d0) & + + om_lim_low + else if (a .lt. 41.9d0) then + ifd = 2.477d0 + Omfd = (61.791d0-om_lim_high)/(41.9d0-40.0d0)*(a-40.0d0) & + + om_lim_high + else if (a .lt. 47.74d0) then + ifd = 10.d0**(ci3(1)/(a - ci3(2)) + ci3(3)) + beta = -((co3(1) + co3(3))*a + co3(2) + co3(4)) + gamma = (co3(1)*a + co3(2))*(co3(3)*a + co3(4)) - co3(5) + delta = max(0., beta**2 - 4.d0*alpha*gamma) + Omfd = (-beta - sqrt(delta))/(2.d0*alpha) + else + ifd = epsilon + Omfd = omega + end if + ifd = min(ifd, 40.d0) + ifd = epsilon + (ifd-epsilon)*damp(inc) + Omfd = omega + (Omfd-omega)*damp(inc) + + return + + end subroutine forced_plane_damp + + subroutine ztopi (var) +!-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- +! This subroutine resets var to be between 0 and 2Pi. +!-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*- +!f2py intent(in,out) var + implicit none + + real (kind=8), intent(inout) :: var + real (kind=8), parameter :: Pi = 3.141592653589793238d0, TwoPi = 2.d0*Pi + +1000 continue + if (var .gt. TwoPi) then + var = var - TwoPi + goto 1000 + end if +1100 continue + if (var .lt. 0.d0) then + var = var + TwoPi + goto 1100 + end if + + return + end subroutine ztopi + +end module rot diff --git a/eupl1.1.-en_0_0.pdf b/eupl1.1.-en_0_0.pdf new file mode 100644 index 0000000..d69ea0d Binary files /dev/null and b/eupl1.1.-en_0_0.pdf differ diff --git a/eupl1.1.-licence-en_0.pdf b/eupl1.1.-licence-en_0.pdf new file mode 100644 index 0000000..80b5007 Binary files /dev/null and b/eupl1.1.-licence-en_0.pdf differ