! ORBIT_ELEMENTS
! contains:
!
! OUT OF MOD:
!      appmag     correction to absolute magnitude to get apparent


MODULE orbit_elements
USE output_control
USE fund_const
USE dyn_param
IMPLICIT NONE
PRIVATE

! public routines
PUBLIC coo_cha, equal_orbels, write_elems_old, read_elems_old, ecc_peri
PUBLIC convertunc, propagunc 
PUBLIC cootyp, attelements, eldiff, tpcar_mtpcar, opik_tpcar 
PUBLIC centerbody_mass, att_prelim2 ! ??? used by bizarre
PUBLIC find_elems, wrikep, wromlr,rd_orb !old I/O
PUBLIC write_elems, read_elems           ! new I/O 
PUBLIC wro1lh2_ngr, wromlh2, wro1lh2, wro1lr     ! new I/O 
PUBLIC wro1lr_matlab
PUBLIC undefined_orb_uncert
PUBLIC rdelem

DOUBLE PRECISION :: eccbou ! boundary in e between use of Kepler equation
                           ! and use of f-g series

! Generic orbital elements set
TYPE orbit_elem
CHARACTER(LEN=3)       :: coo     ! elements type; known types are:
!    'CAR' : cartesian positions and velocities 
!    'EQU' : equinoctal elements
!    'KEP' : classical keplerian elements, singular for
!             0 eccentricity, 0 and 180 degrees inclination;
!    'COM' : cometary elements, valid for ecc.ge.1 
!    'ATT' : attributable plus r rdot
!    'OPI' : Opik type elements not yet implemented

!INTEGER :: ndim  ! actual dimension of elements vector

DOUBLE PRECISION, DIMENSION(6) :: coord ! elements vector 
              ! DANGER: some of the coordinates may actually be angles...   
DOUBLE PRECISION :: t  ! epoch time, MJD, TDT

LOGICAL :: mag_set ! true if some absolute magnitude is available
DOUBLE PRECISION :: h_mag, g_mag ! Absolute magnitude H, opposition effect G

INTEGER :: center ! 0 for Sun, code for planet as in JPL ephemerides

INTEGER :: obscode ! observatory code, for ATT type only
 
END TYPE orbit_elem

! default value when undefined
DOUBLE PRECISION, DIMENSION(6), PARAMETER :: zero_6d_vect = &
&    (/ 0.d0, 0.d0, 0.d0, 0.d0, 0.d0, 0.d0 /)

! undefined orbital elements
TYPE(orbit_elem), PARAMETER :: undefined_orbit_elem = ORBIT_ELEM( &
&  'UNK',           & ! no coordynate defined
& zero_6d_vect,     & ! null elements
&   1.d99,          & ! null epoch 
&   .false.,        & ! default no magnitudes
&   -9.99d0,        & ! null magnitude
&    0.15d0,        & ! default opposition effect
&    0,             & ! default is heliocentric
&    500            & ! default is not topocentric correction 
! &  ,gms            & ! if needed
&     )  ! 

TYPE orb_uncert
! normal and covariance matrix (warning: covariance not full if not 
! all solved)
DOUBLE PRECISION, DIMENSION(ndimx,ndimx):: c,g
! norm of residuals, of last correction, RMS of magnitudes
LOGICAL :: succ ! logical for success (convergent differential corrections)
INTEGER :: ndim  ! actual dimension of solve for vector
END TYPE orb_uncert

! default value when undefined
DOUBLE PRECISION, DIMENSION(6,6), PARAMETER :: zero_6x6_matrix = &
& RESHAPE((/ 0.d0, 0.d0, 0.d0, 0.d0, 0.d0, 0.d0, &
& 0.d0, 0.d0, 0.d0, 0.d0, 0.d0, 0.d0, 0.d0, 0.d0, 0.d0, 0.d0, 0.d0, &
& 0.d0, 0.d0, 0.d0, 0.d0, 0.d0, 0.d0, 0.d0, &
& 0.d0, 0.d0, 0.d0, 0.d0, 0.d0, 0.d0, 0.d0, 0.d0, 0.d0, 0.d0, 0.d0, 0.d0 /),(/ 6, 6 /))

TYPE orb_with_uncert
TYPE(orbit_elem) :: orbit ! nominal orbit, solution of least squares
TYPE(orb_uncert) :: uncertainty ! covariance matrix etc.
END TYPE orb_with_uncert

! public entities
PUBLIC orbit_elem,undefined_orbit_elem,orb_uncert
PUBLIC eccbou

CONTAINS

SUBROUTINE undefined_orb_uncert(nd,unc)
 INTEGER, INTENT(IN) :: nd ! dimension of matrices
 TYPE(orb_uncert), INTENT(OUT) :: unc
 unc%ndim=nd
 unc%succ=.false. 
 unc%g(1:nd,1:nd)=0.d0
 unc%c(1:nd,1:nd)=0.d0
END SUBROUTINE undefined_orb_uncert

! ===================================================================   
! COO_CHA                                                         
! ===================================================================   
!   general purpose elements change 
!    
!  fail _flag = 0 conversion OK
!             = 1 zero eccentricity (conv. to KEP/COM fails)
!             = 2 zero inclination               "
!             = 3 zero ecc and zero inclination   "
!             = 5 eccentricity 1 and more (conv to EQU/KEP fails) 
!             = 9 180 deg inclination (conv to KEP/COM/EQU fails)
!             = 10 zero angular momentum (conv to EQU/KEP/COM fails)
!             = 20 inconsistent calling sequence 
!  if fail_flag =1,2,3 the jacobian is undefined for kep_equ, replaced
!                by an arbitrary matrix (to avoid undefined on exit)
!  if fail_flag =5 conversion to KEP, EQU is refused, something else
!                  is provided in output (should be COM)
!  if fail_flag > 8 there is no conversion provided....
!               we should set up the messages in a suitable way,
!               providing for a smooth failure
! =================================================================== 
SUBROUTINE coo_cha(x,cooy,y,fail_flag,del,obscode)
! MODIFICATION 2009/debris: not allowed to change center, by passing argument iplanet
TYPE(orbit_elem), INTENT(IN) :: x
CHARACTER*3, INTENT(IN) :: cooy
TYPE(orbit_elem), INTENT(INOUT) :: y ! needed to allow for x=y
INTEGER, INTENT(OUT) :: fail_flag ! see above
! partial derivatives \partial y/\partial x
DOUBLE PRECISION, INTENT(OUT), OPTIONAL, DIMENSION(6,6) :: del
INTEGER, INTENT(IN), OPTIONAL :: obscode ! output observer
! END INTERFACE
TYPE(orbit_elem) z
LOGICAL kepgr_in,kepgr_ou,cargr_in,cargr_ou,opposite, done
CHARACTER*3 coox
! target_center is the output gravitational center
INTEGER obscod,target_center
DOUBLE PRECISION, DIMENSION(6,6) :: del1,del2 

done=.false.
fail_flag=0
coox=x%coo
IF(rhs.eq.1.or.rhs.eq.3)THEN
   target_center=0
ELSEIF(rhs.eq.2)THEN
   target_center=3
ELSE
   WRITE(*,*)' coo_cha: rhs not initialized, rhs=', rhs
   STOP
ENDIF

IF(coox.eq.cooy)THEN
   y=x
   IF(PRESENT(del)) CALL eye(6,del)
   done=.true.
   RETURN
ENDIF
IF(.not.PRESENT(obscode).and.(cooy.eq.'ATT'))THEN
   obscod=500
ELSEIF(PRESENT(obscode))THEN
   obscod=obscode
ENDIF
kepgr_in=(coox.eq.'KEP'.or.coox.eq.'EQU'.or.coox.eq.'COM'.or.coox.eq.'COT')
kepgr_ou=(cooy.eq.'KEP'.or.cooy.eq.'EQU'.or.cooy.eq.'COM'.or.cooy.eq.'COT')
cargr_in=(coox.eq.'CAR'.or.coox.eq.'ATT')
cargr_ou=(cooy.eq.'CAR'.or.cooy.eq.'ATT')
!opposite=((kepgr_in.and.cargr_ou).or.(cargr_in.and.kepgr_ou))

IF(cargr_in)THEN
   IF(coox.eq.'ATT')THEN
!      IF(planet.eq.0)THEN
! attributable case must be converted to cartesian heliocentric
      IF(PRESENT(del))THEN
         z=car_att(x,target_center,del1) ! convert to CAR
      ELSE
         z=car_att(x,target_center)      ! convert to CAR
      ENDIF
   ELSEIF(coox.eq.'CAR')THEN
      z=x                                 ! already cartesian
      IF(PRESENT(del))THEN
         CALl eye(6,del1)     ! initialize by identity
      ENDIF
   !   ELSEIF(coox.eq.'OPI')THEN
   !      WRITE(*,*)' coo_cha: this conversion not ready ', coox, ' ',cooy
   !      STOP
   ENDIF
   IF(cooy.eq.'CAR')THEN
      y=z ! done
      IF(PRESENT(del))del=del1
   ELSEIF(cooy.eq.'ATT')THEN
      IF(PRESENT(del))THEN
         y=att_car(z,obscod,target_center,del2)
         del=MATMUL(del2,del1)
      ELSE
         y=att_car(z,obscod,target_center)
      ENDIF
   ELSEIF(kepgr_ou)THEN
! convert to elements  
      IF(PRESENT(del))THEN
         y=to_elems(z,cooy,fail_flag,del2)
         del=MATMUL(del2,del1)
      ELSE
         y=to_elems(z,cooy,fail_flag)
      ENDIF
! handle fail_flag
      IF(fail_flag.gt.0.and.fail_flag.le.3)THEN
! zero eccentricity/inclination
! EQU are differentiable
         IF(cooy.eq.'EQU')fail_flag=0 
! KEP, COM are defined but not diff.
         IF(.not.PRESENT(del))fail_flag=0 
      ELSEIF(fail_flag.eq.5)THEN 
! force to COT, leaving fail_flag=5 as warning
         IF(PRESENT(del))THEN
 ! should not be necessary
            IF(y%coo.ne.'COT')THEN
               WRITE(ierrou,*)' coo_cha: fail_flag, outcoo ', fail_flag, y%coo
               numerr=numerr+1
               y=to_elems(z,'COT',fail_flag,del2)
            ENDIF
            del=MATMUL(del2,del1)
         ELSE
            IF(y%coo.ne.'COT')THEN ! should not be necessary
               WRITE(ierrou,*)' coo_cha: fail_flag, outcoo ', fail_flag, y%coo
               numerr=numerr+1
               y=to_elems(z,'COT',fail_flag)
            ENDIF
         ENDIF
      ELSE
! fail_flag=9,10,20: nothing to do
      ENDIF
      done=.true.
      RETURN
   ELSE
      PRINT *, ' coo_cha: this should not happen 1, coox=',coox,' cooy=', cooy
      WRITE(*,*) x
      STOP
   ENDIF
ENDIF
IF(kepgr_in.and.cargr_ou)THEN
   IF(PRESENT(del))THEN
      z=car_elems(x,fail_flag,del1)
   ELSE
      z=car_elems(x,fail_flag)
   ENDIF
   IF(cooy.eq.'ATT')THEN
      IF(PRESENT(del))THEN
         y=att_car(z,obscod,target_center,del2)
         del=MATMUL(del2,del1)
      ELSE
         y=att_car(z,obscod,target_center)
      ENDIF
   ELSEIF(cooy.eq.'CAR')THEN
      y=z ! done
      IF(PRESENT(del))del=del1
   ENDIF
! handle fail_flag
   IF(fail_flag.gt.0.and.fail_flag.le.3)THEN
      IF(coox.eq.'EQU')fail_flag=0 ! EQU are differentiable
      IF(.not.PRESENT(del))fail_flag=0 ! KEP, COM are defined but not diff.
   ELSE
! fail_flag=5 cannot happen, 9,10,20 cannot be fixed
   ENDIF
   done=.true.
   RETURN
ENDIF
IF(kepgr_in.and.kepgr_ou)THEN ! twelfe cases
! starting from KEP
   IF(coox.eq.'KEP')THEN
      IF(cooy.eq.'EQU')THEN
         IF(PRESENT(del))THEN
            y=equ_kep(x,del)
         ELSE
            y=equ_kep(x)
         ENDIF
      ELSEIF(cooy.eq.'COM')THEN
         IF(PRESENT(del))THEN
            y=com_kep(x,fail_flag,del)
         ELSE
            y=com_kep(x,fail_flag)
         ENDIF
      ELSEIF(cooy.eq.'COT')THEN
         IF(PRESENT(del))THEN
            y=cot_kep(x,fail_flag,del)
         ELSE
            y=cot_kep(x,fail_flag)
         ENDIF
      ENDIF
      done=.true.
      RETURN
   ENDIF
! starting from EQU
   IF(coox.eq.'EQU')THEN
      IF(cooy.eq.'KEP')THEN
         IF(PRESENT(del))THEN
            y=kep_equ(x,fail_flag,del) 
            IF(fail_flag.gt.0) CALl eye(6,del)! del non exi, set to identity
         ELSE
            y=kep_equ(x,fail_flag)
         ENDIF
      ELSEIF(cooy.eq.'COM')THEN
         IF(PRESENT(del))THEN
            z=kep_equ(x,fail_flag,del1)
            IF(fail_flag.eq.0)THEN
               y=com_kep(z,fail_flag,del2)
               del=MATMUL(del2,del1)
            ELSE
               y=com_kep(z,fail_flag)
               del=del1 ! del2 not existent, set to identity 
            ENDIF
         ELSE
            z=kep_equ(x,fail_flag) ! always existent, maybe not C1
            y=com_kep(z,fail_flag) 
         ENDIF
      ELSEIF(cooy.eq.'COT')THEN
         IF(PRESENT(del))THEN
            z=kep_equ(x,fail_flag,del1)
            IF(fail_flag.eq.0)THEN
               y=cot_kep(z,fail_flag,del2)
               del=MATMUL(del2,del1)
            ELSE
               y=cot_kep(z,fail_flag)
               del=del1 ! del2 not existent, set to identity 
            ENDIF
         ELSE
            z=kep_equ(x,fail_flag) ! always existent, maybe not C1
            y=cot_kep(z,fail_flag) 
         ENDIF
      ENDIF
      done=.true.
      RETURN
   ENDIF
! starting from COM
   IF(coox.eq.'COM')THEN
      IF(cooy.eq.'KEP')THEN
         IF(PRESENT(del))THEN
            y=kep_com(x,fail_flag,del)
         ELSE
            y=kep_com(x,fail_flag)
         ENDIF
      ELSEIF(cooy.eq.'EQU')THEN
         IF(PRESENT(del))THEN
            z=kep_com(x,fail_flag,del1)
            IF(fail_flag.lt.4)THEN
               y=equ_kep(z,del2)
               del=MATMUL(del2,del1)
            ELSE
               y=z
               del=del1 ! del 2 non existent, set to identity
            ENDIF
         ELSE
            z=kep_com(x,fail_flag)
            IF(fail_flag.lt.4)THEN
               y=equ_kep(z)
            ELSE
               y=z
            ENDIF
         ENDIF
      ELSEIF(cooy.eq.'COT')THEN
         IF(PRESENT(del))THEN
            y=cot_com(x,fail_flag,del)
         ELSE
            y=cot_com(x,fail_flag)
         ENDIF
      ENDIF
      done=.true.
      RETURN
   ENDIF
! starting from COT
   IF(coox.eq.'COT')THEN
      IF(cooy.eq.'KEP')THEN
         IF(PRESENT(del))THEN
            y=kep_cot(x,fail_flag,del)
         ELSE
            y=kep_cot(x,fail_flag)
         ENDIF
      ELSEIF(cooy.eq.'COM')THEN
         IF(PRESENT(del))THEN
            y=com_cot(x,fail_flag,del) 
            IF(fail_flag.gt.0) CALl eye(6,del)! del non exi, set to identity
         ELSE
            y=com_cot(x,fail_flag)
         ENDIF
      ELSEIF(cooy.eq.'EQU')THEN
         IF(PRESENT(del))THEN
            z=kep_cot(x,fail_flag,del1)
            IF(fail_flag.lt.4)THEN
               y=equ_kep(z,del2)
               del=MATMUL(del2,del1)
            ELSE
               y=z
               del=del1 ! del 2 non existent, set to identity
            ENDIF
         ELSE
            z=kep_cot(x,fail_flag)
            IF(fail_flag.lt.4)THEN
               y=equ_kep(z)
            ELSE
               y=z
            ENDIF
         ENDIF
      ENDIF
      done=.true.
      RETURN
   ENDIF
   IF(.not.done)THEN
      PRINT *, ' coo_cha: this should not happen 2, coox=',coox,' cooy=', cooy
      WRITE(*,*) x
      STOP
   ENDIF
ENDIF
END SUBROUTINE coo_cha
! ====================================================================
! UTILITIES
! ====================================================================
LOGICAL FUNCTION equal_orbels(el1,el2)
  TYPE(orbit_elem), INTENT(IN) :: el1,el2
  INTEGER i
  equal_orbels=.false.
  IF(el1%coo.ne.el2%coo.or.el1%center.ne.el2%center)RETURN
  DO i=1,6
     IF(el1%coord(i).ne.el2%coord(i))RETURN
  ENDDO
  equal_orbels=.true.
!magnitude does not matter!
END FUNCTION equal_orbels

SUBROUTINE eldiff(el1,el0,del)
  TYPE(orbit_elem), INTENT(IN) :: el1,el0
  DOUBLE PRECISION, DIMENSION(6) :: del
  DOUBLE PRECISION pridif
  IF(el1%coo.ne.el0%coo)THEN
     WRITE(ierrou,*)'eldiff: different coordinates ',el1%coo, el0%coo
     numerr=numerr+1
     del=0.d0
     RETURN
  ENDIF
  IF(el0%coo.eq.'CAR')THEN
     del=el1%coord-el0%coord
  ELSEIF(el0%coo.eq.'EQU')THEN
     del(1:5)=el1%coord(1:5)-el0%coord(1:5)
     del(6)=pridif(el1%coord(6),el0%coord(6))
  ELSEIF(el0%coo.eq.'KEP')THEN
     del(1:3)=el1%coord(1:3)-el0%coord(1:3)
     del(4)=pridif(el1%coord(4),el0%coord(4))
     del(5)=pridif(el1%coord(5),el0%coord(5))
     del(6)=pridif(el1%coord(6),el0%coord(6))
  ELSEIF(el0%coo.eq.'COM')THEN
     del(1:3)=el1%coord(1:3)-el0%coord(1:3)
     del(4)=pridif(el1%coord(4),el0%coord(4))
     del(5)=pridif(el1%coord(5),el0%coord(5))
     del(6)=el1%coord(6)-el0%coord(6)
  ELSEIF(el0%coo.eq.'ATT')THEN
     del(2:6)=el1%coord(2:6)-el0%coord(2:6)
     del(1)=pridif(el1%coord(1),el0%coord(1))
  ELSE
     WRITE(*,*)'eldiff: unknown coordinates ',el1%coo, el0%coo
     STOP
  ENDIF
END SUBROUTINE eldiff

! ====================================================================
! ECC_PERI
! computes eccentricty, pericenter, apocenter, mean motion
! if energy>0 then enne=0, qg=1.d+50
SUBROUTINE ecc_peri(el,ecc,q,qg,enne)
TYPE(orbit_elem), INTENT(IN) :: el
DOUBLE PRECISION, INTENT(OUT) :: ecc,q,qg,enne
! ===================================================================
DOUBLE PRECISION gm,eps,a
DOUBLE precision x(3),y(3),ang(3),vlenz(3),gei,gei2,r0,vel2,alpha,prscal,vsize
TYPE(orbit_elem) el1
INTEGER fail_flag
!
eps=100*epsilon(1.d0)
IF(rhs.eq.1)THEN
   gm=centerbody_mass(0) ! we want always eccentricity heliocentric
ELSEIF(rhs.eq.2)THEN
   gm=centerbody_mass(3) ! we want always eccentricity geocentric
ELSE
! WARNING: automatic choice of method/dynamical model not supported in ORBIT9
   WRITE(*,*)'ecc_peri: ORBIT9 not supported here '
   STOP
ENDIF 
IF(el%coo.eq.'EQU')THEN
   ecc=sqrt(el%coord(2)**2+el%coord(3)**2)
   q=el%coord(1)*(1.d0-ecc)
   qg=el%coord(1)*(1.d0+ecc)
   enne=sqrt(gm/el%coord(1)**3)
ELSEIF(el%coo.eq.'KEP')THEN
   ecc=el%coord(2)
   q=el%coord(1)*(1.d0-ecc)
   qg=el%coord(1)*(1.d0+ecc)
   enne=sqrt(gm/el%coord(1)**3)
ELSEIF(el%coo.eq.'COM'.or.el%coo.eq.'COT')THEN
   ecc=el%coord(2)
   q=el%coord(1)
   IF(ecc.gt.1.d0-eps)THEN
      qg=1.d+50
      enne=0.d0
   ELSE
      a=el%coord(1)/(1.d0-ecc)
      qg=a*(1.d0+ecc)
      enne=sqrt(gm/a**3)
   ENDIF
ELSE
   IF(el%coo.eq.'ATT')THEN
! convert to cartesian if it is ATT
      CALL coo_cha(el,'CAR',el1,fail_flag)
   ELSEIF(el%coo.eq.'CAR')THEN
      el1=el
   ELSE
      WRITE(*,*)' ecc_peri: unknown coord type, coo= ',el%coo
      WRITE(*,*) el
      STOP
   ENDIF
! compute eccentricity
   x=el1%coord(1:3)
   y=el1%coord(4:6)
!  radius and velocity squared                                          
   vel2=prscal(y,y)          ! velocity
   r0=vsize(x)               ! distance from center
   call prvec(x,y,ang)       !  angular momentum 
! non singular first element
   gei2=prscal(ang,ang)
   gei=sqrt(gei2)
   IF(gei.eq.0.d0)THEN
      WRITE(*,*) ' ecc_peri: zero angular momentum ',el
   ENDIF
   call prvec(y,ang,vlenz) ! Lenz vector
   vlenz=vlenz*(1.d0/gm)-x(1:3)*(1.d0/r0)
   ecc=vsize(vlenz)
   q=gei2/(gm*(1.d0+ecc))
   alpha=vel2-2*gm/r0 ! 2* energy 
   IF(ecc.gt.1.d0-eps)THEN
      qg=1.d+50
      enne=0.d0
   ELSE
      qg=gei2/(gm*(1.d0-ecc))
      a=-gm/alpha
      enne=sqrt(gm/a**3)
   ENDIF
ENDIF
END SUBROUTINE ecc_peri

DOUBLE PRECISION FUNCTION centerbody_mass(iplanet)
USE planet_masses
INTEGER,INTENT(IN) :: iplanet
IF(iplanet.eq.0)THEN
! WARNING: does not work in ORBIT9
   IF(rhs.gt.2)THEN
      WRITE(*,*)' centerbody_mass: not ready for rhs=', rhs
      STOP
   ENDIF
   centerbody_mass=gms
ELSEIF(iplanet.eq.3)THEN
   centerbody_mass=gmearth
ELSEIF(iplanet.le.npla+iatrue.and.iplanet.ne.3)THEN
   centerbody_mass=gm(iplanet)
ELSE
   WRITE(*,*)' centerbody_mass: center not available ', iplanet
   STOP
ENDIF 
END FUNCTION centerbody_mass
!
TYPE(orbit_elem) FUNCTION coord_translate(el,iplanet)
TYPE(orbit_elem), INTENT(IN) :: el
INTEGER, INTENT(IN) :: iplanet ! planet; JPL ephem codes
DOUBLE PRECISION :: xpla(6)
IF(el%coo.ne.'CAR')THEN
   WRITE(*,*)' coord_translate: wrong input coord ',el
   STOP
ENDIF
coord_translate=el
! check planet, observers are on Earth! 
IF(iplanet.eq.3.and.el%center.eq.0)THEN
   CALL earcar(el%t,xpla,1)
   coord_translate%coord=el%coord-xpla
   coord_translate%center=3
ELSEIF(iplanet.eq.0.and.el%center.eq.3)THEN
   CALL earcar(el%t,xpla,1)
   coord_translate%coord=el%coord+xpla
   coord_translate%center=0
ELSE
   WRITE(*,*)' coord_translate: not yet available for ', iplanet, el%center
   STOP
ENDIF
END FUNCTION coord_translate
! ====================================================================
! END UTILITIES
! ====================================================================
! CARTESIAN TO/FROM ELEMENTS
! ====================================================================
! TO_ELEMS
! ====================================================================
TYPE(orbit_elem) FUNCTION to_elems(el,cooy,fail_flag,del)
USE ever_pitkin
TYPE(orbit_elem),INTENT(IN):: el ! must be cartesian
CHARACTER*3, INTENT(IN) :: cooy ! can be 'KEP', 'EQU', 'COM', 'COT'
INTEGER, INTENT(OUT) :: fail_flag ! error flag
! optional jacobian matrix. in output: derivatives d(to_elems)/d(el)
DOUBLE PRECISION, DIMENSION(6,6),INTENT(OUT),OPTIONAL :: del
DOUBLE PRECISION gm
DOUBLE PRECISION, DIMENSION(3) :: x0, y0
DOUBLE PRECISION, DIMENSION(3) :: x, y, ang ! position, velocity, ang. momentum
DOUBLE PRECISION, DIMENSION(3) :: vlenz, fb, gb, ww ! Lenz vector, broucke f and g, pos x vlenz
DOUBLE PRECISION :: gei2, gei, ecc2, ecc, varpi, div, tgi2, omnod
DOUBLE PRECISION :: prscal, vsize, princ, eps
DOUBLE PRECISION :: r0, vel2, alpha, sig0, peri, psi, dt
DOUBLE PRECISION :: enne, rad, beta, chk, ch, ck, fe, x2, y2, cosf, sinf
DOUBLE PRECISION :: eq(6), cosv,sinv,dd
TYPE(orbit_elem) :: elback ! for inversion
DOUBLE PRECISION, DIMENSION(6,6) :: del1
INTEGER :: fail_flag1, ising
DOUBLE PRECISION :: det
! paranoia check of the calling arguments
IF(el%coo.ne.'CAR')THEN
   WRITE(*,*)' to_elems: wrong coordinates in input',x
   fail_flag=20
   RETURN
ENDIF
IF(cooy.ne.'COM'.and.cooy.ne.'KEP'.and.cooy.ne.'EQU'.and.cooy.ne.'COT')THEN
   WRITE(*,*)' to_elems: wrong output coordinates', cooy
   fail_flag=20
   RETURN
ENDIF
! initializations
fail_flag=0
gm=centerbody_mass(el%center) ! mass of central body
eps=100*epsilon(1.d0)         ! small number for rounding off
to_elems=el                   ! copy defaults
! =====================================================================
! portion of the algorithm used for all output coordinates
! =====================================================================
x=el%coord(1:3)               ! cartesian position
y=el%coord(4:6)               ! cartesian velocity
vel2=prscal(y,y)              ! velocity squared
r0=vsize(x)                   ! distance from center
CALL prvec(x,y,ang)           ! angular momentum vector 
gei2=prscal(ang,ang)
gei=sqrt(gei2)                ! angular momentum scalar
IF(gei.eq.0.d0)THEN 
   WRITE(*,*) ' to_elems: zero angular momentum ',el 
   fail_flag=10
   RETURN
ENDIF
CALL prvec(y,ang,vlenz)       ! Lenz vector 
vlenz=vlenz*(1.d0/gm)-x(1:3)*(1.d0/r0) 
ang=ang*(1.d0)/gei            ! orbit normal unit vector
!   zero divide occurs for inclination of 180 degrees                   
div=1.d0+ang(3) 
IF(div.eq.0.d0)THEN 
   WRITE(*,*) ' to_elems: 180 deg. inclination ',el
   fail_flag=9
   RETURN
ENDIF   
!  unit vectors of the equinoctal reference system (Broucke and         
!  Cefola 1972, CM 5, 303--310) are fb, gb, ang                           
fb(1)=1.d0-ang(1)**2/div 
fb(2)=-ang(1)*ang(2)/div 
fb(3)=-ang(1) 
CALL prvec(ang,fb,gb) 
!  elements related to eccentricity and inclination                     
eq(2)=prscal(vlenz,gb) 
eq(3)=prscal(vlenz,fb) 
eq(4)=ang(1)/div 
eq(5)=-ang(2)/div 
ecc2=eq(2)**2+eq(3)**2 
ecc=sqrt(ecc2)                ! eccentricity
IF(ecc.lt.eps)THEN
   varpi=0.d0
   fail_flag=1
ELSE
   varpi=atan2(eq(2),eq(3))   ! longitude of pericenter
ENDIF                                          
tgi2=sqrt(eq(4)**2+eq(5)**2)  ! tanget of inclination/2
if(tgi2.lt.eps)then 
   omnod=0.d0 
   fail_flag=fail_flag+2
else 
   omnod=princ(atan2(eq(4),eq(5))) ! Omega (longitude of node) 
endif
alpha=vel2-2*gm/r0            ! 2* energy
IF(alpha.gt.-eps)THEN
  IF(cooy.eq.'EQU'.or.cooy.eq.'KEP') fail_flag=5 ! hyperbolic orbit; is a failure
! if requested coordinates in output are KEP, EQU
ENDIF
sig0=prscal(x,y)              ! scalar product as used in ever_pitkin
! =====================================================================
! part depending upon output coordinates
! =====================================================================
IF(cooy.eq.'COM'.or.cooy.eq.'COT'.or.fail_flag.eq.5)THEN
!         coord(1)=q=a(1-e)
!         coord(2)=e
!         coord(3)=I
!         coord(4)=Omega 
!         coord(5)=omega 
! if cooy=COM
!         coord(6)=t_0 (time at pericenter)
! elseif cooy=COT
!         coord(6)=v (true anomaly)
   eq(1)=gei2/(gm*(1.d0+ecc))
   eq(2)=ecc
   eq(5)=princ(varpi-omnod)
   eq(3)=2*atan(tgi2)
   eq(4)=omnod
   peri=eq(1)
   IF(cooy.eq.'COM')THEN
      IF(ecc.gt.eps)THEN
! WARNING: use of solve_peri may fail, add error flag which affects
!          the output fail_flag
         CALL solve_peri(r0,sig0,peri,gm,alpha,psi,dt)
         eq(6)=el%t+dt
         fail_flag=0
      ELSE
         fail_flag=1  ! numerically circular orbit
         eq(6)=el%t   ! always at perihelion
      ENDIF
      to_elems%coo='COM' 
   ELSEIF(cooy.eq.'COT'.or.fail_flag.eq.5)THEN
      dd=r0*ecc
      cosv=prscal(x,vlenz)/dd  ! cosine of true anomaly
      CALL prvec(vlenz,x,ww)
      sinv=prscal(ang,ww)/dd   ! sine of true anomaly   
      eq(6)=atan2(sinv,cosv)
      IF(eq(6).lt.0.d0) eq(6)=eq(6)+dpig
      to_elems%coo='COT'
   ENDIF
ELSE
   eq(1)=-gm/alpha
! selection of different methods
   IF(ecc.gt.eccbou)THEN
      peri=eq(1)*(1-ecc)
      enne=dsqrt(gm/eq(1)**3)
      CALL solve_peri(r0,sig0,peri,gm,alpha,psi,dt)
      eq(6)=princ(-dt*enne+varpi)
   ELSE
! mean longitude from non--singular Kepler equation                   
      rad=dsqrt(1.d0-ecc2) 
      beta=1.d0/(1.d0+rad) 
      chk=eq(2)*eq(3)*beta 
      ch=1.d0-eq(2)**2*beta 
      ck=1.d0-eq(3)**2*beta 
      x2=prscal(x,fb) 
      y2=prscal(x,gb) 
      cosf=eq(3)+(ck*x2-chk*y2)/(eq(1)*rad) 
      sinf=eq(2)+(ch*y2-chk*x2)/(eq(1)*rad) 
! eccentric longitude                                                 
      fe=datan2(sinf,cosf) 
      eq(6)=fe+eq(2)*cosf-eq(3)*sinf 
! reduction to principal value  
      eq(6)=princ(eq(6))
   ENDIF
   IF(cooy.eq.'KEP'.and.fail_flag.eq.0)THEN
      eq(2)=ecc
      eq(4)=omnod
      eq(3)=2*atan(tgi2)
      eq(5)=princ(varpi-omnod)  
      eq(6)=princ(eq(6)-varpi)
      to_elems%coo='KEP'
   ELSEIF(cooy.eq.'EQU'.or.(fail_flag.ge.1.and.fail_flag.le.3))THEN
!          coord(2)= h = e sin(varpi) ; varpi = Omega + omega
!          coord(3)= k = e cos(varpi)
!          coord(4)= p = tg(I/2) sin(Omega)
!          coord(5)= q = tg(I/2) cos(Omega)
!          coord(6)= mean longitude, measured from fb vector (mean an + varpi)
      fail_flag=0
      to_elems%coo='EQU'
   ELSE
      PRINT *, ' to_elems: this should not happen ', cooy, fail_flag
   ENDIF
ENDIF
to_elems%coord=eq
IF(PRESENT(DEL))THEN
! compute jacobian matrix by using the inverse transformation
! this is dirty, but the derivatives of elems w.r. to cartesian
! are not as important as the inverse
elback=car_elems(to_elems,fail_flag1,del1)
! WARNING: how to use fail_flag1??????
! inversion
CALL matin(del1,det,6,0,6,ising,1)
IF(ising.eq.0)THEN
   del=del1
ELSE
   del=0.d0
! WARNING: how to set fail_flag??? as above????
ENDIF
ENDIF
END FUNCTION to_elems
! ====================================================================
! CAR_ELEMS
! ====================================================================
TYPE(orbit_elem) FUNCTION car_elems(el,fail_flag,del)
USE ever_pitkin
TYPE(orbit_elem),INTENT(IN):: el ! can be 'KEP', 'EQU', 'COM', 'COT'
INTEGER, INTENT(OUT) :: fail_flag
! optional jacobian matrix. in output: derivatives d(car)/d(el)
DOUBLE PRECISION, DIMENSION(6,6),INTENT(OUT),OPTIONAL :: del
! END INTERFACE
DOUBLE PRECISION :: gm ! center body mass
CHARACTER*(3) :: coox
TYPE(orbit_elem) :: equ ! used for equinoctal, to remove the KEP case 
DOUBLE PRECISION, DIMENSION(6,6) :: del1,del0,del2,del3
DOUBLE PRECISION :: eq(6) ! coordinate 6d-vectors
DOUBLE PRECISION, DIMENSION(3) :: x, y, x0, y0, acc
DOUBLE PRECISION :: princ, vsize, eps ! functions, small no
DOUBLE PRECISION :: r0, v0, dt
DOUBLE PRECISION, DIMENSION(3) :: fb,gb,dfdp,dfdq,dgdp,dgdq ! Brouke's vectors
DOUBLE PRECISION :: tgim2, tgim, ecc, upq, ecc2, varpi, enne, mean_an
! vectors used in computing derivatives (not the complete derivatives)
DOUBLE PRECISION :: dedh, dedk, ddigdh, ddigdk, fact, rsvperi,q, upecv2
DOUBLE PRECISION cosn,sinn,coso,sino,cosi,sini, cosv,sinv,rsq,fgfact,eye3(3,3)
! paranoia check of input arguments
coox=el%coo
IF(coox.ne.'COM'.and.coox.ne.'KEP'.and.coox.ne.'EQU'.and.coox.ne.'COT')THEN
   WRITE(*,*)' car_elems: wrong input coordinates', coox
   fail_flag=20
   RETURN
ENDIF
! initializations
fail_flag=0
gm=centerbody_mass(el%center)
eps=100*epsilon(1.d0) ! small number
car_elems=el ! copy other data (e.g., magnitude)
! ================================================================
! the position and velocity at pericenter are computed for all
!  output coordiantes
! ================================================================ 
IF(coox.eq.'COM'.or.coox.eq.'COT')THEN
   ecc=el%coord(2)
   ecc2=ecc**2
   IF(ecc.lt.eps)THEN
      varpi=0.d0
      fail_flag=1
   ELSE
      varpi=el%coord(4)+el%coord(5)
   ENDIF   
   IF(el%coord(3).lt.eps) fail_flag=fail_flag+2
! pericenter  
   r0=el%coord(1)
   v0=sqrt(gm*(1.d0+ecc)/el%coord(1))
!   WRITE(*,*) ' r0,v0 ',r0,v0
   cosn=cos(el%coord(4))
   sinn=sin(el%coord(4))
   coso=cos(el%coord(5))
   sino=sin(el%coord(5))
   cosi=cos(el%coord(3))
   sini=sin(el%coord(3))
   x0=(/coso*cosn-sino*sinn*cosi, coso*sinn+sino*cosn*cosi, sino*sini/)*r0
   y0=(/-sino*cosn-coso*sinn*cosi, -sino*sinn+coso*cosn*cosi, coso*sini/)*v0
   IF(coox.eq.'COM')THEN
      dt=el%t-el%coord(6) ! time between peihelion and epoch, for 2-body propagation
   ELSEIF(coox.eq.'COT')THEN
      cosv=cos(el%coord(6)) ! cosine and sine of true anomaly, for 2-body formulae
      sinv=sin(el%coord(6))
   ENDIF
ELSE
   IF(el%coo.eq.'KEP')THEN
      IF(PRESENT(del))THEN
         equ=equ_kep(el,del1)
      ELSE
         equ=equ_kep(el)
      ENDIF
   ELSEIF(el%coo.eq.'EQU')THEN
      equ=el
      IF(PRESENT(del))CALL eye(6,del1)
   ENDIF
   eq=equ%coord
   ecc2=eq(2)**2+eq(3)**2
   ecc=sqrt(ecc2)
! pericenter
   r0=eq(1)*(1.d0-ecc)
   v0=sqrt(gm*(1.d0+ecc)/(eq(1)*(1.d0-ecc)))
   enne=sqrt(gm/eq(1)**3)
!   WRITE(*,*) ' r0,v0, enne ',r0,v0,enne
! varpi= Omega+ omega
   IF(ecc.lt.eps)THEN
      varpi=0.d0
      fail_flag=1
   ELSE
      varpi=atan2(eq(2),eq(3))
   ENDIF
   mean_an=princ(eq(6)-varpi)
   IF(mean_an.gt.pig)mean_an=mean_an-dpig
   dt=mean_an/enne
   tgim2=eq(4)**2+eq(5)**2 
   tgim=dsqrt(tgim2)
   IF(tgim.lt.eps) fail_flag=fail_flag+2                                     
   upq=1.d0+tgim2 
!  unit vectors of the equinoctal reference system (Broucke and         
!  Cefola 1972, CM 5, 303--310) are f, g and the angular momentum       
!  unit vector                                                          
   fb(1)=(1.d0-eq(4)**2+eq(5)**2)/upq 
   fb(2)=2*eq(4)*eq(5)/upq 
   fb(3)=-2*eq(4)/upq 
   gb(1)=2*eq(4)*eq(5)/upq 
   gb(2)=(1.d0+eq(4)**2-eq(5)**2)/upq 
   gb(3)=2*eq(5)/upq    
! pericenter vector, parallel to Lenz vector
   x0=r0*(fb*cos(varpi)+gb*sin(varpi))
! velocity at pericenter
   y0=v0*(-fb*sin(varpi)+gb*cos(varpi))
ENDIF
! WARNING: for now the use of f-g series is avoided only for coox='COT'
IF(coox.eq.'COT')THEN
   q=el%coord(1)
   rsq=(1.d0+ecc)/(1.d0+ecc*cosv) ! r/q
   rsvperi=sqrt((1.d0+ecc)*q**3/gm)/(1+ecc*cosv)
   x=x0*(rsq*cosv)+y0*(rsvperi*sinv)
   fgfact=sqrt(gm/(q**3*(1.d0+ecc)))
   y=(-sinv*fgfact)*x0+((cosv+ecc)/(1.d0+ecc))*y0
ELSE
   IF(PRESENT(del))THEN
      CALL fser_propag_der(x0,y0,0.d0,dt,gm,x,y,del3)
   ELSE
      CALL fser_propag(x0,y0,0.d0,dt,gm,x,y)
   ENDIF
ENDIF
car_elems%coord(1:3)=x
car_elems%coord(4:6)=y
car_elems%coo='CAR'
IF(PRESENT(del))THEN

! compute partials of x0, y0 w.r. to elements eq(1:5)

   IF(coox.eq.'COM'.or.coox.eq.'COT')THEN
      del0=0.d0
      del0(1:3,1)=  x0*(1/el%coord(1))
      del0(4:6,1)=  y0*(-0.5/el%coord(1))
      del0(1:3,2)=  0.d0
      del0(4:6,2)=  y0*0.5/(1.d0+ecc)
      del0(1:3,3)=  (/ sinn*sino*sini, -cosn*sino*sini, sino*cosi/)*r0
      del0(4:6,3)=  (/ sinn*coso*sini, -cosn*coso*sini, coso*cosi/)*v0
      del0(1:3,4)=  (/-sinn*coso-cosn*sino*cosi, cosn*coso-sinn*sino*cosi, 0.d0/)*r0
      del0(4:6,4)=  (/sinn*sino-cosn*coso*cosi, -cosn*sino-sinn*coso*cosi, 0.d0/)*v0
      fact= r0/v0
      del0(1:3,5)=  fact*y0
      del0(4:6,5)= -(1.d0/fact)*x0  
      IF(coox.eq.'COM')THEN    
! now use Goodyear to compute jacobian of propagation from t=eq(6) to t=eqn%t
         del=MATMUL(del3,del0)
! then fix derivatives with respect to time of pericenter
         acc=-x*gm/vsize(x)**3
         del(1:3,6)=  -y
         del(4:6,6)=  -acc
      ELSEIF(coox.eq.'COT')THEN
         del3=0.d0 ! del3=\partial (x,y)/\partial (x0,y0)
         CALL eye(3,eye3)
         del3(1:3,1:3)=rsq*cosv*eye3
         del3(1:3,4:6)=rsvperi*sinv*eye3
         del3(4:6,1:3)=-fgfact*sinv*eye3
         del3(4:6,4:6)=((cosv+ecc)/(1.d0+ecc))*eye3
         del2=0.d0 ! del2= \partial (x,y)/\partial el%coord (not through x0,y0)
         upecv2=(1.d0+ecc*cosv)**(-2)
         del2(1:3,1)=1.5d0*sqrt(q*(1.d0+ecc)/gm)*sinv/(1+ecc*cosv) * y0 ! \partial x/\partial q
         del2(1:3,2)=cosv*(1.d0-cosv)*upecv2 * x0    &! \partial x/\partial ecc
&            +0.5d0*sqrt(q**3/((1.d0+ecc)*gm))*sinv*(1.d0-ecc*cosv-2.d0*cosv)*upecv2 * y0
         del2(1:3,6)=-sinv*(1.d0+ecc)*upecv2 * x0  & ! \partial x/\partial v
&             +sqrt(q**3*(1.d0+ecc)/gm)*(cosv+ecc)*upecv2* y0
         del2(4:6,1)=1.5d0*sqrt(gm/(q**5*(1.d0+ecc)))*sinv * x0 ! \partial y/\partial q
         del2(4:6,2)=0.5d0*sqrt(gm/(q**3*(1.d0+ecc)**3))*sinv * x0  & ! \partial y/\partial e
&             + (1.d0-cosv)/(1.d0+ecc)**2 * y0 
         del2(4:6,6)= -sqrt(gm/(q**3*(1.d0+ecc)))*cosv * x0 - sinv/(1.d0+ecc) * y0 ! \partial y/\partial v
         del=MATMUL(del3,del0)+del2 ! chain rule
      ENDIF
   ELSEIF(coox.eq.'EQU'.or.coox.eq.'KEP')THEN
      del2=0.d0
      del2(1:3,1)=  x0*(1/eq(1))
      del2(4:6,1)=  y0*(-0.5/eq(1))
      dedh=eq(2)/ecc
      dedk=eq(3)/ecc
      ddigdh=eq(3)/ecc2
      ddigdk=-eq(2)/ecc2
      fact=r0/v0
      del2(1:3,2)= -dedh/(1.d0-ecc)*x0+fact*ddigdh*y0
      del2(1:3,3)= -dedk/(1.d0-ecc)*x0+fact*ddigdk*y0
!      del2(1:3,2)=  eq(1)*(-gb+ddigdh*(-sin(varpi)*fb+cos(varpi)*gb))
      del2(1:3,3)= -dedk/(1.d0-ecc)*x0+fact*ddigdk*y0
!      del2(1:3,3)=  eq(1)*(-fb+ddigdk*(-sin(varpi)*fb+cos(varpi)*gb))
      del2(4:6,2)=  dedh/(1.d0-ecc2)*y0-(1.d0/fact)*ddigdh*x0
      del2(4:6,3)=  dedk/(1.d0-ecc2)*y0-(1.d0/fact)*ddigdk*x0
      dfdp(1)=-2*eq(4)
      dfdp(2)=2*eq(5)
      dfdp(3)=-2.d0
      dfdp=(dfdp-fb*(2*eq(4)))/upq
      dfdq(1)=2*eq(5)
      dfdq(2)=2*eq(4)
      dfdq(3)=0.d0
      dfdq=(dfdq-fb*(2*eq(5)))/upq
      dgdp(1)=2*eq(5)
      dgdp(2)=2*eq(4)
      dgdp(3)=0.d0
      dgdp=(dgdp-gb*(2*eq(4)))/upq
      dgdq(1)=2*eq(4)
      dgdq(2)=-2*eq(5)
      dgdq(3)=2.d0
      dgdq=(dgdq-gb*(2*eq(5)))/upq
      del2(1:3,4)=  r0*(cos(varpi)*dfdp+sin(varpi)*dgdp)
      del2(1:3,5)=  r0*(cos(varpi)*dfdq+sin(varpi)*dgdq)
      del2(4:6,4)=  v0*(-sin(varpi)*dfdp+cos(varpi)*dgdp)
      del2(4:6,5)=  v0*(-sin(varpi)*dfdq+cos(varpi)*dgdq)
      del2(1:6,6)=  0.d0
      del0=MATMUL(del2,del1)
! now use Goodyear to compute jacobian of propagation from t=eq(6) to t=eqn%t
      del=MATMUL(del3,del0)
! then fix derivatives with respect to mean anomaly
      acc=-x*gm/vsize(x)**3
      del(1:3,6)=y/enne
      del(4:6,6)=acc/enne
! then correct derivatives with respect to a to account for enne in dt
      del(1:3,1)= del(1:3,1)+(1.5d0*mean_an/(enne*eq(1)))*y
      del(4:6,1)= del(4:6,1)+(1.5d0*mean_an/(enne*eq(1)))*acc
! then correct derivatives with respedct to h,k to account for varpi in dt
      IF(coox.eq.'EQU')THEN
         del(1:3,2)=del(1:3,2)-(ddigdh/enne)*y
         del(1:3,3)=del(1:3,3)-(ddigdk/enne)*y
         del(4:6,2)=del(4:6,2)-(ddigdh/enne)*acc
         del(4:6,3)=del(4:6,3)-(ddigdk/enne)*acc
      ENDIF
   ENDIF
ENDIF

END FUNCTION car_elems
! ====================================================================
! ELEMENTS INTERNAL CONVERSIONS
! ====================================================================
! KEP_EQU
! keplerian elements from equinoctal
! =======================================================
TYPE(orbit_elem) FUNCTION kep_equ(el,fail_flag,del)
TYPE(orbit_elem),INTENT(IN):: el ! must be equinoctal
INTEGER,INTENT(OUT):: fail_flag
! optional jacobian matrix. in output: derivatives d(el)/d(eqn)
DOUBLE PRECISION, DIMENSION(6,6),INTENT(OUT),OPTIONAL :: del
! end interface
DOUBLE PRECISION dig,ecc,ecc2,tgi2,tgi2s,princ,eps
DOUBLE PRECISION,DIMENSION(6) :: eq
!
IF(el%coo.ne.'EQU')THEN
   WRITE(*,*)' kep_equ: wrong input coordinates',el
   STOP
ENDIF
fail_flag=0
eps=100*epsilon(1.d0) ! negligible number
eq=el%coord
kep_equ=el  !default
kep_equ%coord(1)=eq(1) !semimajor axis
!  test on eccentricity                                                 
ecc2=eq(2)**2+eq(3)**2
ecc=sqrt(ecc2) 
if(ecc.lt.eps)then 
   fail_flag=1 ! but computation is completed
   dig=0.d0 
else 
   dig=atan2(eq(2),eq(3)) 
endif
kep_equ%coord(2)=ecc 
!   test on tangent of half inclination                                 
tgi2s=eq(4)**2+eq(5)**2
tgi2=sqrt(tgi2s) 
if(tgi2.lt.eps)then 
   fail_flag=fail_flag+2 ! but computation is completed
   kep_equ%coord(4)=0.d0 
else 
   kep_equ%coord(4)=atan2(eq(4),eq(5)) 
endif
kep_equ%coord(3)=2.d0*atan(tgi2) 
!   angular variables                                                   
kep_equ%coord(5)=dig-kep_equ%coord(4) 
kep_equ%coord(6)=eq(6)-dig 
kep_equ%coord(4)=princ(kep_equ%coord(4)) 
kep_equ%coord(5)=princ(kep_equ%coord(5)) 
kep_equ%coord(6)=princ(kep_equ%coord(6)) 
kep_equ%coo='KEP'
IF(PRESENT(del).and.fail_flag.eq.0)THEN
   del=0.d0 ! jacobian is not usable in the nearly singular case
   del(1,1)=  1.d0 
   del(2,2)=  eq(2)/ecc ! ecc=h^2+k^2
   del(2,3)=  eq(3)/ecc
   del(3,4)=  2*eq(4)/((1+tgi2**2)*tgi2) 
   del(3,5)=  2*eq(5)/((1+tgi2**2)*tgi2)
   del(4,4)=  eq(5)/tgi2s ! Omega=atan2(p,q)
   del(4,5)= -eq(4)/tgi2s
   del(5,2)=  eq(3)/ecc2  ! dig=atan2(h,k)
   del(5,3)= -eq(2)/ecc2  ! omega=dig-Omega 
   del(5,4)= -eq(5)/tgi2s
   del(5,5)=  eq(4)/tgi2s
   del(6,2)= -eq(3)/ecc2  ! ell=lambda-dig
   del(6,3)=  eq(2)/ecc2
   del(6,6)=  1.d0
! else the jacobian is not defined!!! 
ENDIF

RETURN
END FUNCTION kep_equ

! =======================================================
! EQU_KEP
! equinoctal elements from keplerian
! =======================================================
TYPE(orbit_elem) FUNCTION equ_kep(el,del)
TYPE(orbit_elem),INTENT(IN):: el ! must be keplerian
! optional jacobian matrix. in output: derivatives d(el)/d(eqn)
DOUBLE PRECISION, DIMENSION(6,6),INTENT(OUT),OPTIONAL :: del
! end interface
DOUBLE PRECISION dig,ecc,tgim,princ, inc
DOUBLE PRECISION sindig,cosdig,sinom,cosom,dti2
!
IF(el%coo.ne.'KEP')THEN
   WRITE(*,*)' equ_kep: wrong input coordinates',el
   STOP
ENDIF 
equ_kep=el  !default
equ_kep%coord(1)=el%coord(1) 
dig=el%coord(4)+el%coord(5) 
ecc=el%coord(2) 
sindig=sin(dig)
cosdig=cos(dig)
equ_kep%coord(2)=ecc*sindig 
equ_kep%coord(3)=ecc*cosdig
IF(el%coord(3).ge.pig)THEN
   inc=pig-1.d-4
ELSE
   inc= el%coord(3)
ENDIF
tgim=tan(inc/2.d0) 
sinom=sin(el%coord(4))
cosom=cos(el%coord(4))
equ_kep%coord(4)=tgim*sinom
equ_kep%coord(5)=tgim*cosom 
equ_kep%coord(6)=dig+el%coord(6) 
equ_kep%coord(6)=princ(equ_kep%coord(6))
equ_kep%coo='EQU' 
IF(PRESENT(del))THEN
   del=0.d0
   del(1,1)=1.d0 
   del(2,2)=sindig
   del(2,5)=equ_kep%coord(3)
   del(2,4)=equ_kep%coord(3) 
   del(3,2)=cosdig
   del(3,5)=-equ_kep%coord(2)
   del(3,4)=-equ_kep%coord(2)
   dti2=1.d0/(2*cos(el%coord(3)/2)**2) 
   del(4,3)=dti2*sinom 
   del(4,4)= equ_kep%coord(5)
   del(5,3)=dti2*cosom 
   del(5,4)=- equ_kep%coord(4)
   del(6,4)=1.d0 
   del(6,5)=1.d0 
   del(6,6)=1.d0 
ENDIF
RETURN
END FUNCTION equ_kep
! =======================================================
! COM_KEP
! cometary elelements from keplerian
! =======================================================
TYPE(orbit_elem) FUNCTION com_kep(el,fail_flag,del)
TYPE(orbit_elem),INTENT(IN):: el ! must be keplerian
INTEGER,INTENT(OUT):: fail_flag
! optional jacobian matrix. in output: derivatives d(el)/d(eqn)
DOUBLE PRECISION, DIMENSION(6,6),INTENT(OUT),OPTIONAL :: del
! end interface
DOUBLE PRECISION ecc, eps, gm, enne, princ
DOUBLE PRECISION,DIMENSION(6) :: eq
!
IF(el%coo.ne.'KEP')THEN
   WRITE(*,*)' com_kep: wrong input coordinates',el
   STOP
ENDIF
fail_flag=0
gm=centerbody_mass(el%center)
eps=100*epsilon(1.d0) ! negligible number
eq=el%coord
com_kep=el  !default
ecc=eq(2)
com_kep%coord(1)=eq(1)*(1.d0-ecc)
!  test on eccentricity                                                 
if(ecc.lt.eps) fail_flag= 1 ! but computation is completed
!   test on tangent of half inclination                                 
if(eq(3).lt.eps) fail_flag=fail_flag+2 ! but computation is completed
! time of pericenter from mean anomaly
enne=sqrt(gm/eq(1)**3)
eq(6)=princ(eq(6))
com_kep%coord(6)=el%t-eq(6)/enne 
com_kep%coo='COM'
IF(PRESENT(del).and.fail_flag.eq.0)THEN
   del=0.d0 ! jacobian is not usable in the nearly singular case
   del(1,1)=  1.d0-ecc
   del(1,2)=  -eq(1) 
   del(2,2)=  1.d0
   del(3,3)=  1.d0
   del(4,4)=  1.d0
   del(5,5)=  1.d0
   del(6,1)= -1.5d0*eq(6)/(enne*eq(1))
   del(6,6)=  -1.d0/enne
ENDIF
END FUNCTION com_kep

! =======================================================
! COT_KEP
! cometary elements with true anomaly from keplerian
! =======================================================
TYPE(orbit_elem) FUNCTION cot_kep(el,fail_flag,del)
TYPE(orbit_elem),INTENT(IN):: el ! must be keplerian
INTEGER,INTENT(OUT):: fail_flag
! optional jacobian matrix. in output: derivatives d(el)/d(eqn)
DOUBLE PRECISION, DIMENSION(6,6),INTENT(OUT),OPTIONAL :: del
! end interface
DOUBLE PRECISION ecc, eps, gm, enne, princ
DOUBLE PRECISION, DIMENSION(6) :: eq
DOUBLE PRECISION ell,u,du,cosu,sinu,cosv,sinv,dvdu,dvde,dudell,dude 
INTEGER, PARAMETER :: itmax=10
INTEGER i
!
IF(el%coo.ne.'KEP')THEN
   WRITE(*,*)' cot_kep: wrong input coordinates',el
   STOP
ENDIF
fail_flag=0
gm=centerbody_mass(el%center)
eps=100*epsilon(1.d0) ! negligible number
eq=el%coord
cot_kep=el  !default
ecc=eq(2)
cot_kep%coord(1)=eq(1)*(1.d0-ecc)
!  test on eccentricity                                                 
if(ecc.lt.eps) fail_flag= 1 ! but computation is completed
!   test on tangent of half inclination                                 
if(eq(3).lt.eps) fail_flag=fail_flag+2 ! but computation is completed
! time of pericenter from mean anomaly
!enne=sqrt(gm/eq(1)**3)
! Kepler's equation to find eccentric anomaly
u=pig
DO i=1,itmax
   cosu=cos(u);sinu=sin(u)
   ell=u-ecc*sinu
   du=(eq(6)-ell)/(1.d0-ecc*cosu)
   u=u+du
   IF(abs(du).lt.eps*10) GOTO 3
ENDDO
WRITE(ierrou,*)' cot_kep: non convergent Kepler equation, u,e,du:', u,ecc,du
3 CONTINUE
cosv=(cosu-ecc)/(1.d0-ecc*cosu)
sinv=sqrt(1.d0-ecc**2)*sinu/(1.d0-ecc*cosu)
cot_kep%coord(6)=atan2(sinv,cosv)
cot_kep%coord(6)=princ(cot_kep%coord(6))
cot_kep%coo='COT'
IF(PRESENT(del).and.fail_flag.eq.0)THEN
   del=0.d0 ! jacobian is not usable in the nearly singular case
   del(1,1)=  1.d0-ecc
   del(1,2)=  -eq(1) 
   del(2,2)=  1.d0
   del(3,3)=  1.d0
   del(4,4)=  1.d0
   del(5,5)=  1.d0
   dvdu=sqrt((1.d0+ecc)/(1.d0-ecc))*(cosv+1.d0)/(cosu+1.d0)
   dvde=tan(u/2.d0)*(cosv+1.d0)/sqrt((1.d0+ecc)*(1.d0-ecc)**3)
   dudell=1.d0/(1.d0-ecc*cosu)
   dude=sinu/(1.d0-ecc*cosu)
   del(6,6)=dvdu*dudell
   del(6,2)=dvde+dvdu*dude
ENDIF
END FUNCTION cot_kep
! =======================================================
! KEP_COT
! keplerian elelements from cometary with true anomaly
! =======================================================
TYPE(orbit_elem) FUNCTION kep_cot(el,fail_flag,del)
TYPE(orbit_elem),INTENT(IN):: el ! must be cometary
INTEGER,INTENT(OUT):: fail_flag
! optional jacobian matrix. in output: derivatives d(el)/d(eqn)
DOUBLE PRECISION, DIMENSION(6,6),INTENT(OUT),OPTIONAL :: del
! end interface
DOUBLE PRECISION ecc, eps, gm, enne, a, dldn
DOUBLE PRECISION,DIMENSION(6) :: eq
DOUBLE PRECISION r,x1,y1,sinu,cosu,u,delldu,dellde,dudv,dude,v2
!
IF(el%coo.ne.'COT')THEN
   WRITE(*,*)' kep_cot: wrong input coordinates',el
   STOP
ENDIF
fail_flag=0
gm=centerbody_mass(el%center)
eps=100*epsilon(1.d0) ! negligible number
eq=el%coord
kep_cot=el  !default
ecc=eq(2)
if(ecc.ge.1.d0-eps)THEN
! this conversion is impossible
   fail_flag=5
! in output the elements are left as 'COT', del= identity
   IF(PRESENT(del)) CALL eye(6,del)
   RETURN
ENDIF
a=eq(1)/(1.d0-ecc)
kep_cot%coord(1)=a
!  test on eccentricity                                                 
if(ecc.lt.eps) fail_flag= 1 ! but computation is completed
!   test on tangent of half inclination                                 
if(eq(3).lt.eps) fail_flag=fail_flag+2 ! but computation is completed
! eccentric anomaly from true anomaly
r=a*(1.d0-ecc**2)/(1.d0+ecc*cos(eq(6)))
x1=r*cos(eq(6)); y1=r*sin(eq(6))
sinu=y1/(a*sqrt(1.d0-ecc**2)); cosu=x1/a+ecc
u=atan2(sinu,cosu)
! kepler's equation (easy way)
kep_cot%coord(6)=u-ecc*sinu
kep_cot%coo='KEP' ! successful conversion
IF(PRESENT(del).and.fail_flag.eq.0)THEN
   del=0.d0 ! jacobian is not usable in the nearly singular case
   del(1,1)=  1.d0/(1.d0-ecc)
   del(1,2)=  +eq(1)/(1.d0-ecc)**2 
   del(2,2)=  1.d0
   del(3,3)=  1.d0
   del(4,4)=  1.d0
   del(5,5)=  1.d0
! derivatives of mean anomaly (composite through eccentric anomaly)
   delldu=1.d0-ecc*cosu
   dellde=-sinu
   v2=eq(6)/2.d0
   dudv=sqrt((1.d0-ecc)/(1.d0+ecc))*cos(u/2.d0)**2/cos(v2)**2
   dude=(-2*tan(v2)*cos(u/2.d0)**2)/sqrt((1.d0-ecc)*(1.d0+ecc)**3)
   del(6,6)=delldu*dudv
   del(6,2)=dellde + delldu*dude
ENDIF
END FUNCTION kep_cot
! =======================================================
! COT_COM
! cometary elements with true anomaly from cometary
! =======================================================
TYPE(orbit_elem) FUNCTION cot_com(el,fail_flag,del)
TYPE(orbit_elem),INTENT(IN):: el ! must be keplerian
INTEGER,INTENT(OUT):: fail_flag
! optional jacobian matrix. in output: derivatives d(el)/d(eqn)
DOUBLE PRECISION, DIMENSION(6,6),INTENT(OUT),OPTIONAL :: del
! end interface
DOUBLE PRECISION ecc, eps , gm
DOUBLE PRECISION delcar(6,6),dcarel(6,6)
TYPE(orbit_elem) elcar,elcot
DOUBLE PRECISION dt, alpha,v,dvdq,dvde,dvdt0, vel2
INTEGER fail_flag1, fail_flag2
! ======================================================
IF(el%coo.ne.'COM')THEN
   WRITE(*,*)' cot_com: wrong input coordinates',el
   STOP
ENDIF
fail_flag=0
gm=centerbody_mass(el%center)
eps=100*epsilon(1.d0) ! negligible number
! eq=el%coord
cot_com=el  !default
ecc=el%coord(2)
!  test on eccentricity                                                 
if(ecc.lt.eps) fail_flag= 1 ! but computation is completed
!   test on tangent of half inclination                                 
if(el%coord(3).lt.eps) fail_flag=fail_flag+2 ! but computation is completed
! find true anomaly from dt
dt=el%t-el%coord(6)
alpha=gm*(el%coord(2)-1.d0)/el%coord(1) ! 2* energy
! CALL ta_from_t0(el%coord(1),ecc,dt,gm,alpha,v,dvdq,dvde,dvdt0)
! cot_com%coord(6)=v
cot_com%coo='COT'
! use brutal composition of two available changes
IF(PRESENT(del))THEN
   CALL coo_cha(el,'CAR',elcar,fail_flag1,delcar)
   CALL coo_cha(elcar,'COT',elcot, fail_flag2,dcarel)
   del=MATMUL(dcarel,delcar)
ELSE
   CALL coo_cha(el,'CAR',elcar,fail_flag1)
   CALL coo_cha(elcar,'COT',elcot, fail_flag2)
ENDIF 
cot_com%coord(6)=elcot%coord(6)
! cot_com%coord(6)=v
IF(PRESENT(del))THEN
  CALL eye(5,del(1:5,1:5))
!  del(6,1)=dvdq ! NOT READY
!  del(6,2)=dvde
  del(6,3:5)=0.d0
!  del(6,6)=dvdt0
  del(1:5,6)=0.d0
ENDIF
RETURN
! WARNING: what to do with fail_flag1 and fail_flag2 ???
END FUNCTION cot_com
! ======================================================
! TA_FROM_T0 true anomaly from time of perihelion
! ======================================================
SUBROUTINE ta_from_t0(q,ecc,dt,mu,alpha,v,dvdq,dvde,dvdt0)
  USE ever_pitkin
  DOUBLE PRECISION,INTENT(IN) :: q,ecc ! perihelion distance, eccentricity
  DOUBLE PRECISION,INTENT(IN) :: dt ! dt=t-t0 difference between 
                       ! present time and time of passage at perihelion
  DOUBLE PRECISION,INTENT(IN) :: mu, alpha ! grav. constant, 2*energy
  DOUBLE PRECISION,INTENT(OUT) :: v ! true anomaly
  DOUBLE PRECISION,INTENT(OUT) :: dvdq,dvde,dvdt0 ! derivatives of true anomaly
! end interface                                   ! w.r.t. q,e,t0
  DOUBLE PRECISION :: sig0, rdot, psi, s0,s1,s2,s3, r
  DOUBLE PRECISION :: cosv,dvdr, dpsidt0, drdt0,drdpsi, ds2da,ds0da, drda
  DOUBLE PRECISION s0v,s1v,s2v,s3v,dd, psiv, eccv,qv, vv,cosvv
! *******************************
! compute psi at time t (psi0=0)
! *******************************
  sig0 = 0.d0
! solve universal Kepler's equation and compute Stumpff's functions
  CALL solve_kepuniv2(dt,q,sig0,mu,alpha,ecc,psi,s0,s1,s2,s3,DS2DA=ds2da,DS0DA=ds0da)
  dd=1.d-12
!  CALL s_funct(psi,alpha+dd,s0v,s1v,s2v,s3v)
!  CALL solve_kepuniv2(dt,q,sig0,mu,alpha+dd,ecc,psiv,s0v,s1v,s2v,s3v)
!  WRITE(*,*)' ds2da=',(s2v-s2)/dd, ds2da, s2
!  WRITE(*,*)' ds0da=',(s0v-s0)/dd, ds0da, s0
!  WRITE(*,*)' psi=', psi, psiv,psi-psiv
!  ds2da=(s2v-s2)/dd
!  ds0da=(s0v-s0)/dd
! *******************************
! compute r at time t 
! *******************************
  r = q*s0 + sig0*s1 + mu*s2
  drdpsi = sig0*s0 + (mu+alpha*q)*s1
  cosv =(-1.d0 + q/r*(1.d0+ecc))/ecc 
  v=acos(cosv)*sign(1.d0,drdpsi)

  qv=q+dd
  cosvv =(-1.d0 + qv/r*(1.d0+ecc))/ecc 
  vv=acos(cosvv)*sign(1.d0,drdpsi)
  dvdq=(vv-v)/dd
  WRITE(*,*)' dvdq direct inc=', dvdq
  eccv=ecc+dd
  cosvv =(-1.d0 + q/r*(1.d0+eccv))/eccv 
  vv=acos(cosvv)*sign(1.d0,drdpsi)
  dvde=(vv-v)/dd
  WRITE(*,*)' dvde direct inc=', dvde

! derivatives of v w.r. to q,e 
! ***********WARNING: WRONG***************
! missing dependence of r from q,e 
! dr/dq=s0+q*d(s0)/dq+ mu*d(s2)/dq 
! dr/d(alpha)=q*d(s0)/d(alpha)+ mu*d(s2)/d(alpha)
! d(sj)/dq=d(sj)/d(alpha) d(alpha)/dq ; d(sj)/de=d(sj)/d(alpha) d(alpha)/de
! d(alpha)/dq=-alpha/q ; d(alpha)/de= mu/q
  dvdq = -1.d0/sqrt(1.d0-cosv**2)*(1.d0+ecc)/(r*ecc) *sign(1.d0,drdpsi)
  WRITE(*,*)' dvdq direct=', dvdq
  dvdr =  1.d0/sqrt(1.d0-cosv**2)*q*(1.d0+ecc)/(r**2*ecc)*sign(1.d0,drdpsi)! OK
  drda= q*ds0da + mu*ds2da
!  drdq=s0+(q*ds0da+mu*ds2da)*(-alpha/q)
  dvdq=dvdq + dvdr*(drda*(-alpha/q)+s0)
  dvde = -1.d0/sqrt(1.d0-cosv**2)*(1.d0-q/r)/ecc**2*sign(1.d0,drdpsi)
  WRITE(*,*)' dvde direct=', dvde
  dvde= dvde + dvdr*drda*mu/q
! ************END WARNING*****************
! derivative of v w.r. to t0
  dpsidt0 = -1.d0/(q*s0 + sig0*s1 + mu*s2)
  drdt0 = drdpsi*dpsidt0
  dvdt0 = dvdr*drdt0
END SUBROUTINE ta_from_t0

! =======================================================
! COM_COT
! cometary elelements from cometary with true anomaly
! =======================================================
TYPE(orbit_elem) FUNCTION com_cot(el,fail_flag,del)
TYPE(orbit_elem),INTENT(IN):: el ! must be cometary
INTEGER,INTENT(OUT):: fail_flag
! optional jacobian matrix. in output: derivatives d(el)/d(eqn)
DOUBLE PRECISION, DIMENSION(6,6),INTENT(OUT),OPTIONAL :: del
! end interface
DOUBLE PRECISION ecc, eps, gm !, enne, a, dldn, dt
!DOUBLE PRECISION,DIMENSION(6) :: eq
!DOUBLE PRECISION r,x1,y1,sinu,cosu,u,delldu,dellde,dudv,dude,v2,ell,delldv
DOUBLE PRECISION delcar(6,6),dcarel(6,6)
TYPE(orbit_elem) elcar,elcom
INTEGER fail_flag1, fail_flag2
! =======================================================
IF(el%coo.ne.'COT')THEN
   WRITE(*,*)' com_cot: wrong input coordinates',el
   STOP
ENDIF
fail_flag=0
eps=100*epsilon(1.d0) ! negligible number
com_cot=el  !default
ecc=el%coord(2) ! eccentricity
!  test on eccentricity                                                 
IF(ecc.lt.eps) fail_flag= 1 ! but computation is completed
!   test on inclination                                 
IF(el%coord(3).lt.eps) fail_flag=fail_flag+2 ! but computation is completed
! use brutal composition of two available changes
IF(PRESENT(del))THEN
   CALL coo_cha(el,'CAR',elcar,fail_flag1,delcar)
   CALL coo_cha(elcar,'COM',elcom, fail_flag2,dcarel)
   del=MATMUL(dcarel,delcar)
ELSE
   CALL coo_cha(el,'CAR',elcar,fail_flag1)
   CALL coo_cha(elcar,'COM',elcom, fail_flag2)
ENDIF 
! WARNING: what to do with fail_flag1 and fail_flag2 ???
com_cot=elcom
END FUNCTION com_cot
! =======================================================
! KEP_COM
! keplerian elelements from cometary
! =======================================================
TYPE(orbit_elem) FUNCTION kep_com(el,fail_flag,del)
TYPE(orbit_elem),INTENT(IN):: el ! must be cometary
INTEGER,INTENT(OUT):: fail_flag
! optional jacobian matrix. in output: derivatives d(el)/d(eqn)
DOUBLE PRECISION, DIMENSION(6,6),INTENT(OUT),OPTIONAL :: del
! end interface
DOUBLE PRECISION ecc, eps, gm, enne, a, dldn
DOUBLE PRECISION,DIMENSION(6) :: eq
!
IF(el%coo.ne.'COM')THEN
   WRITE(*,*)' kep_com: wrong input coordinates',el
   STOP
ENDIF
fail_flag=0
gm=centerbody_mass(el%center)
eps=100*epsilon(1.d0) ! negligible number
eq=el%coord
kep_com=el  !default
ecc=eq(2)
if(ecc.ge.1.d0-eps)THEN
! this conversion is impossible
   fail_flag=5
! in output the elements are left as 'COM', del= identity
   IF(PRESENT(del)) CALL eye(6,del)
   RETURN
ENDIF
a=eq(1)/(1.d0-ecc)
kep_com%coord(1)=a
!  test on eccentricity                                                 
if(ecc.lt.eps) fail_flag= 1 ! but computation is completed
!   test on tangent of half inclination                                 
if(eq(3).lt.eps) fail_flag=fail_flag+2 ! but computation is completed
! time of pericenter from mean anomaly
enne=sqrt(gm/a**3)
kep_com%coord(6)=(el%t-eq(6))*enne 
kep_com%coo='KEP'
IF(PRESENT(del).and.fail_flag.eq.0)THEN
   del=0.d0 ! jacobian is not usable in the nearly singular case
   del(1,1)=  1.d0/(1.d0-ecc)
   del(1,2)=  +eq(1)/(1.d0-ecc)**2 
   del(2,2)=  1.d0
   del(3,3)=  1.d0
   del(4,4)=  1.d0
   del(5,5)=  1.d0
   dldn=      el%t-eq(6)
   del(6,1)=  -1.5d0*enne/eq(1)*dldn
   del(6,2)=  -1.5d0*enne/(1.d0-ecc)*dldn
   del(6,6)=  -enne
ENDIF

RETURN
END FUNCTION kep_com

! ====================================================================
! END ELEMENTS INTERNAL CONVERSIONS
! ======================================================================
! ATT submodule
! ====================================================================
! =======================================================
! ATT_CAR
! attributable elements from cartesian position and velocity, 
! planetocentric  
! attributable is:
!     %coord(1) = alpha
!     %coord(2) = delta
!     %coord(3) = d(alpha)/dt
!     %coord(4) = d(delta)/dt
!     %coord(5) = range
!     %coord(6) = range rate
! WARNING: att in equatorial, car in ecliptic coordinates 
! if the center is different from the Earth
! =======================================================
TYPE(orbit_elem) FUNCTION att_car(el,obscode,target_center,del)
  TYPE(orbit_elem),INTENT(IN):: el
  INTEGER, INTENT(IN) :: obscode
  INTEGER, INTENT(IN) ::target_center
! optional jacobian matrix. in output: derivatives d(el)/d(eqn)
  DOUBLE PRECISION, DIMENSION(6,6),INTENT(OUT),OPTIONAL :: del 
! end interface
  DOUBLE PRECISION, DIMENSION(6) :: att ! attributable
  DOUBLE PRECISION, DIMENSION(3) :: dadx,dddx,d,dv !first derivatives, 
  DOUBLE PRECISION, DIMENSION(3,3) :: ddadx,ddddx ! second derivatives
  DOUBLE PRECISION :: dz,dis0,vsize,prscal,tobs,princ ! distances, functions
  DOUBLE PRECISION, DIMENSION(6) :: xtop,xobs
! aux. var for computation of second derivatives
  DOUBLE PRECISION den,x2,y2,z2,x4,y4 
  DOUBLE PRECISION x2y2,x2z2,y2z2, ddd(3)
! =======================================================               
  IF(el%coo.ne.'CAR')THEN
     WRITE(*,*)' att_car: wrong source coordinates ', el%coo
     STOP
  ENDIF
  IF(target_center.ne.0.and.target_center.ne.3)THEN
     WRITE(*,*)'att_car: Center different from Earth or Sun'
     STOP
  END IF
! compute observer position, taking into account light time
  CALL iter_obs(el%t,el%coord,obscode,xtop,tobs,dis0,target_center)
  d=xtop(1:3)
  dv=xtop(4:6)
! Computation of observation: right ascension (radians)                 
  dz=d(1)**2+d(2)**2 
  if (dz.le.100*epsilon(1.d0)) then
! remove singularity at poles 
     att(1)=0.d0 
  else 
     att(1)=atan2(d(2),d(1)) 
  endif
! Computation of observation: declination (radians)                     
  att(2)=asin(d(3)/dis0) 
! Computation of first derivatives of $\alpha$ and $\delta$ w.r. to posi
! (if required): we derived eq. (2.20)                                  
  dadx(1)=-d(2)/dz 
  dadx(2)=d(1)/dz 
  dadx(3)=0.d0 
  dddx(1)=-d(3)*(d(1)/(sqrt(dz)*dis0**2)) 
  dddx(2)=-d(3)*(d(2)/(sqrt(dz)*dis0**2)) 
  dddx(3)=sqrt(dz)/dis0**2 
! Apparent motion:                                                      
  att(3)=prscal(dadx,dv) 
  att(4)=prscal(dddx,dv)
! range, range rate
  att(5)=dis0
  att(6)=prscal(d,dv)/dis0
! assign attributable type elements
  att_car=el !defaults
  att_car%center=3
  att_car%coord=att
  att_car%coo='ATT'
! aberration correction NO!!!
  att_car%obscode=obscode
! h magnitude is already the absolute one
! derivatives, if required
  IF(PRESENT(del))THEN
     IF(target_center.eq.0)THEN
        del=0.d0
        del(1,1:3)= MATMUL(roteqec,dadx) ! = MATMUL(dadx,roteceq)
        del(2,1:3)=MATMUL(roteqec,dddx)
        del(3,4:6)=del(1,1:3) ! commutation of d/dt
        del(4,4:6)=del(2,1:3) ! commutation of d/dt
        del(5,1:3)=MATMUL(roteqec,d*(1.d0/dis0))
        del(6,1:3)=MATMUL(roteqec,d*(-att(6)/dis0**2)+dv*(1.d0/dis0))
        del(6,4:6)=del(5,1:3) ! commutation of d/dt
! ===================================================================== 
! partials of adot,ddot with respect to positions require the           
! second derivatives of alpha, delta with respect to the                
! equatorial reference system                                           
! ===================================================================== 
! Computation of second derivatives of $\alpha$ w.r. to positions       
        ddadx(1,1)=2.d0*d(1)*d(2)/dz**2 
        ddadx(1,2)=(d(2)**2-d(1)**2)/dz**2 
        ddadx(2,1)=ddadx(1,2) 
        ddadx(2,2)=-ddadx(1,1) 
        ddadx(3,1)=0.d0 
        ddadx(3,2)=0.d0 
        ddadx(3,3)=0.d0 
        ddadx(2,3)=0.d0 
        ddadx(1,3)=0.d0 
! chain rule for derivatives of adot                                    
        ddd=MATMUL(ddadx,dv) 
! rotation to the equatorial reference system   
        del(3,1:3)=MATMUL(roteqec,ddd)
! =======================================================               
! Computation of second derivatives of $\delta$ w.r. to positions       
        den=1.d0/(dis0**4*dz*sqrt(dz)) 
        x2=d(1)**2
        y2=d(2)**2 
        z2=d(3)**2
        x4=x2*x2 
        y4=y2*y2                                    
        x2y2=x2*y2 
        x2z2=x2*z2 
        y2z2=y2*z2 
        ddddx(1,1)=d(3)*(2.d0*x4+x2y2-y2z2-y4)*den 
        ddddx(2,2)=d(3)*(2.d0*y4+x2y2-x2z2-x4)*den 
        ddddx(1,2)=d(1)*d(2)*d(3)*(z2+3.d0*x2+3.d0*y2)*den 
        ddddx(2,1)=ddddx(1,2) 
        ddddx(3,3)=-2.d0*d(3)*dz**2*den 
        ddddx(1,3)=d(1)*dz*(z2-x2-y2)*den 
        ddddx(3,1)=ddddx(1,3) 
        ddddx(2,3)=d(2)*dz*(z2-x2-y2)*den 
        ddddx(3,2)=ddddx(2,3) 
! chain rule for derivatives of ddot                                    
        CALL prodmv(ddd,ddddx,dv)
! rotation to the equatorial reference system 
        CALL prodmv(del(4,1:3),roteqec,ddd)
! the same for the earth, but without rotation
     ELSEIF(target_center.eq.3)THEN
        del=0.d0
        del(1,1:3)=dadx
        del(2,1:3)=dddx
        del(3,4:6)=del(1,1:3)
        del(4,4:6)=del(2,1:3)
        del(5,1:3)=d*(1.d0/dis0)
        del(6,1:3)=d*(-att(6)/dis0**2)+dv*(1.d0/dis0)
        del(6,4:6)=del(5,1:3) ! commutation of d/dt
        ddadx(1,1)=2.d0*d(1)*d(2)/dz**2 
        ddadx(1,2)=(d(2)**2-d(1)**2)/dz**2 
        ddadx(2,1)=ddadx(1,2) 
        ddadx(2,2)=-ddadx(1,1) 
        ddadx(3,1)=0.d0 
        ddadx(3,2)=0.d0 
        ddadx(3,3)=0.d0 
        ddadx(2,3)=0.d0 
        ddadx(1,3)=0.d0 
        CALL prodmv(ddd,ddadx,dv) 
        del(3,1:3)=ddd
        den=1.d0/(dis0**4*dz*sqrt(dz)) 
        x2=d(1)**2
        y2=d(2)**2 
        z2=d(3)**2
        x4=x2*x2 
        y4=y2*y2                                    
        x2y2=x2*y2 
        x2z2=x2*z2 
        y2z2=y2*z2 
        ddddx(1,1)=d(3)*(2.d0*x4+x2y2-y2z2-y4)*den 
        ddddx(2,2)=d(3)*(2.d0*y4+x2y2-x2z2-x4)*den 
        ddddx(1,2)=d(1)*d(2)*d(3)*(z2+3.d0*x2+3.d0*y2)*den 
        ddddx(2,1)=ddddx(1,2) 
        ddddx(3,3)=-2.d0*d(3)*dz**2*den 
        ddddx(1,3)=d(1)*dz*(z2-x2-y2)*den 
        ddddx(3,1)=ddddx(1,3) 
        ddddx(2,3)=d(2)*dz*(z2-x2-y2)*den 
        ddddx(3,2)=ddddx(2,3) 
        CALL prodmv(ddd,ddddx,dv)
        del(4,1:3)=ddd
     ELSE
        WRITE(*,*)'att_car: Center different from Earth or Sun'
        STOP
     END IF
  ENDIF
END FUNCTION att_car

! =======================================================
! ITER-OBS
! =======================================================
SUBROUTINE iter_obs(t,x,obscode,xtop,tobs,r,target_center)
USE reference_systems, ONLY: observer_position
! x is heliocentric ecliptic (tc=0), geocentric equatorial (tc=3)
DOUBLE PRECISION, INTENT(IN) :: t, x(6)
INTEGER, INTENT(IN) :: obscode, target_center
DOUBLE PRECISION, DIMENSION(6), INTENT(OUT):: xtop ! topocentric
! pos/vel with respect to observer obscode at appropriate time
! also ecliptic heliocentric position of observer
DOUBLE PRECISION, INTENT(OUT) :: tobs, r !obs time, range
DOUBLE PRECISION, DIMENSION(6) :: xpla
DOUBLE PRECISION vsize,tc, xobs(6)
INTEGER, PARAMETER :: jmax=5
DOUBLE PRECISION, PARAMETER :: tcont=1.d-8
INTEGER j
! change center: get ecliptic coordinates of the observer
tc=t
DO j=1,jmax
   IF(target_center.eq.0)THEN 
      CALL observer_position(tc,xobs(1:3),xobs(4:6),OBSCODE=obscode)
      CALL earcar(tc,xpla,1) ! xpla=heliocentric ecliptic of geocenter
      xobs=xpla+xobs ! now xobs= heliocentric ecliptic of topos
      xtop=x-xobs  ! xtop=topocentric ecliptic of asteroid 
      CALL prodmv(xtop(1:3),roteceq,xtop(1:3)) ! xtop=topocentric equatorial
      CALL prodmv(xtop(4:6),roteceq,xtop(4:6))
      IF(rhs.gt.1) STOP '**** noniter_obs: target_center=0 and rhs>1 not allowed***' 
   ELSEIF(target_center.eq.3)THEN
      CALL observer_position(tc,xobs(1:3),xobs(4:6),OBSCODE=obscode)
! rotation to equatorial is handled by observer_position if rhs=2
      IF(rhs.ne.2) STOP '**** noniter_obs: target_center=3 and rhs=1 not allowed***'
      xtop=x-xobs ! topocentric equatorial
   ELSE
      STOP
   ENDIF
! ecliptic topocentric vector of the asteroid
   r=vsize(xtop)
   tobs=t +r/vlight !aberration correction
   IF(abs(tc-tobs).lt.tcont)RETURN
   tc=tobs
ENDDO
WRITE(*,*) 'iter_obs: slow convergence ', j-1, tobs, t, r
END SUBROUTINE iter_obs

! =======================================================
! NONITER_OBS
! =======================================================
SUBROUTINE noniter_obs(tobs,obscode,xobs, target_center)
  USE reference_systems, ONLY: observer_position
  DOUBLE PRECISION, INTENT(IN) :: tobs
  INTEGER, INTENT(IN) :: obscode, target_center
  DOUBLE PRECISION, DIMENSION(6), INTENT(OUT):: xobs ! topocentric ecliptic
! pos/vel with respect to observer obscode at appropriate time
  DOUBLE PRECISION, DIMENSION(6) :: xpla
! change center: get ecliptic coordinates of the observer
  IF(target_center.eq.0)THEN
     CALL observer_position(tobs,xobs(1:3),xobs(4:6),OBSCODE=obscode)
     CALL earcar(tobs,xpla,1) ! xpla=heliocentric ecliptic of geocenter
     xobs=xpla+xobs ! now xobs= heliocentric ecliptic of topos
     IF(rhs.gt.1) STOP '**** noniter_obs: target_center=0 and rhs>1 not allowed***' 
  ELSEIF(target_center.eq.3)THEN
     CALL observer_position(tobs,xobs(1:3),xobs(4:6),OBSCODE=obscode)
! rotation to equatorial is handled by observer_position if rhs=2
     IF(rhs.ne.2) STOP '**** noniter_obs: target_center=3 and rhs=1 not allowed***' 
  ELSE
     STOP
  ENDIF
END SUBROUTINE noniter_obs

! =======================================================
! CAR_ATT
! cartesian position and velocity, planetocentric, from 
! attributable elements 
! WARNING:  car in ecliptic coordinates, att in equatorial
! =======================================================
TYPE(orbit_elem) FUNCTION car_att(el,target_center,del)
  TYPE(orbit_elem),INTENT(IN) :: el
! New variable to know the center
  INTEGER, INTENT(in) ::target_center
! optional jacobian matrix. in output: derivatives d(el)/d(eqn)
  DOUBLE PRECISION, DIMENSION(6,6),INTENT(OUT),OPTIONAL :: del 
! end interface
  DOUBLE PRECISION, DIMENSION(3) :: rhat,ralphat,rdelhat,ralpalp, &
     & rdeldel,ralpdel
  DOUBLE PRECISION, DIMENSION(6) :: att,x
  DOUBLE PRECISION cosa,cosd,sina,sind,r,rdot,tobs
  DOUBLE PRECISION, DIMENSION(6) :: xobs
!
  IF(el%coo.ne.'ATT')THEN
     WRITE(*,*)' car_att: wrong source coordinates ', el%coo
     STOP
  ENDIF
  IF(target_center.ne.3.and.target_center.ne.0)THEN
     WRITE(*,*)'car_att: The center is different from Earth and Sun'
     STOP
  ENDIF
  att=el%coord
! UNIT vectors of the r alpha delta ref. system
  cosa=cos(att(1))
  cosd=cos(att(2))
  sina=sin(att(1))
  sind=sin(att(2))
  r=att(5)
  rhat(1)=cosa*cosd
  rhat(2)=sina*cosd
  rhat(3)=sind
  ralphat(1)=-sina*cosd ! note: contains cosd and is not an unit vector
  ralphat(2)= cosa*cosd ! that is d(rhat)/d(alpha)=ralphat
  ralphat(3)=0.d0
  rdelhat(1)=-cosa*sind
  rdelhat(2)=-sina*sind
  rdelhat(3)=cosd
! second derivatives
  ralpalp(1)=-cosa*cosd
  ralpalp(2)=-sina*cosd
  ralpalp(3)=0.d0
  ralpdel(1)=sina*sind
  ralpdel(2)=-cosa*sind
  ralpdel(3)=0.d0
  rdeldel(1)=-cosa*cosd
  rdeldel(2)=-sina*cosd
  rdeldel(3)=-sind
  IF(target_center.eq.0)THEN
! planetocentric position and velocity, equatorial, converted to ecliptic
     CALL prodmv(x(1:3),roteqec,att(5)*rhat)
     CALL prodmv(x(4:6),roteqec,att(6)*rhat+r*(att(3)*ralphat+att(4)*rdelhat))
  ELSEIF(target_center.eq.3)THEN
     x(1:3)=att(5)*rhat
     x(4:6)=att(6)*rhat+r*(att(3)*ralphat+att(4)*rdelhat)
  ELSE
     STOP
  END IF
! assign output
  car_att=el
! aberration correction
  tobs=el%t+r/vlight
  CALL noniter_obs(tobs,el%obscode,xobs,target_center)
! change center: heliocentric, ecliptic/ geocentric, equatorial pos. at epoch
  car_att%center=target_center
  car_att%coord=x+xobs ! 
  car_att%coo='CAR'
! Obscode does not change for the Earth case
  car_att%obscode=500
! h magnitude is already absolute
  IF(PRESENT(del))THEN
! for the earth case
     IF(target_center.eq.3)THEN
        del=0.d0
        del(1:3,1)=ralphat*r
        del(1:3,2)=rdelhat*r
        del(1:3,5)=rhat
        del(4:6,1)=att(6)*ralphat+att(5)*(att(3)*ralpalp+att(4)*ralpdel)
        del(4:6,2)=att(6)*rdelhat+att(5)*(att(3)*ralpdel+att(4)*rdeldel)
        del(4:6,3)=del(1:3,1) ! commutation of d/dt
        del(4:6,4)=del(1:3,2) ! commutation of d/dt
        del(4:6,5)=att(3)*ralphat+att(4)*rdelhat
        del(4:6,6)=del(1:3,5) ! commutation of d/dt
     ELSE
        del=0.d0
        CALL prodmv(del(1:3,1),roteqec,ralphat*r)
        CALL prodmv(del(1:3,2),roteqec,rdelhat*r)
        CALL prodmv(del(1:3,5),roteqec,rhat)
        CALL prodmv(del(4:6,1),roteqec,att(6)*ralphat+att(5)*(att(3)*ralpalp+att(4)*ralpdel))
        CALL prodmv(del(4:6,2),roteqec,att(6)*rdelhat+att(5)*(att(3)*ralpdel+att(4)*rdeldel))
        del(4:6,3)=del(1:3,1) ! commutation of d/dt
        del(4:6,4)=del(1:3,2) ! commutation of d/dt
        CALL prodmv(del(4:6,5),roteqec,att(3)*ralphat+att(4)*rdelhat)
        del(4:6,6)=del(1:3,5) ! commutation of d/dt
     END IF
  ENDIF
END FUNCTION car_att
!=============================
! ATTELEMENTS 
! ATT elements from attributable and r,rdot
!=============================
SUBROUTINE attelements(att,r,rdot,elatt,unc)
  USE station_coordinates, ONLY: statcode
  USE attributable
  USE reference_systems, ONLY: observer_position
  TYPE(attrib), INTENT(IN):: att ! 4-dim attributable
  DOUBLE PRECISION, INTENT(IN) :: r,rdot ! range, range rate
  TYPE(orbit_elem), INTENT(OUT):: elatt ! elements of type ATT
  TYPE(orb_uncert), INTENT(OUT), OPTIONAL :: unc ! covar/normal matrices
! end interface
  TYPE(orbit_elem) xcar
  INTEGER fail_flag,indp
  DOUBLE PRECISION rsun,phase,cosph, xx(3),xpla(6),vsize,prscal,appmag,ws(4)
  DOUBLE PRECISION pos(3),vel(3),xx1(3),xpla1(6)
!=============================
  elatt=undefined_orbit_elem
  elatt%center=3 ! ???????????????????????????
  elatt%coord(1:4)=att%angles
  elatt%coord(5)=r
  elatt%coord(6)=rdot
  elatt%coo='ATT'
  elatt%t=att%tdtobs-r/vlight
  CALL statcode(att%obscod, elatt%obscode)
! cartesian coordinates needed to get phase right
  CALL coo_cha(elatt,'CAR',xcar,fail_flag)
  xx=xcar%coord(1:3)
  CALL earcar(att%tdtobs,xpla,1)

  IF(rhs.eq.1)THEN
! xx1 is heliocentric ecliptic position of asteroid?????
     xx1=xx
  ELSEIF(rhs.eq.2)THEN
! xx1 is heliocentric ecliptic position of debris??????????
     CALL prodmv(xx1,roteqec,xx)
     xx1=xx1+xpla(1:3)
  ELSE
     WRITE(*,*)'attelements: rhs=', rhs
     STOP
  END IF
  rsun=vsize(xx1)
! xpla1 is the heliocentric/ecliptic observer position
  IF(rhs.eq.1)THEN
     xpla1=xpla
  ELSEIF(rhs.eq.2)THEN
     CALL observer_position(att%tdtobs,pos,vel,OBSCODE=elatt%obscode)
     pos=MATMUL(roteqec,pos)
     vel=MATMUL(roteqec,vel)
     xpla1(1:3)=pos+xpla(1:3)
     xpla1(4:6)=vel+xpla(4:6)
  ELSE
     WRITE(*,*)'attelements: rhs=', rhs
     STOP
  END IF
  cosph=prscal(xx1,xx1-xpla1(1:3))/(rsun*r)
  IF(cosph.gt.1.d0)cosph=1.d0
  IF(cosph.lt.-1.d0)cosph=-1.d0
  phase=acos(cosph)
  IF(att%apm.gt.0.d0)THEN
     elatt%mag_set=.true.
     elatt%h_mag=att%apm-appmag(0.d0,elatt%g_mag,rsun,r,phase)
  ELSE
     elatt%mag_set=.false.
  ENDIF
  IF(PRESENT(unc))THEN
     CALL undefined_orb_uncert(6,unc)
     unc%g=0.d0
     unc%c=0.d0
     unc%g(1:4,1:4)=att%g
     unc%succ=.true.
     CALL tchinv(att%g,4,unc%c(1:4,1:4),ws,indp)
  ENDIF
END SUBROUTINE attelements

! ================================================
! ATT_PRELIM
! attributable preliminary orbit
! ================================================
SUBROUTINE att_prelim2(name0,attr,elk,uncatt,fail)
  USE attributable
  CHARACTER*(*), INTENT(IN) :: name0 ! asteroid designation
  TYPE(attrib), INTENT(IN) :: attr  ! attributable as computed 
  TYPE(orbit_elem), INTENT(OUT) :: elk ! new elements
  TYPE(orb_uncert), INTENT(OUT) :: uncatt ! its covariance
  LOGICAL, INTENT(OUT) :: fail
! END INTERFACE
  INTEGER iunmat,nroots,nrootsd
  DOUBLE PRECISION ::  c(0:5)! polynomial coefficients
  DOUBLE PRECISION, DIMENSION(3) :: roots ! real roots
  DOUBLE PRECISION, DIMENSION(8) :: rootsd ! roots of the derivative
  DOUBLE PRECISION :: a_orb,E_bound ! bound for the energy
! matrix of the r eps theta ref. system, topocentric equatorial
  DOUBLE PRECISION :: refs(3,3)
 ! observer pos/vel, asteroid pos/vel, equatorial heliocentric
  DOUBLE PRECISION :: xo(3), vo(3)
  DOUBLE PRECISION rr, rdot,ecc

  LOGICAL bizarre
  a_orb = 100.d0 ! comet boundary
!   a_orb = 2.5d0 ! main belt
  fail=.true.
  iunmat=0
  CALL admis_reg(name0,iunmat,attr%tdtobs,attr%angles,attr%obscod,  &
   &  nroots,roots,c,refs,xo,vo,a_orb,E_bound,nrootsd,rootsd)
! handle rare cases
  IF(nroots.eq.3)THEN
     write(iun_log,*)'att_prelim: ',name0,' two c.c.',roots(1:nroots)
     rr=(roots(3)+roots(2))/2.d0  ! privilege given to farther component
  ELSEIF(nroots.eq.2)THEN
     WRITE(iun_log,*)'att_prelim: ',name0, ' rare, nroots= ',nroots,roots(1:nroots)
     rr=roots(1)*0.8d0 ! somewhat inside the admissible region
  ELSEIF(nroots.eq.1)THEN
     rr=roots(1)*0.8d0 ! somewhat inside the admissible region
  ENDIF
  rdot=-c(1)/2.d0
  CALL attelements(attr,rr,rdot,elk,uncatt)
  
  fail=bizarre(elk,ecc)
  IF(verb_prelim.gt.10)THEN
     WRITE(iun_log,100)name0,rr,ecc
100  FORMAT(A,1X,'rho=',f8.3,1X,'ecc=',f8.4)
  ENDIF
END SUBROUTINE att_prelim2


! ====================================================================
! END ATT submodule
! ====================================================================
! ====================================================================
! OPIK submodule
! ====================================================================
! TPCAR_MTPCAR from perihelion position and velocity, obtain position on the 
! asymptote with rotation by -gamma/2 and rescaling by b/d of position,
! d/b of velocity
! ====================================================================
TYPE(orbit_elem) FUNCTION tpcar_mtpcar(el,fail_flag,bd,del)
  TYPE(orbit_elem), INTENT(IN) :: el ! must be cartesian
  INTEGER, INTENT(OUT) :: fail_flag 
  DOUBLE PRECISION, INTENT(OUT), OPTIONAL :: bd ! b/d rescaling factor
  DOUBLE PRECISION, DIMENSION(6,6), OPTIONAL :: del ! jacobian matrix
! =============END INTERFACE==========================================
  DOUBLE PRECISION gm, v, d, eps, sig0, vsize, prscal
  DOUBLE PRECISION bsd, bsd2, sinhalfgam,coshalfgam
  DOUBLE PRECISION dbsddd, dbsddv, dgdd, dgdv ! derivatives of b/d, sinhalfgamma w.r. to d,v
  DOUBLE PRECISION, DIMENSION(3) :: x,y, ang, xtp,ytp, xmtp,xhat,yhat, xp,yp, xphat,yphat
  DOUBLE PRECISION, DIMENSION(3,3) :: rottp,drotpdg2 ! R, dr/d(sin(gam/2)
  DOUBLE PRECISION, DIMENSION(6,6) :: dlbd, dxpdx ! derivatives of rescaling,
                          ! of rotation by -halfgamma 
  DOUBLE PRECISION, DIMENSION(3,3,3) :: dtaudJ ! d(rotation)/d(ang)
  DOUBLE PRECISION, DIMENSION(3,3,3) :: ericci ! Ricci tensor
  DOUBLE PRECISION, DIMENSION(3,3)   :: dbJdX
  DOUBLE PRECISION, DIMENSION(3,3)   :: dbJdY
  DOUBLE PRECISION                   :: normJ
  INTEGER i,j,k,h ! =1,3
  IF(el%coo.ne.'CAR')THEN
     WRITE(*,*)' tpcar_mtpcar: wrong coords; center ',el%coo,' ', el%center
     STOP
  ENDIF
  fail_flag=0
  gm=centerbody_mass(el%center)
  tpcar_mtpcar=el ! copy defaults
  x=el%coord(1:3) ! position vector 
  y=el%coord(4:6) ! velocity vector
!  radius and velocity                                          
  v=vsize(y)      ! velocity
  xmtp=x-prscal(x,y)*y/v**2
  d=vsize(xmtp)      ! distance from center
  eps=2.d-6
! check on orthogonality
  sig0=prscal(x,y)/(v*d)
  IF(abs(sig0).gt.eps)THEN
     WRITE(*,199) sig0,eps
 199 FORMAT(' tpcar_mtpcar: pos, vel not orthogonal ',1P,2d10.3)
  ENDIF
  call prvec(x,y,ang)           !  angular momentum 
! conformal change in the TP
  bsd2=(d*v**2)/(d*v**2-2*gm)
  IF(bsd2.gt.0.d0)THEN
     bsd=sqrt(bsd2)
  ELSE
     fail_flag=6 ! means the orbit is not hyperbolic
     RETURN
  ENDIF
  sinhalfgam=gm/(v**2*d-gm)
  coshalfgam=sqrt(1.d0-sinhalfgam**2)
! halfgamma=asin(sinhalfgam)
  normJ=vsize(ang)
  ang=ang/normJ
! coordinates rotated by -halfgamma around ang
  CALL rotmhalfgam(sinhalfgam,ang,rottp,drotpdg2,dtaudJ)
  xp=MATMUL(rottp,x)
  yp=MATMUL(rottp,y)
!  xp=x
!  yp=y
! rescaling to asymptote
  xtp=xp*bsd
  ytp=yp*(1.d0/bsd)  
  tpcar_mtpcar%coord(1:3)=xtp
  tpcar_mtpcar%coord(4:6)=ytp
!  tpcar_mtpcar%coord(1:3)=xp
!  tpcar_mtpcar%coord(4:6)=yp
  IF(PRESENT(bd))THEN
    bd=bsd
  ENDIF
  tpcar_mtpcar%coo='TPC'
! partial derivatives
  IF(PRESENT(del))THEN
! derivatives d(xp,yp)/d(x,y) 
! neglecting derivatives with respect to ang
     xhat=x*(1.d0/d)
     yhat=y*(1.d0/v)
     dgdd=-v**2/gm*sinhalfgam**2
     dgdv=-2.d0*v*d/gm*sinhalfgam**2
     dxpdx=0.d0
     dxpdx(1:3,1:3)=rottp
     dxpdx(4:6,4:6)=rottp
     DO i=1,3
       DO j=1,3 
         DO k=1,3
             dxpdx(i,j)=dxpdx(i,j)+drotpdg2(i,k)*x(k)*dgdd*xhat(j)
             dxpdx(i+3,j+3)=dxpdx(i+3,j+3)+drotpdg2(i,k)*y(k)*dgdv*yhat(j)
             dxpdx(i,j+3)=dxpdx(i,j+3)+drotpdg2(i,k)*x(k)*dgdv*yhat(j)
             dxpdx(i+3,j)=dxpdx(i+3,j)+drotpdg2(i,k)*y(k)*dgdd*xhat(j)
         ENDDO
       ENDDO
     ENDDO
! derivatives of rottp w.r. to ang
     ericci(1:3,1:3,1:3)=0.d0
     ericci(1,2,3)=1.d0
     ericci(2,3,1)=1.d0
     ericci(3,1,2)=1.d0
     ericci(1,3,2)=-1.d0
     ericci(3,2,1)=-1.d0
     ericci(2,1,3)=-1.d0
     DO j=1,3
        DO h=1,3
         dbJdX(h,j)=1.d0/normJ*(DOT_PRODUCT(ericci(j,1:3,h),y)-v/d*x(j)*ang(h))
         dbJdY(h,j)=1.d0/normJ*(DOT_PRODUCT(ericci(1:3,j,h),x)-d/v*y(j)*ang(h))
      ENDDO
     ENDDO
! add term depending upon d(rottp0/d(ang)
      DO i=1,3
        DO j=1,3
          DO k=1,3
             dxpdx(i,j)=dxpdx(i,j)+DOT_PRODUCT(dtaudJ(i,k,1:3),dbJdX(1:3,j))*x(k)
             dxpdx(i,j+3)=dxpdx(i,j+3)+DOT_PRODUCT(dtaudJ(i,k,1:3),dbJdY(1:3,j))*x(k)
             dxpdx(i+3,j)=dxpdx(i+3,j)+DOT_PRODUCT(dtaudJ(i,k,1:3),dbJdX(1:3,j))*y(k)
             dxpdx(i+3,j+3)=dxpdx(i+3,j+3)+DOT_PRODUCT(dtaudJ(i,k,1:3),dbJdY(1:3,j))*y(k)
          ENDDO
        ENDDO
     ENDDO
! derivatives d(txp,typ)/d(xp,yp) 
     xphat=xp*(1.d0/d)
     yphat=yp*(1.d0/v)
     dbsddd=1.d0/(2.d0*bsd)*(-2.d0*gm*v**2)/(v**2*d-2.d0*gm)**2
     dbsddv=1.d0/(2.d0*bsd)*(-4.d0*gm*v*d)/(v**2*d-2.d0*gm)**2
     dlbd=0.d0
     DO i=1,3
       dlbd(i,i)=bsd
       dlbd(3+i,3+i)=1.d0/bsd
       DO j=1,3
         dlbd(i,j)=dlbd(i,j)+dbsddd*xphat(j)*xp(i)
         dlbd(i,j+3)=dbsddv*yphat(j)*xp(i)
         dlbd(i+3,j)=(-1.d0/bsd**2)*dbsddd*xphat(j)*yp(i)
         dlbd(i+3,j+3)=dlbd(i+3,j+3)+(-1.d0/bsd**2)*dbsddv*yphat(j)*yp(i)
       ENDDO
     ENDDO
     del=MATMUL(dlbd,dxpdx)
!     del=dlbd
!     del=dxpdx
  ENDIF
CONTAINS
SUBROUTINE rotmhalfgam(sgmez,bJ,tau,dtaudsgm,dtaudJ)
  DOUBLE PRECISION,INTENT(IN) :: sgmez !sinus(gamma/2)
  DOUBLE PRECISION,DIMENSION(3),INTENT(IN) :: bJ !ang. momentum unit vector
! rotation of -gamma/2 around bJ and its derivative w.r.t. sinus(gamma/2)
  DOUBLE PRECISION,DIMENSION(3,3),INTENT(OUT) :: tau,dtaudsgm 
! derivative of tau w.r.t. the components of bJ
! dtaudJ(*,*,k) is the derivative w.r.t. the k-th component  
  DOUBLE PRECISION,DIMENSION(3,3,3),INTENT(OUT) :: dtaudJ 
! ------------- end interface -------------------------------------------
  DOUBLE PRECISION :: normJ,HRoy
  DOUBLE PRECISION :: cI,sI,cO,sO, Incl, Omeg, gmez
  DOUBLE PRECISION :: cgmez,tgmez ! cosinus and tangent of gamma/2
! rotation matrices
  DOUBLE PRECISION, DIMENSION(3,3) :: ROmeg
  DOUBLE PRECISION, DIMENSION(3,3) :: RmOmeg
  DOUBLE PRECISION, DIMENSION(3,3) :: RIncl
  DOUBLE PRECISION, DIMENSION(3,3) :: RmIncl
  DOUBLE PRECISION, DIMENSION(3,3) :: Rmhgam
! for the derivatives
  DOUBLE PRECISION,DIMENSION(3) :: dHRoydJ
  DOUBLE PRECISION, DIMENSION(3,3) :: delta
  DOUBLE PRECISION, DIMENSION(3,3) :: dRmhgam
  DOUBLE PRECISION, DIMENSION(3,3) :: RORI,RmIRmO,tauauxr,tauauxl !auxiliary
  DOUBLE PRECISION, DIMENSION(3,3,3) :: dROmegdJ,dRmOmegdJ
  DOUBLE PRECISION, DIMENSION(3,3,3) :: dRIncldJ,dRmIncldJ
  DOUBLE PRECISION, DIMENSION(3,3,3) :: tensor1,tensor2,tensor3,tensor4
  DOUBLE PRECISION :: vsize !functions
  INTEGER :: h,k ! loop indexes
! ---------------------------------------------
  normJ=vsize(bJ)
  IF(abs(normJ-1.d0).gt.100*epsilon(1.d0))THEN
     WRITE(*,*) 'rotmhalfgam: angular momentum not normalized! normJ=',normJ
  ENDIF
  HRoy=sqrt(bJ(1)**2+bJ(2)**2)
  IF(HRoy/normJ.lt.epsilon(1.d0)*100)THEN
     cI=1.d0
     sI=0.d0
     cO=1.d0
     sO=0.d0
  ELSE
     cI = bJ(3) ! bJ(3)/normJ (but J is unitary)
     sI = HRoy  ! HRoy/normJ
     cO = -bJ(2)/HRoy ! singular for zero I
     sO = bJ(1)/HRoy ! singular for zero I
  ENDIF
  cgmez = sqrt(1.d0-sgmez**2)
  tgmez = sgmez/cgmez

  Incl=atan2(sI,cI)
  Omeg=atan2(sO,cO)
  gmez=atan2(sgmez,cgmez)

! *******************************************************
! rotation matrix tau = ROmeg*RIncl*Rmhgam*RmIncl*RmOmeg
! *******************************************************
  CALL rotmt(-Omeg,ROmeg,3)
  CALL rotmt(Omeg,RmOmeg,3)
  CALL rotmt(-Incl,RIncl,1)
  CALL rotmt(Incl,RmIncl,1)
  CALL rotmt(gmez,Rmhgam,3)
  RORI=matmul(ROmeg,RIncl)
  RmIRmO=matmul(RmIncl,RmOmeg)
  tauauxr = matmul(Rmhgam,RmIRmO)
  tauauxl = matmul(RORI,Rmhgam) ! will be used later
  tau = matmul(RORI,tauauxr)

! *********************************************
! derivative of Rmhgam w.r.t. sgmez 
! dtaudsgm = ROmeg*RIncl*dRmhgam*RmIncl*RmOmeg
! *********************************************
  dRmhgam(1:3,1:3)=0.d0! initialization
  dRmhgam(1,1)=-tgmez
  dRmhgam(1,2)=1.d0
  dRmhgam(2,1)=-1.d0
  dRmhgam(2,2)=-tgmez
  dtaudsgm = matmul(dRmhgam, RmIRmO) ! temporary step
  dtaudsgm = matmul(RORI, dtaudsgm)

! *************************************************
! derivative of Rmhgam w.r.t. the components of bJ
! *************************************************
!  IF(HRoy/normJ.lt.epsilon(1.d0)*100)THEN
     dHRoydJ(1) = sO
     dHRoydJ(2) = -cO
!  ELSE
!     dHRoydJ(1) = bJ(1)/HRoy ! singular for zero I
!     dHRoydJ(2) = bJ(2)/HRoy ! singular for zero I
!  ENDIF
  dHRoydJ(3) = 0.d0

  CALL eye(3,delta) ! Kronecker's delta

  dROmegdJ(1:3,1:3,1:3)=0.d0 ! initialization
  IF(HRoy/normJ.ge.epsilon(1.d0)*100)THEN
     dROmegdJ(1,1,1:3)= -(1.d0/HRoy)*delta(2,1:3) + (bJ(2)/HRoy**2)*dHRoydJ(1:3) ! singular for zero I
     dROmegdJ(1,2,1:3)= -(1.d0/HRoy)*delta(1,1:3) + (bJ(1)/HRoy**2)*dHRoydJ(1:3) ! singular for zero I
  ELSE
!     dROmegdJ(1,1,1:3)= -(1.d0/HRoy)*delta(2,1:3) +  (-cO/HRoy)*dHRoydJ(1:3) ! singular for zero I
!     dROmegdJ(1,2,1:3)= -(1.d0/HRoy)*delta(1,1:3) +  (sO/HRoy)*dHRoydJ(1:3) ! singular for zero I
  ENDIF
  dROmegdJ(2,1,1:3)= -dROmegdJ(1,2,1:3)
  dROmegdJ(2,2,1:3)=  dROmegdJ(1,1,1:3)
  dRmOmegdJ(1:3,1:3,1:3)= dROmegdJ(1:3,1:3,1:3)
  dRmOmegdJ(1,2,1:3)= -dROmegdJ(1,2,1:3)
  dRmOmegdJ(2,1,1:3)= -dRmOmegdJ(1,2,1:3)
  dRIncldJ(1:3,1:3,1:3)=0.d0 ! initialization
  dRIncldJ(2,2,1:3)=  delta(3,1:3)
  dRIncldJ(2,3,1:3)= -dHRoydJ(1:3)
  dRIncldJ(3,2,1:3)= -dRIncldJ(2,3,1:3)
  dRIncldJ(3,3,1:3)=  dRIncldJ(2,2,1:3)

  dRmIncldJ(1:3,1:3,1:3)= dRIncldJ(1:3,1:3,1:3)
  dRmIncldJ(2,3,1:3)= -dRIncldJ(2,3,1:3)
  dRmIncldJ(3,2,1:3)= -dRmIncldJ(2,3,1:3)

! **********************************************
! tensor1 = dROmegdJ*RIncl*Rmhgam*RmIncl*RmOmeg
! **********************************************
  DO k=1,3
     tensor1(1:3,1:3,k) = matmul(RIncl,tauauxr)
     tensor1(1:3,1:3,k) = matmul(dROmegdJ(1:3,1:3,k),tensor1(1:3,1:3,k))
  ENDDO
! **********************************************
! tensor2 = ROmeg*dRIncldJ*Rmhgam*RmIncl*RmOmeg
! **********************************************
  DO k=1,3
     tensor2(1:3,1:3,k) = matmul(dRIncldJ(1:3,1:3,k),tauauxr)
     tensor2(1:3,1:3,k) = matmul(ROmeg,tensor2(1:3,1:3,k))
  ENDDO
! **********************************************
! tensor3 = ROmeg*RIncl*Rmhgam*dRmIncldJ*RmOmeg
! **********************************************
  DO k=1,3
     tensor3(1:3,1:3,k) = matmul(tauauxl,dRmIncldJ(1:3,1:3,k))
     tensor3(1:3,1:3,k) = matmul(tensor3(1:3,1:3,k),RmOmeg)
  ENDDO
! **********************************************
! tensor4 = ROmeg*RIncl*Rmhgam*RmIncl*dRmOmegdJ
! **********************************************
  DO k=1,3
     tensor4(1:3,1:3,k) = matmul(tauauxl,RmIncl)
     tensor4(1:3,1:3,k) = matmul(tensor4(1:3,1:3,k),dRmOmegdJ(1:3,1:3,k))
  ENDDO
  
  dtaudJ = tensor1 + tensor2 + tensor3 + tensor4

END SUBROUTINE rotmhalfgam

END FUNCTION tpcar_mtpcar
! ====================================================================
! OPIK_TPCAR from perihelion position and velocity, obtain Opik elements 
! ====================================================================
! opik= U, long, lat, csi, zeta, eta
! U= unperturbed v_infty
! long, lat  = antiradians (in ecliptical coordinates)
! csi, zeta = TP coordinates with axis zeta along projection of normal to 
!                   ecliptic
! eta = U(t-t0), t0 time of pericenter passage
! WARNING: in .clo file we need also the time dependent theta, phi 
! to be computed using the Earth velocity at t0
! ====================================================================
TYPE(orbit_elem) FUNCTION opik_tpcar(el,fail_flag,del,vt3)
  TYPE(orbit_elem), INTENT(IN) :: el ! must be cartesian
  INTEGER, INTENT(OUT) :: fail_flag 
  DOUBLE PRECISION, DIMENSION(6,6), OPTIONAL :: del ! jacobian matrix
  DOUBLE PRECISION, DIMENSION(3,3), OPTIONAL :: vt3 ! rotation of TP to plane eta=0
! =============END INTERFACE==========================================
  DOUBLE PRECISION opik(6), xytp(6), vsize, b, dz, prscal, cosa, dt
  DOUBLE PRECISION v1(3), v2(3), rv(3,3),xx(3),drv(3,3,3)!,rvt(3,3), rv1(3,3)
  INTEGER i,j
  IF(el%coo.ne.'TPC')THEN
     WRITE(*,*)' opik_tpcar: wrong coordinates ', el
     STOP
  ENDIF
  fail_flag=0
  opik_tpcar=el ! set default; note that epoch remains the one of pericenter
  xytp=el%coord
  opik(1)=vsize(xytp(4:6))
  v1=xytp(4:6)*(1.d0/opik(1))
  b=vsize(xytp(1:3))
  cosa=prscal(xytp(1:3),xytp(4:6))/(b*opik(1))
  dt=-b*cosa/opik(1)
  IF(abs(cosa).gt.2.d-6)THEN
!     xytp(1:3)=xytp(1:3)-b*cosa*v1
!     b=vsize(xytp(1:3))
     WRITE(*,199)cosa, b, opik(1)
 199 FORMAT('opik_tpcar: pos not orthogonal to vel ',1p,3D10.3)
!     fail_flag=7
  ENDIF
! polar coordinates                 
  dz=xytp(4)**2+xytp(5)**2 
  if (dz.le.100*epsilon(1.d0)) then
! remove singularity at poles 
     opik(2)=0.d0 
  else 
! ecliptic longitude of velocity (radians)
     opik(2)=atan2(xytp(5),xytp(4)) 
     if (opik(2).lt.0.d0) then 
        opik(2)=opik(2)+dpig 
     endif
  endif
! ecliptic latitude of velocity (radians)                     
  opik(3)=asin(xytp(6)/opik(1))
! target plane 
  CALL drvdytp(xytp(4:6),drv,rv)
  IF(PRESENT(vt3))THEN
     vt3=rv
  ENDIF
  v1=xytp(4:6)*(1.d0/opik(1))
!  v2=(/0.d0, 0.d0, 1.d0/) 
!  CALL mtp_ref(v1,v2,rvt,rv) 
  xx=MATMUL(rv,xytp(1:3))
  opik(4)=xx(2)
  opik(5)=xx(3)
  IF(abs(xx(1)).gt.2.d-7)THEN
     WRITE(*,*)' opik_tpcar: pos, vel not orthogonal ', xx(1), xx(1)/b
  ENDIF
! last element is the time of pericenter passage 
  opik(6)=opik(1)*dt
! resulting elements
  opik_tpcar%t=el%t
  opik_tpcar%coord=opik
  IF(PRESENT(del))THEN
! partial derivatives
     del(1:3,1:3)=0.d0
     del(1,4:6)=v1
     del(2,4)=-xytp(5)/dz 
     del(2,5)=xytp(4)/dz 
     del(2,6)=0.d0 
     del(3,4)=-xytp(6)*(xytp(4)/(sqrt(dz)*opik(1)**2)) 
     del(3,5)=-xytp(6)*(xytp(5)/(sqrt(dz)*opik(1)**2)) 
     del(3,6)=sqrt(dz)/opik(1)**2
     del(4:5,4:6)=0.d0 ! neglecting derivatives of rv
     del(4:5,1:3)=rv(2:3,1:3) 
     del(6,1:3)=-(1.d0/opik(1))*xytp(4:6)
     del(6,4:6)=-(1.d0/opik(1))*xytp(1:3)+b/opik(1)*cosa*v1
! contribution from derivatives of rv
     DO i=4,5
       DO j=4,5
         del(i,j)=del(i,j)+DOT_PRODUCT(drv(i-2,1:3,j-3),xytp(1:3))
       ENDDO
     ENDDO
  ENDIF
  CONTAINS

SUBROUTINE drvdytp(ytp,drv,rv)
DOUBLE PRECISION, DIMENSION(3), INTENT(IN) :: ytp
                  ! Velocity vector in TP reference frame
DOUBLE PRECISION, DIMENSION(3,3), INTENT(OUT) :: rv
                  ! rotation matrix        
DOUBLE PRECISION, DIMENSION(3,3,3), INTENT(OUT) :: drv 
                  ! Derivatives of rotation matrix rv
DOUBLE PRECISION :: U,sqrx2y2,trhalf,d1,d2, x,y,z
                  ! Auxiliary quantities
! ==================================================================
 U=sqrt(ytp(1)**2+ytp(2)**2+ytp(3)**2)          ! Module of velocity
 x=ytp(1)
 y=ytp(2)
 z=ytp(3)
! Matrix rv
 rv(1,1)=x/U  
 rv(1,2)=y/U  
 rv(1,3)=z/U
 rv(3,1)=-x*z/(U*sqrt(x**2+y**2))  
 rv(3,2)=-y*z/(U*sqrt(x**2+y**2))
 rv(3,3)=(x**2+y**2)/(U*sqrt(x**2+y**2))
 rv(2,1)=y/sqrt(x**2+y**2) 
 rv(2,2)=-x/sqrt(x**2+y**2)
 rv(2,3)=0.d0

! ==================================================================

sqrx2y2=sqrt(ytp(1)**2+ytp(2)**2)
trhalf=sqrx2y2*(ytp(1)**2+ytp(2)**2)
d1=U**3*trhalf
d2=U**3*sqrx2y2
! ===================================================
! Derivatives of first row of matrix R_v
drv(1,1,1)=(ytp(2)**2+ytp(3)**2)/U**3
drv(1,1,2)=-ytp(1)*ytp(2)/U**3
drv(1,1,3)=-ytp(1)*ytp(3)/U**3

drv(1,2,1)=-ytp(1)*ytp(2)/U**3
drv(1,2,2)=(ytp(1)**2+ytp(3)**2)/U**3
drv(1,2,3)=-ytp(2)*ytp(3)/U**3

drv(1,3,1)=-ytp(1)*ytp(3)/U**3
drv(1,3,2)=-ytp(2)*ytp(3)/U**3
drv(1,3,3)=(ytp(1)**2+ytp(2)**2)/U**3
! ===================================================
! Derivatives of third row of matrix R_v
drv(3,1,1)=ytp(3)*((ytp(1)**2-U**2)*(ytp(1)**2+ytp(2)**2)+U**2*ytp(1)**2)/d1
drv(3,1,2)=ytp(1)*ytp(2)*ytp(3)*(ytp(1)**2+ytp(2)**2+U**2)/d1
drv(3,1,3)=ytp(1)*(ytp(3)**2-U**2)/d2

drv(3,2,1)=ytp(1)*ytp(2)*ytp(3)*(ytp(1)**2+ytp(2)**2+U**2)/d1
drv(3,2,2)=ytp(3)*((ytp(2)**2-U**2)*(ytp(1)**2+ytp(2)**2)+U**2*ytp(2)**2)/d1
drv(3,2,3)=ytp(2)*(ytp(3)**2-U**2)/d2

drv(3,3,1)=ytp(1)*ytp(3)**2/d2
drv(3,3,2)=ytp(2)*ytp(3)**2/d2
drv(3,3,3)=-ytp(3)*sqrx2y2/U**3
! ====================================================
! Derivatives of second row of matrix R_v
drv(2,1,1)=-ytp(1)*ytp(2)/trhalf
drv(2,1,2)=ytp(1)**2/trhalf
drv(2,1,3)=0.d0

drv(2,2,1)=-ytp(2)**2/trhalf
drv(2,2,2)=ytp(1)*ytp(2)/trhalf
drv(2,2,3)=0.d0

drv(2,3,1)=0.d0
drv(2,3,2)=0.d0
drv(2,3,3)=0.d0
! =====================================================
END SUBROUTINE drvdytp


END FUNCTION opik_tpcar
! ====================================================================
! END OPIK submodule
! ====================================================================
! ====================================================================
! I/O submodule
! ====================================================================
SUBROUTINE read_elems_old(el,name,eof,file,unit,covar,dp,nd)
  USE io_elems, ONLY: obscod
  USE station_coordinates, ONLY: statcode
  TYPE(orbit_elem), INTENT(OUT) :: el
  CHARACTER*(*),INTENT(OUT) :: name
  LOGICAL, INTENT(OUT) :: eof
  CHARACTER*(*), INTENT(IN) :: file ! input file
  INTEGER, INTENT(IN), OPTIONAL :: unit ! input unit, only if file already opened
  DOUBLE PRECISION,INTENT(OUT),DIMENSION(ndyx),OPTIONAL :: dp ! non grav parameters
  INTEGER,INTENT(IN),OPTIONAL :: nd ! number of solve for parameters
  TYPE(orb_uncert), INTENT(OUT), OPTIONAL :: covar 
! end interface
  DOUBLE PRECISION elem(6),t0,h,g,cove(ndimx,ndimx),nore(ndimx,ndimx),mass !,cnv(6)
  CHARACTER*3 eltype
  CHARACTER*10 rsys,epoch 
  LOGICAL defcov,defnor,defcn
  INTEGER kr, iobs 
  INTEGER nd_cur ! number of solve for parameters (either nd if present, or 6)

  IF(PRESENT(nd))THEN
     nd_cur=nd
  ELSE
     nd_cur=6
  ENDIF

! if present dp, then nd is required
  IF(PRESENT(dp).and..not.PRESENT(nd))THEN
     WRITE(*,*)'read_elems_old: nd not present'
     STOP
  ENDIF

  IF(.not.PRESENT(unit))THEN
     CALL oporbf(file,0)
  ENDIF
  el=undefined_orbit_elem
  IF(rhs.EQ.2)THEN
     el%center=3
  END IF
  IF(PRESENT(dp))THEN
     dp=0.d0
     CALL rd_orb(name,elem,eltype,t0,cove,defcov,nore,defnor,     &
          &                 h,g,mass,rsys,epoch,kr,eof,dp,nd)
  ELSE
     CALL rd_orb(name,elem,eltype,t0,cove,defcov,nore,defnor,     &
          &                 h,g,mass,rsys,epoch,kr,eof)
  ENDIF
  IF(eof)THEN
     IF(.not.PRESENT(unit)) CALL clorbf
     RETURN
  ENDIF
! unit conversions already done by rdorb
  IF(eltype.eq.'ATT')THEN
! input obscod
     CALL statcode(obscod,iobs)
  END IF
! what to do with rsys,epoch????
  IF(rsys.ne.'ECLM'.or.epoch.ne.'J2000')THEN
     WRITE(*,*)' read_elems_old: other reference systems not handled ', rsys,'  ', epoch
     STOP
  ENDIF
! copy into output elements
  el%coo=eltype
  el%coord=elem
  el%t=t0
  IF(h.lt.-10.d0)THEN
     el%mag_set=.false.
  ELSE
     el%mag_set=.true.
     el%h_mag=h
     el%g_mag=g
  ENDIF
  IF(el%coo.eq.'ATT')el%obscode=iobs

  IF(PRESENT(covar))THEN
     CALL undefined_orb_uncert(nd_cur,covar)
     IF(defcov) covar%g(1:nd_cur,1:nd_cur)=cove(1:nd_cur,1:nd_cur)
     IF(defnor) covar%c(1:nd_cur,1:nd_cur)=nore(1:nd_cur,1:nd_cur)
     CALL fixcnm(defcov,defnor,defcn,cove,nore)
     IF(defcn)THEN
        covar%succ=.true.
     ELSE
        covar%succ=.false.
     ENDIF
  ENDIF

  IF(.not.PRESENT(unit)) CALL clorbf
END SUBROUTINE read_elems_old

SUBROUTINE write_elems_old(el,name,form,file,unit,covar,incfit,dp)
  USE name_rules
  USE io_elems, ONLY: obscod
  USE station_coordinates, ONLY: codestat
  USE dyn_param
  TYPE(orbit_elem), INTENT(IN) :: el
  CHARACTER*(*),INTENT(IN) :: name
  CHARACTER*(*), INTENT(IN) :: form ! can be either '1L' or 'ML'
  CHARACTER*(*), INTENT(IN), OPTIONAL :: file ! output file 
                       ! (not used if PRESENT(unit))
  INTEGER, INTENT(IN), OPTIONAL :: unit ! input unit, only if already opened
! WARNING: if file already opened, it must also have the header already written
  TYPE(orb_uncert), INTENT(IN), OPTIONAL :: covar 
  INTEGER, INTENT(IN), OPTIONAL :: incfit ! result from incomplete fit;
    ! warn the users of this!!!!
  DOUBLE PRECISION, OPTIONAL, INTENT(IN) :: dp(ndyx)
! =========== end interface================
  INTEGER :: nd
  INTEGER uniout,le
  CHARACTER*(idnamvir_len) namloc
  CHARACTER*10 rsys,epoch
  LOGICAL defcov,defnor
  DOUBLE PRECISION mass, cove(ndimx,ndimx), nore(ndimx,ndimx) 
! number of parameters
  nd=6+nls
! open output file if not opened already
  IF(.not.PRESENT(unit))THEN
     CALL filopn(uniout,file,'unknown')
! write file header
     rsys='ECLM'
     epoch='J2000'
     IF(form.eq.'1L')THEN
! what to do when el%coo=ATT ?????
        CALL codestat(el%obscode,obscod)
        CALL wro1lh(uniout,rsys,epoch,el%coo)
     ELSEIF(form.eq.'ML')THEN
        CALL wromlh(uniout,rsys,epoch)
     ENDIF
  ELSE
     uniout=unit
  ENDIF
  namloc=name
  CALL rmsp(namloc,le)
  IF(PRESENT(incfit))THEN
     IF(incfit.eq.5)THEN
        WRITE(uniout,200)namloc(1:le)
200     FORMAT('! CONSTRAINED SOLUTION for ',a)
     ELSEIF(incfit.eq.4)THEN
        WRITE(uniout,201)namloc(1:le)
201     FORMAT('! 4_PARAMETERS SOLUTION for ',a)
     ENDIF
  ENDIF
! COVARIANCE; HOWEVER NOT WRITTEN IF FORMAT= '1L'
  IF(PRESENT(covar))THEN
     IF(covar%succ)THEN
        defcov=.true.
        defnor=.true.
        cove(1:nd,1:nd)=covar%g(1:nd,1:nd)
        nore(1:nd,1:nd)=covar%c(1:nd,1:nd)
     ELSE
        defcov=.false.
        defnor=.false.
        cove=0.d0
        nore=0.d0
     ENDIF
  ELSE
     defcov=.false.
     defnor=.false.
     cove=0.d0
     nore=0.d0
  ENDIF
  mass=0.d0
  IF(form.eq.'1L')THEN
     CALL wro1lr(uniout,namloc(1:le),el%coord,el%coo,el%t,el%h_mag,el%g_mag)
  ELSEIF(form.eq.'ML')THEN
     CALL codestat(el%obscode,obscod)
     IF(PRESENT(dp))THEN
        CALL wromlr(uniout,namloc(1:le),el%coord,el%coo,el%t,cove(1:nd,1:nd), &
       &     defcov,nore(1:nd,1:nd),defnor,el%h_mag,el%g_mag,mass,nd,dp)
     ELSE
        CALL wromlr(uniout,namloc(1:le),el%coord,el%coo,el%t,cove(1:nd,1:nd), &
       &     defcov,nore(1:nd,1:nd),defnor,el%h_mag,el%g_mag,mass,nd)
     ENDIF
  ENDIF

  IF(.not.PRESENT(unit)) CALL filclo(uniout,' ')

END SUBROUTINE write_elems_old

! propagation to covariance/normal matrix of conversion of coordinates,
! with a 6x6 jacobian
SUBROUTINE convertunc(unc1,dee,unc2)
! unc1 is input, unc2 is output, but they can be on the same memory location
  TYPE(orb_uncert), INTENT(IN) :: unc1
  TYPE(orb_uncert), INTENT(INOUT) :: unc2
  DOUBLE PRECISION, DIMENSION(6,6), INTENT(IN) :: dee
!*******************************************************
  LOGICAl error
  INTEGER nd

! Dimension
  nd=unc1%ndim
! covariance matrix is propagated by similarity transformation
  CALL convertcovdp(nd,unc1%g(1:nd,1:nd),dee,unc2%g(1:nd,1:nd))
! normal matrix is propagated by similarity transformation  
  CALL norprsdp(unc1%c(1:nd,1:nd),dee,nd,unc2%c(1:nd,1:nd),error)
  IF(error)THEN
     unc2%succ=.false.
     unc2%ndim=unc1%ndim
  ELSE
     unc2%succ=.true.
     unc2%ndim=unc1%ndim
  ENDIF
END SUBROUTINE convertunc

! propagation of covariance/normal to a different time,
! with a nd x nd jacobian
SUBROUTINE propagunc(nd,unc1,depdep,unc2)
! unc1 is input, unc2 is output, but they can be on the same memory location
  INTEGER, INTENT(IN) :: nd
  TYPE(orb_uncert), INTENT(IN) :: unc1
  TYPE(orb_uncert), INTENT(INOUT) :: unc2
  DOUBLE PRECISION, DIMENSION(nd,nd), INTENT(IN) :: depdep
  LOGICAl error
! covariance matrix is propagated by similarity transformation
  CALL convertcov(nd,unc1%g(1:nd,1:nd),depdep(1:nd,1:nd),unc2%g(1:nd,1:nd))
! normal matrix is propagated by similarity transformation  
  CALL norprs(unc1%c(1:nd,1:nd),depdep(1:nd,1:nd),nd,unc2%c(1:nd,1:nd),error)
  IF(error)THEN
     unc2%succ=.false.
     unc2%ndim=unc1%ndim
  ELSE
     unc2%succ=.true.
     unc2%ndim=unc1%ndim
  ENDIF
END SUBROUTINE propagunc

CHARACTER*3 FUNCTION cootyp(ityp)
  INTEGER, INTENT(IN) :: ityp
  if(ityp.eq.1)then 
     cootyp='KEP' 
  elseif(ityp.eq.2)then 
     cootyp='EQU' 
  elseif(ityp.eq.3)then 
     cootyp='CAR'
  elseif(ityp.eq.4)then 
     cootyp='COM'
  elseif(ityp.eq.5)then 
     cootyp='COT'
  elseif(ityp.eq.6)then 
     cootyp='ATT' 
  else
     WRITE(*,*)' cootyp: coord code not known ',ityp
     STOP
  endif
END FUNCTION cootyp

! Copyright (C) 1997-1998 by Orbfit Consortium                          
! Version: December 16, 1998 Steven Chesley chesley@dm.unipi.it         
! Version 3.1 September 2003 A.Milani                               
! --------------------------------------------------------------------- 
!                                                                       
!  *****************************************************************    
!  *                                                               *    
!  *                         W R I K E P                           *    
!  *                                                               *    
!  *  Writes an orbital element record in an orbital element file  *    
!  *                     (multi-line format)                       *    
!  *                                                               *    
!  *****************************************************************    
!                                                                       
! INTPUT:   FILE      -  Output file name                               
!           NAME      -  Name of planet/asteroid/comet                  
!           ELEM      -  Orbital elements
!           COVAR     -  Uncertainty of orbital elements          
!           UNIT      -  If output file already opened
!
! OUPUT:    MOID      -  MOID                                           
!           PHA       -  0 or 1 PHA flag                                
!                                                                       
! WARNING: the routine does write the header of the file: this          
!          must not be generated by calling subroutine wromlh           
!                                                                       
SUBROUTINE wrikep(file,name,elem,covar,moid,pha,unit0)
  CHARACTER*(*), INTENT(INOUT) :: name
  TYPE(orbit_elem), INTENT(IN) :: elem
  TYPE(orb_uncert), INTENT(IN) :: covar
  DOUBLE PRECISION, INTENT(OUT) ::  moid 
  CHARACTER*(*), INTENT(INOUT) :: file 
  LOGICAL, INTENT(OUT) :: pha 
   INTEGER, OPTIONAL, INTENT(IN) :: unit0 
  TYPE(orbit_elem) :: elk
  TYPE(orb_uncert) :: covk,cove, unc6, unc6k
  INCLUDE 'parcmc.h90' 
  INTEGER lnn,lnf,i,j,k,unit, fail_flag
  DOUBLE PRECISION cnv(ndimx),std(ndimx)
  DOUBLE PRECISION enne, dkde(6,6), correl(ndimx,ndimx) 
  INTEGER :: nd
                                                                        
! "Other useful data"                                                   
  DOUBLE PRECISION perihe,aphe,period,dnode,anode,vinfty
  INTEGER iconv 

  nd=covar%ndim
  IF(PRESENT(unit0))THEN
     unit=unit0
  ELSE
! Open FIle    
     call rmsp(file,lnf) 
     call filopn(unit,file(1:lnf),'unknown') 
     CALL wromlh (unit,'ECLM','J2000') 
  ENDIF                                                                      
! Object Name     
  call rmsp(name,lnn) 
  WRITE(unit,100) name(1:lnn) 
100 FORMAT(A) 
! Get Keplerian Elements
  CALL coo_cha(elem,'KEP',elk,fail_flag,dkde) 
! WARNING: not working for cometary elements
! Write Orbital elements
  cnv(1:2)=1.d0
  cnv(3:6)=degrad 
  IF(nd.GT.6)THEN
     cnv(7:6+dyn%ndp)=1.d0
  END IF
  WRITE(unit,201) comcha 
201 FORMAT(A,' Keplerian elements: a, e, i, long. node,',  &
     &         ' arg. peric., mean anomaly')                            
  WRITE(unit,101) (elk%coord(i)*cnv(i),i=1,6) 
101 FORMAT(' KEP ',2F12.6,4F12.3) 
! Epoch
  WRITE(unit,104) elk%t 
104 FORMAT(' MJD ',F12.4,' TDT') 
! Mass                                                                  
!      IF(mass.NE.0.d0) WRITE(unit,105) mass 
!  105 FORMAT(' MAS ',1P,E20.12) 
                                                                        
! Magnitudes                                                            
  IF(elk%mag_set) WRITE(unit,106) elk%h_mag,elk%g_mag 
106 FORMAT(' MAG ',2F7.3) 
! nongrav parameters
  IF(dyn%nmod.ge.0)THEN
! ndp.gt.0 implies PRESENT(nmod) nmod.gt.0
! record giving the matrix dimension and list of solved parameters
     IF(nls.GT.0)THEN
        WRITE(unit,149) comcha
149     FORMAT(A,' Non-grav parameters: model used, actual number in use, dimension, ', &
             & 'list of solve for parameters to be currently determined')
        WRITE(unit,150) dyn%nmod, dyn%ndp, nd, ls(1:nls)
150     FORMAT(' LSP ',1X,I2,1X,I2,3X,I2,2X,4(1X,I2))
     ELSE
        WRITE(unit,'(A)') '! Non-grav parameters: model used, actual number in use, dimension'             
        WRITE(unit,155) dyn%nmod, dyn%ndp, nd
155     FORMAT(' LSP ',1X,I2,1X,I2,3X,I2)
     END IF
  END IF
! NGR
  IF(nd.GT.6)THEN
     WRITE(unit,112) ((dyn%dp(i))*cnv(6+i),i=1,dyn%ndp)
112  FORMAT(' NGR ',1P,4E12.3) !4E22.12
  END IF

! A-priori for non-gravitational parameter     
! ========================================
! Dynamical model nmod=1
     IF(dyn%nmod.EQ.1)THEN
! A/M and A2 to be determined
        IF(ls(1).EQ.1.AND.ls(2).EQ.2)THEN
! A/M and A2 a-priori available
           IF(dyn%apriori(2).AND.dyn%apriori(1))THEN
              WRITE(unit,'(A)') '! A priori value and standard deviation (A/M and Yarkovsky effect)'
              WRITE(unit,352) dyn%apriori_nominal(1), dyn%apriori_std(1)
              WRITE(unit,351) dyn%apriori_nominal(2), dyn%apriori_std(2)
! only A2 a-priori available
           ELSE IF(dyn%apriori(2))THEN
              WRITE(unit,'(A)') '! A priori value and standard deviation (only Yarkovsky effect)'
              WRITE(unit,351) dyn%apriori_nominal(2), dyn%apriori_std(2)
! only A/M a-priori available  
           ELSE IF(dyn%apriori(1))THEN
              WRITE(unit,'(A)') '! A priori value and standard deviation (only A/M)'
              WRITE(unit,352) dyn%apriori_nominal(1), dyn%apriori_std(1)
           END IF
! only A2 to be determined with a-priori
        ELSE IF(ls(1).EQ.2.AND.dyn%apriori(2))THEN 
           WRITE(unit,'(A)') '! A priori value and standard deviation (only Yarkovsky effect)'
           WRITE(unit,351) dyn%apriori_nominal(2), dyn%apriori_std(2)
! only A/M to be determined with a-priori
        ELSE IF(ls(1).EQ.1.AND.dyn%apriori(1))THEN    
           WRITE(unit,'(A)') '! A priori value and standard deviation (only A/M)'
           WRITE(unit,352) dyn%apriori_nominal(1), dyn%apriori_std(1)
        END IF

! Dynamical model nmod=2
     ELSEIF(dyn%nmod.EQ.2)THEN
        IF(ls(1).EQ.1.AND.ls(2).EQ.2.AND.ls(3).EQ.3.AND.ls(4).EQ.4)THEN
           IF(dyn%apriori(1).AND.dyn%apriori(2).AND.dyn%apriori(3).AND.dyn%apriori(4))THEN
              WRITE(unit,'(A)') '! A priori value and standard deviation (A1, A2, A3 and dtdelay)'
              WRITE(unit,353) dyn%apriori_nominal(1), dyn%apriori_std(1)
              WRITE(unit,354) dyn%apriori_nominal(2), dyn%apriori_std(2)
              WRITE(unit,355) dyn%apriori_nominal(3), dyn%apriori_std(3)
              WRITE(unit,356) dyn%apriori_nominal(4), dyn%apriori_std(4)
           END IF
        ELSE IF(ls(1).EQ.1.AND.ls(2).EQ.2.AND.ls(3).EQ.3)THEN
           IF(dyn%apriori(1).AND.dyn%apriori(2).AND.dyn%apriori(3))THEN
              WRITE(unit,'(A)') '! A priori value and standard deviation (A1, A2, A3)'
              WRITE(unit,353) dyn%apriori_nominal(1), dyn%apriori_std(1)
              WRITE(unit,354) dyn%apriori_nominal(2), dyn%apriori_std(2)
              WRITE(unit,355) dyn%apriori_nominal(3), dyn%apriori_std(3)
           END IF
        ELSE IF(ls(1).EQ.1.AND.ls(2).EQ.2)THEN
           IF(dyn%apriori(1).AND.dyn%apriori(2))THEN
              WRITE(unit,'(A)') '! A priori value and standard deviation (A1, A2)'
              WRITE(unit,353) dyn%apriori_nominal(1), dyn%apriori_std(1)
              WRITE(unit,354) dyn%apriori_nominal(2), dyn%apriori_std(2)
           END IF
        END IF
     END IF
351  FORMAT(' APR YAR ',1P,2E22.14)          
352  FORMAT(' APR A/M ',1P,2E22.14)
353  FORMAT(' APR A1C ',1P,2E22.14)
354  FORMAT(' APR A2C ',1P,2E22.14)
355  FORMAT(' APR A3C ',1P,2E22.14)
356  FORMAT(' APR DTD ',1P,2E22.14)

! Other useful data:                                                    
! PERIHELION 1.2345                                                     
! APHELION 2.3456                                                       
! ANODE 1.01234                                                         
! DNODE 1.04321                                                         
! MOID 0.00001                                                          
! PERIOD 1.11111                                                        
! PHA F                                                                 
! VINFTY 14.5678 (in km/s)
  perihe=elk%coord(1)*(1.d0-elk%coord(2)) 
  aphe=elK%coord(1)*(1.d0+elk%coord(2)) 
  enne=sqrt(gms/elk%coord(1)**3)
  period=dpig/enne 
  CALL nomoid(elem%t,elem,moid,anode,dnode) 
  pha = moid.le.0.05d0.and.elk%h_mag.lt.22d0.and.elk%mag_set
  vinfty= v_infty0(elk)
  WRITE(unit,102)comcha,perihe,comcha,aphe,comcha,                  &
     &     anode,comcha,dnode,comcha,moid,comcha,period,comcha,pha, &
     &     comcha,vinfty      
  102 FORMAT(A1,' PERIHELION ',F9.4/                                    &
     &       A1,' APHELION ',F9.4/                                      &
     &       A1,' ANODE ',F10.5/                                        &
     &       A1,' DNODE ',F10.5/                                        &
     &       A1,' MOID ',F10.5/                                         &
     &       A1,' PERIOD ',F16.4/                                       &
     &       A1,' PHA ',L1/                                             &
     &       A1,' VINFTY ',F9.4)
! Get Keplerian Covariance  
  IF(covar%succ) THEN 
     cove=covar
     IF(nd.GT.6)THEN
        unc6%c(1:6,1:6)=cove%c(1:6,1:6)
        unc6%g(1:6,1:6)=cove%g(1:6,1:6)
        unc6%succ=cove%succ
        unc6%ndim=6
        CALL convertunc(unc6,dkde,unc6k)
        covk%c(1:6,1:6)=unc6k%c(1:6,1:6)
        covk%g(1:6,1:6)=unc6k%g(1:6,1:6)
        covk%succ=unc6k%succ
        covk%g(1:nd,7:nd)=cove%g(1:nd,7:nd)
        covk%g(7:nd,1:nd)=cove%g(7:nd,1:nd)
        covk%c(1:nd,7:nd)=cove%c(1:nd,7:nd)
        covk%c(7:nd,1:nd)=cove%c(7:nd,1:nd)
     ELSE
        CALL convertunc(cove,dkde,covk)
     END IF
! Get Correlation Matrix                                                
     DO i=1,nd 
        std(i) = sqrt(covk%g(i,i)) 
     ENDDO
     DO i=1,nd 
        DO j=1,nd 
           correl(i,j) = covk%g(i,j)/std(i)/std(j) 
        ENDDO
     ENDDO                                                                   
! Covariance matrix                                                     
     WRITE(unit,107) comcha,(std(i)*cnv(i),i=1,nd) 
     WRITE(unit,108) ((covk%g(i,k)*cnv(i)*cnv(k),k=i,nd),i=1,nd) 
107  FORMAT(A1,' RMS ',1P,10E12.3) 
108  FORMAT(' COV ',1P,3E23.8) 
! Correlation matrix                                                    
     WRITE(unit,109)  ((correl(i,k),k=i,nd),i=1,nd) 
109  FORMAT(' COR ',3F23.8) 
  ENDIF
  IF(.not.PRESENT(unit0))THEN
     call filclo(unit,' ') 
  ENDIF
                                                                        
END SUBROUTINE wrikep

! ===================================================
! velocity at infinity with repect to Earth
! (circular approx) in au/day
! this is a duplicate of the one in tp_trace, to be used by wrikep
DOUBLE PRECISION FUNCTION v_infty0(el0)
! input elements
  TYPE(orbit_elem), INTENt(IN) :: el0
! END INTERFACE
  TYPE(orbit_elem) eleq
  DOUBLE PRECISION eq0(6) 
  DOUBLE PRECISION cosi,v2
  CHARACTER*3 coo
  INTEGER fail_flag
  coo=el0%coo
  CALL coo_cha(el0,'EQU',eleq,fail_flag)
  eq0=eleq%coord
  IF(fail_flag.eq.0)THEN 
! v_infinity computation in equinoctal
     cosi=(1.d0-eq0(4)**2-eq0(5)**2)/                                  &
     &     (1.d0+eq0(4)**2+eq0(5)**2)                                   
     v2=3.d0-1.d0/eq0(1) -                                             &
     &     2.d0*sqrt(eq0(1)*(1.d0-eq0(2)**2-eq0(3)**2))*cosi            
  ELSEIF(fail_flag.eq.5)THEN
! v_infinity computation in cometary
     cosi=cos(eq0(3))
     v2=3.d0-(1.d0-eq0(2))/eq0(1) - 2.d0*sqrt(eq0(1)*(1.d0+eq0(2)))*cosi
  ELSE
     WRITE(*,*)' v_infty: coordinate conversion failed ', el0,fail_flag  
  ENDIF
! handle non hyperbolic (w.r. to Earth) case
  IF(v2.gt.0.d0)THEN 
     v_infty0=sqrt(v2) 
  ELSE 
     v_infty0=0.d0 
  ENDIF
! normalization in au/day by Gauss constant
  v_infty0=v_infty0*gk
! conversion to km/s required by wrikep
   v_infty0=v_infty0*aukm/86400.d0
END FUNCTION v_infty0

! ===================================
! FIND_ELEMS
! find elements in .eq0 file, if necessary convert into
! required coordinate system, and report on availability
!  ===================================
SUBROUTINE find_elems(nam0,eledir,buildpath,coox,iunout,el0,unc0,ini0,cov0)
  CHARACTER*(*),INTENT(IN) :: nam0 ! asteroid name, could be long
  CHARACTER*60, INTENT(IN) ::  eledir ! root directory for .eq0 files
  LOGICAL, INTENT(IN):: buildpath ! to compose path with eledir, if not use 
                                  ! eledir as complete pathname  
  CHARACTER*3, INTENT(IN) :: coox ! in which coords the elements are required
   ! if coox=' ' coords are not changed, but left as they are
  INTEGER,INTENT(IN) :: iunout ! log output unit  
  TYPE(orbit_elem), INTENT(OUT) :: el0 ! epoch time (MJD), elements, abs. mag. 
  TYPE(orb_uncert), INTENT(OUT) :: unc0 ! covariance and normal matrices 
  LOGICAL,INTENT(OUT) :: ini0,cov0 ! successful input flag, covar. available
! === END INTERFACE =====================
  INCLUDE 'sysdep.h90'
  LOGICAL eof, error, err
  INTEGER nlsloc, lsloc(ndyx)
  INTEGER nd
  CHARACTER*60 elefi0 ! elements  file name
  INTEGER le, fail_flag
  CHARACTER*19 name
  CHARACTER*120 file
  TYPE(orbit_elem) :: elk
  TYPE(orb_uncert) :: unck 
  DOUBLE PRECISION :: dee(6,6)
  
!
  ini0=.false. 
  cov0=.false.
  IF(buildpath)THEN
     CALL fidinam(eledir,nam0,'eq0',elefi0,le)
  ELSE
     file=eledir//dircha//nam0//'.eq0'
     CALL rmsp(file,le)
     elefi0(1:le)=file(1:le)
  ENDIF  
  INQUIRE(file=elefi0(1:le),exist=ini0) 
  IF(ini0)THEN 
     CALL read_elems(elk,name,eof,nlsloc,err,lsloc,UNC=unck,FILE=elefi0(1:le))
     IF(dyn%nmod.GT.0)THEN
        nls=nlsloc
        ls=lsloc
     END IF
     ini0=(.not.eof).and.(name.eq.nam0) 
     IF(ini0)THEN
        cov0=unck%succ
        WRITE(iunout,*)nam0,' OK'
     ELSE
        WRITE(*,*)elefi0(1:le),'not read properly'
         CALL write_err(nam0,iunout,' ELEMENTS NOT FOUND')
     ENDIF
  ELSE 
     WRITE(*,*)elefi0(1:le),'not found'
     CALL write_err(nam0,iunout,' ELEMENT FILE NOT FOUND') 
  ENDIF
  nd=6+nls
  IF(coox.eq.' ')RETURN
  IF(ini0)THEN 
! coordinate change to equinoctal
     IF(cov0)THEN
        CALL coo_cha(elk,coox,el0,fail_flag,dee)
        CALL convertunc(unck,dee,unc0)
        cov0=unc0%succ
     ELSE                                       
        CALL coo_cha(elk,coox,el0,fail_flag)
     ENDIF
     IF(fail_flag.gt.5)THEN
        CALL write_err(nam0,iunout,' FAILED COORD CHANGE ')
        WRITE(ierrou,*) elk%coo,fail_flag,coox
        WRITE(iunout,*) elk%coo,fail_flag,coox
        ini0=.false.
        cov0=.false.
     ELSEIF(fail_flag.eq.5.and.(coox.eq.'EQU'.or.coox.eq.'KEP'))THEN
        CALL write_err(nam0,iunout,' HYPERBOLIC COORD CHANGE ')
        WRITE(ierrou,*) elk%coo,fail_flag,coox
        WRITE(iunout,*) elk%coo,fail_flag,coox
        ini0=.false.
        cov0=.false.
     ENDIF
!  no covariances available 
  ENDIF
END SUBROUTINE find_elems

! Copyright (C) 1997-1998 by Mario Carpino (carpino@brera.mi.astro.it)  
! Modified by Andrea Milani, vers. 1.8.3, January 1999                  
! additional changes for scaling, vers.3.3.1, August 2005
! --------------------------------------------------------------------- 
!                                                                       
!  *****************************************************************    
!  *                                                               *    
!  *                         W R O M L R                           *    
!  *                                                               *    
!  *  Writes an orbital element record in an orbital element file  *    
!  *                     (multi-line format)                       *    
!  *                                                               *    
!  *****************************************************************    
!                                                                       
! INPUT:   UNIT      -  Output FORTRAN unit                            
!           NAME0      -  Name of planet/asteroid/comet                  
!           ELEM(6)   -  Orbital element vector                         
!           ELTYPE    -  Type of orbital elements (KEP/EQU/CAR)         
!           T0        -  Epoch of orbital elements (MJD, TDT)           
!           COVE      -  Covariance matrix of orbital elements          
!           DEFCOV    -  Tells whether the covariance matrix is defined 
!           NORE      -  Normal matrix of orbital elements              
!           DEFNOR    -  Tells whether the normal matrix is defined     
!           H         -  H absolute magnitude (if <-100, missing)       
!           G         -  G slope parameter                              
!           MASS      -  Mass (solar masses)                            
!                                                                       
! WARNING: the routine does not write the header of the file: this      
!          must be generated by calling subroutine wromlh               
!                                                                       
SUBROUTINE wromlr(unit,name0,elem,eltype,t0,cove,defcov,nore,defnor,h,g,mass,n_par,dp)    
  USE fund_const  
  USE output_control
  USE name_rules  
  USE io_elems, ONLY: obscod  
  USE dyn_param
  IMPLICIT NONE
  INTEGER, INTENT(IN) :: unit 
  DOUBLE PRECISION, INTENT(IN) :: elem(6),t0,h,g,cove(n_par,n_par),nore(n_par,n_par),mass
  CHARACTER*(*), INTENT(IN):: name0,eltype 
  LOGICAL, INTENT(IN) ::  defcov,defnor 
  INTEGER, INTENT(IN) :: n_par
  DOUBLE PRECISION, INTENT(IN), OPTIONAL :: dp(ndyx)
! ------------END INTERFACE -------------------------------
  DOUBLE PRECISION gammas(ndimx,ndimx),scales(ndimx),units(ndimx) ! for scaling 
  CHARACTER*(idnamvir_len) name
  INCLUDE 'parcmc.h90'
  INTEGER l1,ln,i,k,iwd,j 
  DOUBLE PRECISION cnv(ndimx),std(ndimx), ele(6), princ
  INTEGER lench 
  EXTERNAL lench 
! eigenvalues, eigenvectors                                             
  DOUBLE PRECISION eigvec(ndimx,ndimx),eigval(ndimx),fv1(ndimx),fv2(ndimx),wdir(ndimx),sdir
  INTEGER ierr 
! Name
  name=' '
  name=name0                                                                  
  CALL rmsp(name,ln)
  WRITE(unit,100) name(1:ln) 
100 FORMAT(A)
! Orbital elements                                                      
  cnv=1.d0 
  ele=elem
  IF(eltype.EQ.'KEP') THEN 
     cnv(3:6)=degrad 
     ele(6)=princ(ele(6))
!     if(ele(6).lt.0.d0)ele(6)=ele(6)+dpig 
     if(ele(5).lt.0.d0)ele(5)=ele(5)+dpig 
     if(ele(4).lt.0.d0)ele(4)=ele(4)+dpig 
     WRITE(unit,201) comcha 
  201 FORMAT(A,' Keplerian elements: a, e, i, long. node,', &
           &         ' arg. peric., mean anomaly')  
     WRITE(unit,101) (ele(i)*cnv(i),i=1,6) 
101  FORMAT(' KEP ',1P,E24.16,0P,F18.15,4F18.13) 
  ELSEIF(eltype.EQ.'CAR') THEN 
     WRITE(unit,202) comcha 
202  FORMAT(A,' Cartesian position and velocity vectors') 
     WRITE(unit,102) ele 
102  FORMAT(' CAR ',1P,6E22.14) 
  ELSEIF(eltype.EQ.'EQU') THEN 
     cnv(6)=degrad 
     ele(6)=princ(ele(6))
     WRITE(unit,203) comcha 
  203 FORMAT(A,' Equinoctial elements: a, e*sin(LP), e*cos(LP),',       &
     &         ' tan(i/2)*sin(LN), tan(i/2)*cos(LN), mean long.') 
     WRITE(unit,103) (ele(i)*cnv(i),i=1,6) 
  103 FORMAT(' EQU ',1P,E22.14,0P,2(1x,f19.15),2(1x,f20.15),1x,F17.13) 
  ELSEIF(eltype.EQ.'COM') THEN 
     cnv(3:5)=degrad 
     if(ele(5).lt.0.d0)ele(5)=ele(5)+dpig 
     if(ele(4).lt.0.d0)ele(4)=ele(4)+dpig 
     WRITE(unit,204) comcha 
  204 FORMAT(A,' Cometary elements: q, e, i, long. node,', &
           &         ' arg. peric., pericenter time')  
     IF(abs(ele(6)).lt.1.d6)THEN
        WRITE(unit,114) (ele(i)*cnv(i),i=1,6) 
114     FORMAT(' COM ',1P,E22.14,0P,1x,F18.15,1x,3(F18.13,1x),F18.10)
     ELSE
        WRITE(unit,215) (ele(i)*cnv(i),i=1,6) 
215     FORMAT(' COM ',1P,E22.14,0P,1x,F18.15,1x,3(F18.13,1x),1P,D18.11)
     ENDIF 
  ELSEIF(eltype.EQ.'COT') THEN 
     cnv(3:6)=degrad 
     if(ele(5).lt.0.d0)ele(5)=ele(5)+dpig 
     if(ele(4).lt.0.d0)ele(4)=ele(4)+dpig
     if(ele(6).lt.0.d0)ele(6)=ele(6)+dpig
     WRITE(unit,206) comcha 
  206 FORMAT(A,' Cometary elements: q, e, i, long. node,', &
           &         ' arg. peric., true anomaly')  
     WRITE(unit,117) (ele(i)*cnv(i),i=1,6) 
117  FORMAT(' COT ',1P,E22.14,0P,1x,F18.15,1x,4(F18.13,1x))
  ELSEIF(eltype.eq.'ATT') THEN
     cnv(1:4)=degrad
     IF(ele(1).lt.0.d0)ele(1)=ele(1)+dpig
          WRITE(unit,205) comcha 
  205 FORMAT(A,' Attributable elements: R.A., DEC, R.A.dot, DECdot, r, rdot,', &
           &         ' auxiliary information on observer needed')  
     WRITE(unit,115) (ele(i)*cnv(i),i=1,6) 
115  FORMAT(' ATT ',2(F18.13,1x),2(F18.13,1x),2(F18.13,1x)) 
     WRITE(unit,116) ' STA ',obscod
116  FORMAT(A5,A3)
  ELSE
     WRITE(*,*)  '**** wromlr: unsupported orbital element type ****', eltype
     STOP '**** wromlr: unsupported orbital element type ****' 
  END IF

! Epoch                                                                 
  WRITE(unit,104) t0 
104 FORMAT(' MJD ',F19.9,' TDT')                                                                     
! Mass                                                                  
  IF(mass.NE.0.d0) WRITE(unit,105) mass 
105 FORMAT(' MAS ',1P,E20.12)                                              
! Magnitudes                                                            
  IF(h.GT.-100.d0) WRITE(unit,106) h,g 
106 FORMAT(' MAG ',2F7.3)                                                 

! nongrav parameters
  IF(n_par.gt.6)THEN
     cnv(7:6+dyn%ndp)=1.d0
     IF(PRESENT(dp))THEN
        WRITE(unit,112) (dp(i)*cnv(6+i),i=1,dyn%ndp)
     ELSE
        WRITE(unit,112) (dp(i)*cnv(6+i),i=1,dyn%ndp)
     ENDIF 
112  FORMAT(' NGR ',1P,4E22.14) 
  ENDIF


! Covariance matrix                                                     
  IF(defcov) THEN 
     DO i=1,n_par 
        std(i)=SQRT(cove(i,i))
     ENDDO 
! eigenvalues
     gammas(1:n_par,1:n_par)=cove(1:n_par,1:n_par)
     CALL weak_dir(gammas(1:n_par,1:n_par),wdir(1:n_par),sdir,-1,eltype,ele,units(1:n_par),n_par)
! Eigenvalues have not to be rescaled
!     scales(1:n_par)=1.d0/units(1:n_par)
!     DO i=1,n_par
!        DO j=1,n_par
!           gammas(i,j)=gammas(i,j)*(scales(i)*scales(j))
!        ENDDO
!     ENDDO
     CALL rs(n_par,n_par,gammas(1:n_par,1:n_par),eigval(1:n_par),1,eigvec(1:n_par,1:n_par),fv1,fv2,ierr) 
     DO i=1,n_par 
        IF(eigval(i).gt.0.d0)THEN 
           eigval(i)=sqrt(eigval(i)) 
        ELSE 
           IF(eigval(i).lt.-1.d-10.and.verb_io.gt.9)THEN
              WRITE(*,*)'wromlr: zero/negative eigenvalue', eigval(i),'  for asteroid ',name(1:ln)
           ENDIF 
           eigval(i)=-sqrt(-eigval(i)) 
        ENDIF
     ENDDO
! RMS, eigenvalues and weak direction are so far commented   
! format to be changed for ndimx>10           
     WRITE(unit,107) comcha,(std(i)*cnv(i),i=1,n_par) 
107  FORMAT(A1,' RMS ',1P,10E14.5) 
     WRITE(unit,111) comcha,eigval(1:n_par)
111  FORMAT(A1,' EIG',1P,10E14.5) 
     WRITE(unit,110) comcha,wdir(1:n_par)
110  FORMAT(A1,' WEA',10F10.5) 
! covariance matrix is given uncommented, to be readable                
     WRITE(unit,108) ((cove(i,k)*cnv(i)*cnv(k),k=i,n_par),i=1,n_par) 
108  FORMAT(' COV ',1P,3E23.15) 
  END IF                                                                      
! Normal matrix                                                         
  IF(defnor) THEN 
     WRITE(unit,109) ((nore(i,k)/(cnv(i)*cnv(k)),k=i,n_par),i=1,n_par) 
109  FORMAT(' NOR ',1P,3E23.15) 
  END IF
END SUBROUTINE wromlr

! --------------------------------------------------------------------- 
!                                                                       
!  *****************************************************************    
!  *                                                               *    
!  *                          R D _ O R B                            *    
!  *                                                               *    
!  *     Read orbital elements from a file opened with OPORBF      *    
!  *                                                               *    
!  *****************************************************************    
!                                                                       
! The routine operates in a sequential way, returning each time         
! the orbital elements of the next object contained in the file         
!                                                                       
! OUTPUT:   NAME      -  Name of planet/asteroid/comet                  
!           ELEM(6)   -  Orbital element vector                         
!           ELTYPE    -  Type of orbital elements (KEP/EQU/CAR/COM/ATT)
!           T0        -  Epoch of orbital elements (MJD, TDT)           
!           COVE      -  Covariance matrix of orbital elements          
!           DEFCOV    -  Tells whether the covariance matrix is defined 
!           NORE      -  Normal matrix of orbital elements              
!           DEFNOR    -  Tells whether the normal matrix is defined     
!           H         -  H absolute magnitude (if <-100, missing)       
!           G         -  G slope parameter                              
!           MASS      -  Mass (solar masses)                            
!           RSYS      -  Reference system type (EQUM/EQUT/ECLM)         
!           EPOCH     -  Epoch specification (J2000/OFDATE)             
!           KR        -  Record number at which object is found         
!           EOF       -  End-of-file flag                               
!                                                                       
 SUBROUTINE rd_orb(name,elem,eltype,t0,cove,defcov,nore,defnor,     &
     &                 h,g,mass,rsys,epoch,kr,eof,dp,nd) 
   USE fund_const
   USE io_elems
   IMPLICIT NONE
   DOUBLE PRECISION elem(6),t0,h,g,cove(ndimx,ndimx),nore(ndimx,ndimx),mass 
   DOUBLE PRECISION,DIMENSION(ndyx),OPTIONAL :: dp ! gravitational parameters
   INTEGER,OPTIONAL :: nd
   CHARACTER*(*) name,eltype,rsys,epoch 
   LOGICAL defcov,defnor,eof 
   INTEGER kr
   CHARACTER rec*200,rest*200,scale*3 
   INTEGER lf,nit,i,k,mjd,mjde,ik 
   DOUBLE PRECISION sec,sece,tmp(ndimx*(ndimx+1)/2),cnv(ndimx) 
   LOGICAL error,end1,noep
   INTEGER lench,nitchs 
! number of elements in upper part of a symmetric matrix
! number of rows to be read
   INTEGER :: nel,nrows
! number of solve for parameters
   INTEGER :: nd_cur
! number of solve for parameters
   EXTERNAL lench,nitchs                                                       
   IF(iicorb.NE.36) STOP '**** rd_orb: internal error (01) ****' 
   IF(orbunt.LE.0) STOP '**** rd_orb: internal error (02) ****' 
! if present dp, then nd is required
  IF(PRESENT(dp).and..not.PRESENT(nd))THEN
     WRITE(*,*)'rd_orb: nd not present'
     STOP
  ENDIF
! number of rows for normal/covariance matrix
   IF(PRESENT(nd))THEN
      nd_cur=nd
   ELSE
      nd_cur=6
   ENDIF
   nel=nd_cur*(nd_cur+1)/2
   nrows=nel/3
   IF(MOD(nel,3).NE.0)THEN
      nrows=nrows+1
   ENDIF
   rsys=dfrsty 
   epoch=dfrsep 
   mass=0.d0 
   DO 1 i=1,ndimx 
      DO 2 k=1,ndimx 
         cove(i,k)=0 
         nore(i,k)=0 
2     END DO
1  END DO
   defcov=.false. 
   defnor=.false. 
   h=-1.d9
   g=-1.d9
   IF(rectyp.EQ.'1L') THEN 
      CALL getrsc(orbunt,rec,orbnr,eof) 
      IF(eof) RETURN 
      CALL strcnt(rec,name,rest,error) 
      IF(error) GOTO 20 
      kr=orbnr 
      nit=nitchs(rest) 
      IF(deft0) THEN 
         t0=dept0 
         IF(nit.LT.8) THEN 
            READ(rest,*,ERR=20) elem 
         ELSE 
            READ(rest,*,ERR=20) elem,h,g 
         END IF
      ELSE 
         IF(nit.LT.9) THEN 
            READ(rest,*,ERR=20) t0,elem 
         ELSE 
            READ(rest,*,ERR=20) t0,elem,h,g 
         END IF
      END IF
      eltype=deltyp 
   ELSEIF(rectyp.EQ.'ML') THEN 
      IF(nxtend) THEN 
         eof=.true. 
         RETURN 
      END IF
      noep=.true. ! no epoch 
! Name                                                                  
      CALL getrsc(orbunt,name,orbnr,eof) 
      IF(eof) RETURN 
      kr=orbnr 
! Orbital elements (mandatory, immediately after the name)              
      CALL getrsc(orbunt,rec,orbnr,end1) 
      IF(end1) GOTO 20
      IF(rec(1:1).ne.' ') GOTO 20 
      eltype=rec(2:4)
      READ(rec(5:),*,ERR=20) elem
! Other keywords                                                        
3     CONTINUE 
      CALL getrsc(orbunt,rec,orbnr,end1) 
      IF(end1) THEN 
         nxtend=.true. 
         GOTO 4 
      END IF
      IF(rec(1:1).NE.' ') THEN 
         BACKSPACE(orbunt) 
         orbnr=orbnr-1 
         GOTO 4 
      END IF
! Epoch of elements                                                     
      IF(rec(1:4).EQ.' MJD' .OR. rec(1:4).EQ.' JD ' .OR.            &
           &       rec(1:4).EQ.' CAL') THEN                                   
         CALL ch2tim(rec,mjd,sec,scale,error) 
         IF(error) GOTO 20 
         CALL cnvtim(mjd,sec,scale,mjde,sece,'TDT') 
         t0=mjde+sece/86400.d0 
         noep=.false. 
      ELSEIF(rec(1:4).EQ.' MAG') THEN 
         READ(rec(5:),*,ERR=20) h,g 
      ELSEIF(rec(1:4).EQ.' STA')THEN ! station code for ATT only
         READ(rec(5:),*,ERR=20) obscod
      ELSEIF(rec(1:4).EQ.' MAS') THEN 
         READ(rec(5:),*,ERR=20) mass 
      ELSEIF(rec(1:4).EQ.' NGR') THEN
         IF(PRESENT(dp))THEN
            READ(rec(5:),*,ERR=20) dp(1:dyn%ndp)
         ELSE
            WRITE(*,*)'rd_orb: dp not present, cannot read non gravitational parameters'
            STOP
         ENDIF
      ELSEIF(rec(1:4).EQ.' COV') THEN 
         READ(rec(5:),*,ERR=40) (tmp(i),i=1,3) 
         DO 17 k=1,nrows-1
            CALL getrsc(orbunt,rec,orbnr,end1) 
            IF(end1) GOTO 40 
            IF(rec(1:4).NE.' COV') GOTO 40
! the last row can be different if the number of elements if MOD(nel,3) NE 0
            IF(k.EQ.nrows-1.AND.MOD(nel,3).NE.0)THEN
               READ(rec(5:),*,ERR=20) (tmp(i),i=3*k+1,3*k+MOD(nel,3)) 
            ELSE
               READ(rec(5:),*,ERR=20) (tmp(i),i=3*k+1,3*k+3) 
            ENDIF
17       CONTINUE 
         ik=0 
         DO 8 i=1,nd_cur
            DO 7 k=i,nd_cur 
               ik=ik+1 
               cove(i,k)=tmp(ik) 
7           ENDDO 
8        ENDDO 
         IF(ik.NE.nel) STOP '**** rd_orb: internal error (03) ****' 
         defcov=.true. 
      ELSEIF(rec(1:4).EQ.' NOR') THEN 
         READ(rec(5:),*,ERR=40) (tmp(i),i=1,3) 
         DO 27 k=1,nrows-1
            CALL getrsc(orbunt,rec,orbnr,end1) 
            IF(end1) GOTO 40 
            IF(rec(1:4).NE.' NOR') GOTO 40 
            IF(k.EQ.nrows-1.AND.MOD(nel,3).NE.0)THEN
               READ(rec(5:),*,ERR=20) (tmp(i),i=3*k+1,3*k+MOD(nel,3))
            ELSE
               READ(rec(5:),*,ERR=40) (tmp(i),i=3*k+1,3*k+3) 
            ENDIF
   27    ENDDO 
         ik=0 
         DO 38 i=1,nd_cur
            DO 37 k=i,nd_cur 
               ik=ik+1 
               nore(i,k)=tmp(ik) 
37          ENDDO
38       ENDDO
         IF(ik.NE.nel) STOP '**** rd_orb: internal error (04) ****' 
         defnor=.true. 
      ELSE 
         GOTO 40 
      END IF 
      GOTO 3 

40    CONTINUE
      lf=lench(orbfn) 
      WRITE(*,240) orbfn(1:lf),orbnr 
240   FORMAT(' ERROR in covariance, file ',A,' at line',I6) 

4     CONTINUE 
      IF(noep) THEN 
         IF(deft0) THEN 
            t0=dept0 
         ELSE 
            GOTO 20 
         END IF
      END IF
   ELSE 
      STOP '**** rd_orb: internal error (05) ****' 
   END IF
! Transformation of angles in orbital elements                          
! and covariance matrix                                                 
   cnv(1:nd_cur)=1 
   IF(eltype.EQ.'KEP') THEN 
      cnv(3:6)=radeg 
   ELSEIF(eltype.EQ.'COM') THEN 
      cnv(3:5)=radeg 
   ELSEIF(eltype.EQ.'COT')THEN
      cnv(3:6)=radeg
   ELSEIF(eltype.EQ.'EQU') THEN 
      cnv(6)=radeg 
   ELSEIF(eltype.EQ.'CAR') THEN 
      CONTINUE 
   ELSEIF(eltype.EQ.'ATT') THEN 
      cnv(1:4)=radeg 
   ELSE 
      STOP '**** rd_orb: internal error (06) ****' 
   END IF
   DO  i=1,6 
      elem(i)=elem(i)*cnv(i) 
   END DO 
   IF(defcov) THEN 
      DO 13 i=1,nd_cur 
         cove(i,i)=cove(i,i)*(cnv(i)**2) 
         DO 12 k=i+1,nd_cur 
            cove(i,k)=cove(i,k)*cnv(i)*cnv(k) 
            cove(k,i)=cove(i,k) 
12       ENDDO
13    ENDDO
   END IF
   IF(defnor) THEN 
      DO 15 i=1,nd_cur
         nore(i,i)=nore(i,i)/(cnv(i)**2) 
         DO 14 k=i+1,nd_cur 
            nore(i,k)=nore(i,k)/(cnv(i)*cnv(k)) 
            nore(k,i)=nore(i,k) 
14       ENDDO
15    ENDDO
   END IF                                                                
   eof=.false. 
   RETURN                                                                      
20 CONTINUE 
   lf=lench(orbfn) 
   WRITE(*,200) orbfn(1:lf),orbnr 
200 FORMAT(' ERROR in file ',A,' at line',I6) 
   STOP '**** rd_orb: abnormal end ****' 
END SUBROUTINE rd_orb

! upgrade to a consistent convention for handling non-gravitational
! perturbations
! ==================================================================
! READ_ELEMS
! New version of read_elemns to handle non-gravitational parameters
! Reads elements and dynamical parameters, if present, by using 
! only information from the file. The size of the covar/normal
! matrices is determined by the LSP record. The dynamical parameters dp are
! output if given in the file. 
! ==================================================================
SUBROUTINE read_elems(el,name,eof,nlsloc,err,lsloc,form,unc,file,unit,mass,rsys,epoch)
  USE station_coordinates, ONLY: statcode
  USE io_elems
  USE planet_masses, ONLY: icrel
! ====BEGIN INTERFACE ==================================
  TYPE(orbit_elem),           INTENT(OUT)   :: el       ! elements
  CHARACTER*(*),              INTENT(OUT)   :: name     ! object name (read from input file)
  LOGICAL,                    INTENT(OUT)   :: eof      ! end of file
  INTEGER,                    INTENT(OUT)   :: nlsloc   ! number of effective solve-for par 
  LOGICAL,                    INTENT(OUT)   :: err      ! error flag
  ! optional
  INTEGER,                    INTENT(OUT)   :: lsloc(ndyx) ! solve-for dynamical parameters
  CHARACTER*(*),    OPTIONAL, INTENT(INOUT) :: form        ! can be either '1L' or 'ML' (single/multi line)
  TYPE(orb_uncert), OPTIONAL, INTENT(OUT)   :: unc         ! covariance/normal matrices 
  CHARACTER*(*),    OPTIONAL, INTENT(IN)    :: file        ! input file (not used if PRESENT(unit))
  INTEGER,          OPTIONAL, INTENT(IN)    :: unit        ! input unit, only if already opened
  DOUBLE PRECISION, OPTIONAL, INTENT(OUT)   :: mass        ! mass of the object 
  CHARACTER*10,     OPTIONAL, INTENT(OUT)   :: rsys,epoch
  !  INTEGER, INTENT(OUT), OPTIONAL :: incfit ! result from incomplete fit;
! WARNING: if file already opened, it must also have the header already written
! ====END INTERFACE ==================================
  INTEGER          :: nd                                      ! dimension of covar/normal matrices 
  INTEGER          :: inunit                                  ! input orbit file (unit) 
  INTEGER          :: inunitngr                               ! input 1-line catalog with non-grav param
  INTEGER          :: le, i, k, j, ik, lench, lf              ! integers: loop indices and length of files
  CHARACTER*2      :: formcur
  CHARACTER*50     :: namloc
  CHARACTER*6      :: vers                                    ! version (we need to switch to the new version 2.0)
  LOGICAL          :: error, eop, noep, end1                  ! various case of failure
  CHARACTER        :: rec*200,rest*200,scale*3                ! record, full and left
  INTEGER          :: nit,nitchs                   
  INTEGER          :: kr
  INTEGER          :: mjd,mjde                                ! time (mjd)
  DOUBLE PRECISION :: sec,sece, tmp(ndimx*(ndimx+1)/2) 
  DOUBLE PRECISION :: cnv(ndimx),std(ndimx), princ, mass_loc
  DOUBLE PRECISION :: cove(ndimx,ndimx),nore(ndimx,ndimx)     ! covariance and normal matrices
  TYPE(dyn_par)    :: dyn_loc
  INTEGER          :: ndloc
  INTEGER          :: sum_active                              ! sum of the active non-grav parameters
! number of elements in upper part of a symmetric matrix
! number of rows to be read
  INTEGER :: nel,nrows
  INCLUDE 'parcmc.h90'
!--------------------------------------------------------
! Initizialing
  err=.false.;eof=.false.
  dyn_loc = undefined_dyn_par
!  IF(PRESENT(lsloc))THEN
  lsloc=0
!  ENDIF
! to avoid unassigned variable
  IF(PRESENT(mass))THEN
     mass=0.d0
  ENDIF
! open output file if not opened already
  IF(.not.PRESENT(unit))THEN
     CALL oporbf(file,0)
     vers=oef_vers
     IF(PRESENT(form))THEN
        form=rectyp
     ENDIF
     inunit=orbunt
     IF(orbunt.LE.0)THEN
        WRITE(*,*) ' read_elems: failed oporbf, orbunt= ',orbunt
        STOP ' read_elems: failed oporbf'
     ENDIF
  ELSE
     inunit=unit 
     vers=oef_vers
  ENDIF
  IF(PRESENT(rsys).AND.PRESENT(epoch))THEN
     rsys=dfrsty 
     epoch=dfrsep 
  END IF
  IF(vers.EQ.'OEF1.1')THEN
     nd=6
     nlsloc=nd-6
!     IF(PRESENT(lsloc))THEN
        lsloc=0
!     ENDIF
  END IF
! which record type?
  IF(rectyp.eq.'1L')THEN
     el=undefined_orbit_elem
     CALL getrsc(orbunt,rec,orbnr,eof) 
     IF(eof) RETURN 
     CALL strcnt(rec,name,rest,error) 
     IF(error) GOTO 20 
     kr=orbnr 
     nit=nitchs(rest) 
     IF(nit.LT.9) THEN 
! no magnitude available
        READ(rest,*,ERR=20) el%t,el%coord 
        el%mag_set=.false.
     ELSEIF(nit.EQ.10)THEN 
! include magnitude H,G
        READ(rest,*,ERR=20) el%t,el%coord,el%h_mag,el%g_mag
        el%mag_set=.true.
! old version, no dyn-param
     ELSEIF(nit.EQ.11.AND.vers.EQ.'OEF1.1')THEN
        ! include magnitude H,G
        ! No dyn parameters
        READ(rest,*,ERR=20) el%t,el%coord,el%h_mag,el%g_mag
        el%mag_set=.TRUE.
! new version: non-parameters
     ELSEIF(nit.EQ.11.AND.vers.EQ.'OEF2.0')THEN
        READ(rest,*,ERR=20) el%t,el%coord,el%h_mag,el%g_mag,sum_active
        el%mag_set=.TRUE.
     ELSE
        WRITE(*,*)'read_elems: one magnitude missing, nit=',nit
        READ(rest,*,ERR=20) el%t,el%coord,el%h_mag
        el%mag_set=.true.
     END IF
     el%coo=deltyp 
  ELSEIF(rectyp.eq.'ML')THEN
      IF(nxtend)THEN
        eof=.true.
        RETURN
     ENDIF
! setup of number of parameters in the default way, without nongrav
     el=undefined_orbit_elem
! initialize unc
     IF(PRESENT(unc))CALL undefined_orb_uncert(ndimx,unc)
     noep=.true. ! no epoch 
! Name                                                                  
     CALL getrsc(orbunt,name,orbnr,eof) 
     IF(eof) GOTO 4 
     kr=orbnr 
! Orbital elements (mandatory, immediately after the name)              
     CALL getrsc(orbunt,rec,orbnr,end1) 
     IF(end1) GOTO 20
     IF(rec(1:1).ne.' ') GOTO 20 
     el%coo=rec(2:4)
     READ(rec(5:),*,ERR=20) el%coord
     ! Other keywords                                                        
     DO j=1,100 
        CALL getrsc(orbunt,rec,orbnr,end1) 
        IF(end1) THEN 
           nxtend=.true. 
           IF(noep) GOTO 20 ! epoch time missing
           GOTO 4 !end of file, but normalize the last
        END IF
        IF(rec(1:1).NE.' ') THEN 
! next object, needs to be unread
           BACKSPACE(orbunt) 
           orbnr=orbnr-1
           IF(noep) GOTO 20 ! epoch time missing
           GOTO 4  !end of one object, normalize
        END IF
! Epoch of elements                                                     
        IF(rec(1:4).EQ.' MJD' .OR. rec(1:4).EQ.' JD ' .OR.            &
             &       rec(1:4).EQ.' CAL') THEN   
           CALL ch2tim(rec,mjd,sec,scale,error) 
           IF(error) GOTO 20 
           CALL cnvtim(mjd,sec,scale,mjde,sece,'TDT') 
           el%t=mjde+sece/86400.d0 
           noep=.false. ! epoch time available
        ELSEIF(rec(1:4).EQ.' MAG') THEN 
           READ(rec(5:),*,ERR=20) el%h_mag,el%g_mag
           el%mag_set=.true. 
        ELSEIF(rec(1:4).EQ.' STA')THEN ! station code for ATT only
           READ(rec(5:),*,ERR=20) obscod
           IF(el%coo.eq.'ATT') THEN 
              CALL statcode(obscod, el%obscode)
           ENDIF
        ELSEIF(rec(1:4).EQ.' MAS') THEN 
! problem: this mass is not used
           READ(rec(5:),*,ERR=20) mass_loc 
           IF(PRESENT(mass))THEN
              mass=mass_loc
           END IF
! nongrav handling 
        ELSEIF(rec(1:4).EQ.' LSP')THEN 
           nit=nitchs(rec(5:))
           IF(nit-3.GE.0.AND.nit-3.LE.ndyx)THEN
              READ(rec(5:),*,ERR=20)dyn_loc%nmod, dyn_loc%ndp, nd, lsloc(1:nit-3)
              nlsloc=nd-6
            ELSE IF(nit-3.LT.0) THEN
              WRITE(*,*) 'read_elems: error in reading LSP field: ', rec(5:)
              STOP
           ENDIF
        ELSEIF(rec(1:4).EQ.' NGR') THEN
           READ(rec(5:),*,ERR=20) dyn_loc%dp(1:dyn_loc%ndp)
           IF(dyn_loc%dp(1).NE.0.d0)THEN
              dyn_loc%active(1)=.TRUE.
           END IF
           IF(dyn_loc%dp(2).NE.0.d0)THEN
              dyn_loc%active(2)=.TRUE.
              IF(dyn_loc%nmod.EQ.1)THEN
                 icrel=2
              END IF
           END IF
           IF(dyn_loc%dp(3).NE.0.d0)THEN
              dyn_loc%active(3)=.TRUE.
           END IF
           IF(dyn_loc%dp(4).NE.0.d0)THEN
              dyn_loc%active(4)=.TRUE.
           END IF
        ELSEIF(rec(1:9).EQ.' APR YAR ') THEN
           dyn_loc%apriori(2)=.TRUE.
           READ(rec(10:),*,ERR=20) dyn_loc%apriori_nominal(2), dyn_loc%apriori_std(2)
        ELSEIF(rec(1:9).EQ.' APR A/M ') THEN
           dyn_loc%apriori(1)=.TRUE.
           READ(rec(10:),*,ERR=20) dyn_loc%apriori_nominal(1), dyn_loc%apriori_std(1)
        ELSEIF(rec(1:9).EQ.' APR A1C ') THEN
           dyn_loc%apriori(1)=.TRUE.
           READ(rec(10:),*,ERR=20) dyn_loc%apriori_nominal(1), dyn_loc%apriori_std(1)
        ELSEIF(rec(1:9).EQ.' APR A2C ') THEN
           dyn_loc%apriori(2)=.TRUE.
           READ(rec(10:),*,ERR=20) dyn_loc%apriori_nominal(2), dyn_loc%apriori_std(2)
        ELSEIF(rec(1:9).EQ.' APR A3C ') THEN
           dyn_loc%apriori(3)=.TRUE.
           READ(rec(10:),*,ERR=20) dyn_loc%apriori_nominal(3), dyn_loc%apriori_std(3)
        ELSEIF(rec(1:9).EQ.' APR DTD ') THEN
           dyn_loc%apriori(4)=.TRUE.
           READ(rec(10:),*,ERR=20) dyn_loc%apriori_nominal(4), dyn_loc%apriori_std(4)
        ELSEIF(rec(1:4).EQ.' COV') THEN 
! covariance and normal matrices: size
           nel=nd*(nd+1)/2
           nrows=nel/3
! initialize cove and nore
           DO  i=1,ndimx 
              DO  k=1,ndimx 
                 cove(i,k)=0 
                 nore(i,k)=0 
              END DO
           END DO
! initialize unc
           IF(PRESENT(unc))CALL undefined_orb_uncert(nd,unc)
           IF(MOD(nel,3).NE.0)THEN
              nrows=nrows+1
           ENDIF
           READ(rec(5:),*,ERR=40) (tmp(i),i=1,3) 
           DO  k=1,nrows-1
              CALL getrsc(orbunt,rec,orbnr,end1) 
              IF(end1) GOTO 40
              IF(rec(1:4).NE.' COV') GOTO 40
! the last row can be different if the number of elements if MOD(nel,3) NE 0
              IF(k.EQ.nrows-1.AND.MOD(nel,3).NE.0)THEN
                 READ(rec(5:),*,ERR=20) (tmp(i),i=3*k+1,3*k+MOD(nel,3)) 
              ELSE
                 READ(rec(5:),*,ERR=20) (tmp(i),i=3*k+1,3*k+3) 
              ENDIF
           ENDDO
           ik=0 
           DO i=1,nd
              DO  k=i,nd 
                 ik=ik+1 
                 cove(i,k)=tmp(ik) 
              ENDDO
           ENDDO
           IF(ik.NE.nel) STOP 'read_elems: insufficient COV records'
        ELSEIF(rec(1:4).EQ.' NOR') THEN
! assumed NOR is after COV, thus size has been computed 
           READ(rec(5:),*,ERR=40) (tmp(i),i=1,3) 
           DO k=1,nrows-1
              CALL getrsc(orbunt,rec,orbnr,end1) 
              IF(end1) GOTO 40 
              IF(rec(1:4).NE.' NOR') GOTO 40 
              IF(k.EQ.nrows-1.AND.MOD(nel,3).NE.0)THEN
                 READ(rec(5:),*,ERR=20) (tmp(i),i=3*k+1,3*k+MOD(nel,3))
              ELSE
                 READ(rec(5:),*,ERR=40) (tmp(i),i=3*k+1,3*k+3) 
              ENDIF
           ENDDO
           ik=0 
           DO i=1,nd
              DO k=i,nd 
                 ik=ik+1 
                 nore(i,k)=tmp(ik) 
              ENDDO
           ENDDO
           IF(ik.NE.nel) STOP 'read_elems: insufficient NOR records'
! covariance is available, but copied lateafter rescaling
           IF(PRESENT(unc)) unc%succ=.true.          
        ELSE 
           WRITE(*,*) 'read_elems: unrecognized keyword =',rec(1:4)
           STOP ' unrecognized record'  
        END IF
! potentially infinite loop removed
     ENDDO

40   CONTINUE
     lf=lench(orbfn) 
     WRITE(*,240) orbfn(1:lf),orbnr 
240  FORMAT(' ERROR in covariance, file ',A,' at line',I6) 
  ELSE
     WRITE(*,*) ' read_elems: unknown format, form=',form
     STOP ' orbit format?? '
  ENDIF
4 CONTINUE
! elements and possibly matrices read, now rescale
! Transformation of angles in orbital elements                          
! and covariance matrix   
  cnv(1:nd)=1 
   IF(el%coo.EQ.'KEP') THEN 
      cnv(3:6)=radeg 
   ELSEIF(el%coo.EQ.'COM') THEN 
      cnv(3:5)=radeg 
   ELSEIF(el%coo.EQ.'COT')THEN
      cnv(3:6)=radeg
   ELSEIF(el%coo.EQ.'EQU') THEN 
      cnv(6)=radeg 
   ELSEIF(el%coo.EQ.'CAR') THEN 
      CONTINUE 
   ELSEIF(el%coo.EQ.'ATT') THEN 
      cnv(1:4)=radeg 
   ELSE 
      WRITE(*,*)'read_elems: unknown coordinate type=',el%coo
      STOP 'unknown coordinate type' 
   ENDIF
   DO  i=1,6 
      el%coord(i)=el%coord(i)*cnv(i) 
   END DO
   IF(PRESENT(unc))THEN 
      IF(unc%succ)THEN 
         DO i=1,nd 
            cove(i,i)=cove(i,i)*(cnv(i)**2) 
            DO k=i+1,nd 
               cove(i,k)=cove(i,k)*cnv(i)*cnv(k) 
               cove(k,i)=cove(i,k) 
            ENDDO
         ENDDO
         DO i=1,nd
            nore(i,i)=nore(i,i)/(cnv(i)**2) 
            DO k=i+1,nd 
               nore(i,k)=nore(i,k)/(cnv(i)*cnv(k)) 
               nore(k,i)=nore(i,k) 
            ENDDO
         ENDDO
! copy into output, if output array available
         unc%g(1:nd,1:nd)=cove(1:nd,1:nd)
         unc%c(1:nd,1:nd)=nore(1:nd,1:nd)
      ENDIF
   ENDIF
   IF(.NOT.ngr_opt)THEN
      dyn=dyn_loc
   ELSE
   END IF
   IF(.not.PRESENT(unit)) CALL clorbf
! succesful termination  
   RETURN
! error condition, reading not succeded
20 CONTINUE 
   WRITE(*,*)'read_elems: error condition '
   err=.true.  
 END SUBROUTINE read_elems


! ==================================================================
! WRITE_ELEMS
! Writes elements and dynamical parameters, if present, by using 
! only information from dummy arguments. The size of the covar/normal
! matrices is determined by unc%ndim.
! A record contains the list ls, the integer nls.
! ==================================================================
SUBROUTINE write_elems(el,name,form,dynam_par,nlsloc,lsloc,unc,incfit,file,unit,unit2)
  USE station_coordinates, ONLY: codestat
  USE io_elems, ONLY: obscod
! ====BEGIN INTERFACE ====================================================================================
  TYPE(orbit_elem),          INTENT(IN) :: el            ! orbital elements
  CHARACTER*(*),             INTENT(IN) :: name          ! name
  CHARACTER*(*),   OPTIONAL, INTENT(IN) :: form          ! can be either '1L' or 'ML',
  TYPE(dyn_par),   OPTIONAL, INTENT(IN) :: dynam_par     ! dynamical non-grav parameters
  INTEGER,         OPTIONAL, INTENT(IN) :: nlsloc        ! number of effective solve-for par
  INTEGER,         OPTIONAL, INTENT(IN) :: lsloc(ndyx)   ! solve-for dynamic parameters
  TYPE(orb_uncert),OPTIONAL, INTENT(IN) :: unc           ! covariance/normal matrices 
  INTEGER,         OPTIONAL, INTENT(IN) :: incfit        ! result from incomplete fit;
  CHARACTER*(*),   OPTIONAL, INTENT(IN) :: file          ! output file 
                                                         ! (not used if PRESENT(unit))
  INTEGER,         OPTIONAL, INTENT(IN) :: unit          ! input unit, only if already opened
                                                         ! WARNING: if file already opened, 
                                                         ! it must also have the header already written
  INTEGER,         OPTIONAL, INTENT(IN) :: unit2         ! input unit, only if already opened for 1L catalog
                                                         ! in case of non-grav
! ====END INTERFACE =======================================================================================
  INTEGER           :: nd                                         ! dimension of covar/normal matrices 
  INTEGER           :: uniout, uniout_ngr                         ! output orbit file and non-grav file
  CHARACTER*10      :: rsys,epoch
  INTEGER           :: le, i, k                                   ! loop indices
! Local version of the optional variables
  CHARACTER*2       :: formcur
  CHARACTER*50      :: namloc
  TYPE(dyn_par)     :: dyn_loc
  INTEGER           :: sum_active                                 ! sum of the active non-grav parameters
  CHARACTER*6       :: vers                                       ! version
  DOUBLE PRECISION  :: cnv(ndimx),std(ndimx), ele(6), princ, mass
  DOUBLE PRECISION  :: units(ndimx)                               ! for scaling
! eigenvalues, eigenvectors                                             
  DOUBLE PRECISION  :: eigvec(ndimx,ndimx),eigval(ndimx),fv1(ndimx),fv2(ndimx),wdir(ndimx),sdir
  INTEGER           :: ierr 
  INCLUDE 'parcmc.h90'
! ==========================================================================================================
! initialization
  sum_active=0
! form
  IF(PRESENT(form))THEN
     formcur=form
  ELSE
     formcur='ML'
  ENDIF
! version
  vers='OEF2.0'
! dynamical parameters
  IF(PRESENT(dynam_par))THEN
     dyn_loc=dynam_par
  ELSE
     dyn_loc=undefined_dyn_par
  END IF
! dyn%nmod must be ge 0
  IF(dyn_loc%nmod.lt.0)THEN
     WRITE(*,*)'write_elems: wrong nongrav model nmod=', dyn_loc%nmod
     STOP
  END IF
! ndyx is the maximum for ndp
  IF(dyn_loc%ndp.gt.ndyx)THEN
     WRITE(*,*)'write_elems: wrong dimension given for dp, ndp= ', dyn_loc%ndp
     STOP
  ENDIF
! lsloc and nlsloc must come together, and with meaningful dimensions
  IF(PRESENT(lsloc).AND..NOT.PRESENT(nlsloc))THEN
     WRITE(*,*)'write_elems: lsloc is present without nlsloc ', lsloc
     STOP
  END IF
  IF(PRESENT(lsloc).AND.nlsloc.GT.0)THEN
     IF(nlsloc.GT.ndyx)THEN
        WRITE(*,*)'write_elems: wrong dimension for lsloc, nlsloc= ',nlsloc
        STOP
     ELSEIF(dyn_loc%ndp.LT.nlsloc)THEN
        WRITE(*,*)'write_elems: matrix larger than state vector ',nlsloc, dyn_loc%ndp
        STOP
     ENDIF
  ELSE IF(.NOT.PRESENT(lsloc))THEN
     IF(PRESENT(nlsloc).AND.nlsloc.gt.0)THEN
        WRITE(*,*)'write_elems: dimension given for non existing lsloc, nlsloc=',nlsloc
        STOP
     ENDIF
  ENDIF
  IF(.NOT.PRESENT(unit))THEN
     IF(PRESENT(file))THEN
        CALL filopn(uniout,file,'unknown')
! write file header
        rsys='ECLM'
        epoch='J2000'
        IF(formcur.eq.'1L')THEN
           CALL wro1lh2(uniout,rsys,epoch,el%coo,vers)
        ELSEIF(formcur.eq.'ML')THEN
           CALL wromlh2(uniout,rsys,epoch,vers)
        ENDIF
     ENDIF
  ELSE
     uniout=unit
  ENDIF
  IF(PRESENT(unit2))THEN
     uniout_ngr=unit2
  ELSE
     IF(PRESENT(file).AND.dyn_loc%ndp.GT.0)THEN
        IF(formcur.EQ.'1L')THEN
           CALL filopn(uniout_ngr,file,'unknown')
        ! write file header
           rsys='ECLM'
           epoch='J2000'
           CALL wro1lh2_ngr(uniout_ngr,rsys,epoch,vers,dyn_loc%ndp)
        ENDIF
     END IF
  END IF
! name without blanks
  namloc=' '
  namloc=name
  CALL rmsp(namloc,le)
! solutions with less than 6 parameters
  IF(PRESENT(incfit))THEN
     IF(incfit.eq.5)THEN
        WRITE(uniout,200)namloc
200     FORMAT('! CONSTRAINED SOLUTION for ',a)
     ELSEIF(incfit.eq.4)THEN
        WRITE(uniout,291)name
291     FORMAT('! 4_PARAMETERS SOLUTION for ',a)
     ENDIF
  ENDIF

! covariance/normal matrices
  IF(PRESENT(nlsloc).AND.nlsloc.gt.0)THEN
     nd=6+nlsloc
  ELSE
     nd=6
  ENDIF

  IF(PRESENT(unc))THEN
     IF(unc%ndim.ne.nd)THEN
        WRITE(*,*)' write_elems: ndim.ne.nd ', unc%ndim, nd
        !STOP
     ENDIF
  ENDIF

! when el%coo=ATT, station numeric code is needed
  CALL codestat(el%obscode,obscod)

! write in format 1L/ML
  IF(formcur.eq.'1L')THEN
! Check active non-grav parameters
     IF(dyn_loc%nmod.GT.0)THEN
        DO i=1,nlsloc
           IF(dyn_loc%active(lsloc(i)))THEN
              sum_active=sum_active+1
           END IF
        END DO
     END IF
     IF(sum_active.GT.0)THEN
        CALL wro1lr(uniout,namloc(1:le),el%coord,el%coo,el%t,el%h_mag,el%g_mag,sum_active)
        CALL wro1lr_ngr(uniout_ngr,namloc(1:le),dyn_loc%nmod,dyn_loc%ndp,dyn_loc%dp)
     ELSE
        CALL wro1lr(uniout,namloc(1:le),el%coord,el%coo,el%t,el%h_mag,el%g_mag)
     END IF
  ELSEIF(formcur.eq.'ML')THEN
! object name 
     WRITE(uniout,'(A)') namloc(1:le) 
! Orbital elements                                                      
     cnv=1.d0 
     ele=el%coord
     IF(el%coo.EQ.'KEP') THEN 
        cnv(3:6)=degrad 
        ele(6)=princ(ele(6))
        IF(ele(5).lt.0.d0)ele(5)=ele(5)+dpig 
        IF(ele(4).lt.0.d0)ele(4)=ele(4)+dpig 
        WRITE(uniout,201) comcha 
201     FORMAT(A,' Keplerian elements: a, e, i, long. node, arg. peric., mean anomaly')  
        WRITE(uniout,101) (ele(i)*cnv(i),i=1,6) 
101     FORMAT(' KEP ',1P,E24.16,0P,F18.15,4F18.13) 
     ELSEIF(el%coo.EQ.'CAR') THEN 
        WRITE(uniout,202) comcha 
202     FORMAT(A,' Cartesian position and velocity vectors') 
        WRITE(uniout,102) ele 
102     FORMAT(' CAR ',1P,6E22.14) 
     ELSEIF(el%coo.EQ.'EQU') THEN 
        cnv(6)=degrad 
        ele(6)=princ(ele(6))
        WRITE(uniout,203) comcha 
203     FORMAT(A,' Equinoctial elements: a, e*sin(LP), e*cos(LP),',       &
             &         ' tan(i/2)*sin(LN), tan(i/2)*cos(LN), mean long.') 
        WRITE(uniout,103) (ele(i)*cnv(i),i=1,6) 
103     FORMAT(' EQU ',1P,E24.16,0P,2(1x,f19.15),2(1x,f20.15),1x,F17.13) 
     ELSEIF(el%coo.EQ.'COM') THEN 
        cnv(3:5)=degrad 
        if(ele(5).lt.0.d0)ele(5)=ele(5)+dpig 
        if(ele(4).lt.0.d0)ele(4)=ele(4)+dpig 
        WRITE(uniout,204) comcha 
204     FORMAT(A,' Cometary elements: q, e, i, long. node,', &
             &         ' arg. peric., pericenter time')  
        IF(abs(ele(6)).lt.1.d6)THEN
           WRITE(uniout,114) (ele(i)*cnv(i),i=1,6) 
114        FORMAT(' COM ',1P,E22.14,0P,1x,F18.15,1x,3(F18.13,1x),F18.10)
        ELSE
           WRITE(uniout,215) (ele(i)*cnv(i),i=1,6) 
215        FORMAT(' COM ',1P,E22.14,0P,1x,F18.15,1x,3(F18.13,1x),1P,D18.11)
        ENDIF
     ELSEIF(el%coo.EQ.'COT') THEN 
        cnv(3:6)=degrad 
        if(ele(5).lt.0.d0)ele(5)=ele(5)+dpig 
        if(ele(4).lt.0.d0)ele(4)=ele(4)+dpig
        if(ele(6).lt.0.d0)ele(6)=ele(6)+dpig
        WRITE(uniout,206) comcha 
206     FORMAT(A,' Cometary elements: q, e, i, long. node,', &
             &         ' arg. peric., true anomaly')  
        WRITE(uniout,117) (ele(i)*cnv(i),i=1,6) 
117     FORMAT(' COT ',1P,E22.14,0P,1x,F18.15,1x,4(F18.13,1x))
     ELSEIF(el%coo.eq.'ATT') THEN
        cnv(1:4)=degrad
        IF(ele(1).lt.0.d0)ele(1)=ele(1)+dpig
        WRITE(uniout,205) comcha 
205     FORMAT(A,' Attributable elements: R.A., DEC, R.A.dot, DECdot, r, rdot,', &
             &         ' auxiliary information on observer needed')  
        WRITE(uniout,115) (ele(i)*cnv(i),i=1,6) 
115     FORMAT(' ATT ',2(F18.13,1x),2(F18.13,1x),2(F18.13,1x)) 
        WRITE(uniout,116) ' STA ',obscod
116     FORMAT(A5,A3)
     ELSE
        WRITE(*,*)  '**** write_elems: unsupported orbital element type ****', el%coo
        STOP '**** write_elems: unsupported orbital element type ****' 
     END IF
! Epoch                                                                 
     WRITE(uniout,104) el%t 
104  FORMAT(' MJD ',F19.9,' TDT')                                                                     
! Mass: not ready                                                                  
!     IF(mass.NE.0.d0) WRITE(uniout,105) mass 
!105  FORMAT(' MAS ',1P,E20.12)                                              
! Magnitudes                                                            
     !IF(el%mag_set) WRITE(uniout,106)el%h_mag,el%g_mag 
     IF(el%h_mag.gt.-100.d0) WRITE(uniout,106)el%h_mag,el%g_mag 
106  FORMAT(' MAG ',2F7.3)                                                 

     !IF(formcur.EQ.'ML')THEN
! nongrav parameters
     IF(dyn_loc%nmod.ge.0)THEN
! ndp.gt.0 implies PRESENT(nmod) nmod.gt.0
! record giving the matrix dimension and list of solved parameters
        IF(PRESENT(nlsloc).AND.nlsloc.GT.0)THEN
           WRITE(uniout,149) comcha
149        FORMAT(A,' Non-grav parameters: model used, actual number in use, dimension, ', &
                & 'list of solve for parameters to be currently determined')
           WRITE(uniout,150) dyn_loc%nmod, dyn_loc%ndp, nd, lsloc(1:nlsloc)
150        FORMAT(' LSP ',1X,I2,1X,I2,3X,I2,2X,4(1X,I2))
        ELSEIF(PRESENT(incfit).AND.incfit.EQ.5) THEN
           WRITE(uniout,'(A)') '! Non-grav parameters: model used, actual number in use, dimension'             
           WRITE(uniout,155) dyn_loc%nmod, dyn_loc%ndp, 5
        ELSE
           WRITE(uniout,'(A)') '! Non-grav parameters: model used, actual number in use, dimension'             
           WRITE(uniout,155) dyn_loc%nmod, dyn_loc%ndp, nd
155        FORMAT(' LSP ',1X,I2,1X,I2,3X,I2)
        END IF
        IF(dyn_loc%nmod.gt.0)THEN
           ! because ndp.gt.nlsloc, weights are available for output matrices
           cnv(7:6+dyn_loc%ndp)=1.d0
           WRITE(uniout,'(A)') '! Vector of dynamical parameters'
           WRITE(uniout,112) ((dyn_loc%dp(i))*cnv(6+i),i=1,dyn_loc%ndp)
112        FORMAT(' NGR ',1P,4E22.14)
        END IF
     ENDIF
     !END IF

! A-priori for non-gravitational parameter     
! ========================================
! Dynamical model nmod=1
     IF(dyn_loc%nmod.EQ.1)THEN
! A/M and A2 to be determined
        IF(lsloc(1).EQ.1.AND.lsloc(2).EQ.2)THEN
! A/M and A2 a-priori available
           IF(dyn_loc%apriori(2).AND.dyn_loc%apriori(1))THEN
              WRITE(uniout,'(A)') '! A priori value and standard deviation (A/M and Yarkovsky effect)'
              WRITE(uniout,352) dyn_loc%apriori_nominal(1), dyn_loc%apriori_std(1)
              WRITE(uniout,351) dyn_loc%apriori_nominal(2), dyn_loc%apriori_std(2)
! only A2 a-priori available
           ELSE IF(dyn_loc%apriori(2))THEN
              WRITE(uniout,'(A)') '! A priori value and standard deviation (only Yarkovsky effect)'
              WRITE(uniout,351) dyn_loc%apriori_nominal(2), dyn_loc%apriori_std(2)
! only A/M a-priori available  
           ELSE IF(dyn_loc%apriori(1))THEN
              WRITE(uniout,'(A)') '! A priori value and standard deviation (only A/M)'
              WRITE(uniout,352) dyn_loc%apriori_nominal(1), dyn_loc%apriori_std(1)
           END IF
! only A2 to be determined with a-priori
        ELSE IF(lsloc(1).EQ.2.AND.dyn_loc%apriori(2))THEN 
           WRITE(uniout,'(A)') '! A priori value and standard deviation (only Yarkovsky effect)'
           WRITE(uniout,351) dyn_loc%apriori_nominal(2), dyn_loc%apriori_std(2)
! only A/M to be determined with a-priori
        ELSE IF(lsloc(1).EQ.1.AND.dyn_loc%apriori(1))THEN    
           WRITE(uniout,'(A)') '! A priori value and standard deviation (only A/M)'
           WRITE(uniout,352) dyn_loc%apriori_nominal(1), dyn_loc%apriori_std(1)
        END IF

! Dynamical model nmod=2
     ELSEIF(dyn_loc%nmod.EQ.2)THEN
        IF(lsloc(1).EQ.1.AND.lsloc(2).EQ.2.AND.lsloc(3).EQ.3.AND.lsloc(4).EQ.4)THEN
           IF(dyn_loc%apriori(1).AND.dyn_loc%apriori(2).AND.dyn_loc%apriori(3).AND.dyn_loc%apriori(4))THEN
              WRITE(uniout,'(A)') '! A priori value and standard deviation (A1, A2, A3 and dtdelay)'
              WRITE(uniout,353) dyn_loc%apriori_nominal(1), dyn_loc%apriori_std(1)
              WRITE(uniout,354) dyn_loc%apriori_nominal(2), dyn_loc%apriori_std(2)
              WRITE(uniout,355) dyn_loc%apriori_nominal(3), dyn_loc%apriori_std(3)
              WRITE(uniout,356) dyn_loc%apriori_nominal(4), dyn_loc%apriori_std(4)
           END IF
        ELSE IF(lsloc(1).EQ.1.AND.lsloc(2).EQ.2.AND.lsloc(3).EQ.3)THEN
           IF(dyn_loc%apriori(1).AND.dyn_loc%apriori(2).AND.dyn_loc%apriori(3))THEN
              WRITE(uniout,'(A)') '! A priori value and standard deviation (A1, A2, A3)'
              WRITE(uniout,353) dyn_loc%apriori_nominal(1), dyn_loc%apriori_std(1)
              WRITE(uniout,354) dyn_loc%apriori_nominal(2), dyn_loc%apriori_std(2)
              WRITE(uniout,355) dyn_loc%apriori_nominal(3), dyn_loc%apriori_std(3)
           END IF
        ELSE IF(lsloc(1).EQ.1.AND.lsloc(2).EQ.2)THEN
           IF(dyn_loc%apriori(1).AND.dyn_loc%apriori(2))THEN
              WRITE(uniout,'(A)') '! A priori value and standard deviation (A1, A2)'
              WRITE(uniout,353) dyn_loc%apriori_nominal(1), dyn_loc%apriori_std(1)
              WRITE(uniout,354) dyn_loc%apriori_nominal(2), dyn_loc%apriori_std(2)
           END IF
        END IF
     END IF
351  FORMAT(' APR YAR ',1P,2E22.14)          
352  FORMAT(' APR A/M ',1P,2E22.14)
353  FORMAT(' APR A1C ',1P,2E22.14)
354  FORMAT(' APR A2C ',1P,2E22.14)
355  FORMAT(' APR A3C ',1P,2E22.14)
356  FORMAT(' APR DTD ',1P,2E22.14)

! covariance matrix                                                     
     IF(PRESENT(unc)) THEN 
        IF(.not.unc%succ)THEN
           WRITE(*,*)' write_elems: covariance not given'
        ELSE
           DO i=1,nd
              std(i)=SQRT(unc%g(i,i))
           ENDDO
! eigenvalues
           CALL weak_dir(unc%g(1:nd,1:nd),wdir(1:nd),sdir,-1,el%coo,ele,units(1:nd),nd)
           CALL rs(nd,nd,unc%g(1:nd,1:nd),eigval(1:nd),1,eigvec(1:nd,1:nd),fv1,fv2,ierr) 
           DO i=1,nd 
              IF(eigval(i).gt.0.d0)THEN 
                 eigval(i)=sqrt(eigval(i)) 
              ELSE 
                 IF(eigval(i).lt.-1.d-10.and.verb_io.gt.9)THEN
                    WRITE(*,*)'wromlr: zero/negative eigenvalue', eigval(i),'  for asteroid ',namloc(1:le)
                 ENDIF
                 eigval(i)=-sqrt(-eigval(i)) 
              ENDIF
           ENDDO
! RMS, eigenvalues and weak direction are so far commented   
! format to be changed for ndimx>10           
           WRITE(uniout,107) comcha,(std(i)*cnv(i),i=1,nd) 
107        FORMAT(A1,' RMS ',1P,10E14.5) 
           WRITE(uniout,111) comcha,eigval(1:nd)
111        FORMAT(A1,' EIG',1P,10E14.5) 
           WRITE(uniout,110) comcha,wdir(1:nd)
110        FORMAT(A1,' WEA',10F10.5) 
! covariance matrix is given uncommented, to be readable                
           WRITE(uniout,108) ((unc%g(i,k)*cnv(i)*cnv(k),k=i,nd),i=1,nd) 
108        FORMAT(' COV ',1P,3E23.15)
           WRITE(uniout,109) ((unc%c(i,k)/(cnv(i)*cnv(k)),k=i,nd),i=1,nd) 
109        FORMAT(' NOR ',1P,3E23.15) 
        ENDIF
     END IF
  ENDIF
  IF(.not.PRESENT(unit)) CALL filclo(uniout,' ')
END SUBROUTINE write_elems

! Copyright (C) 1997-1998 by Mario Carpino (carpino@brera.mi.astro.it)  
! Version: June 21, 1998                                                
! Version December 18, 1998 Steven Chesley (chesley@dm.unipi.it)        
! Added ELTYPE, Changed time field to F13.4 to handle jpleph406 
! Version May 20, 2015 Federica Spoto (spoto@mail.dm.unipi.it)
! Now the 1L catalogs must contain a flag to indicate if the dynamical 
! parameters are active        
! --------------------------------------------------------------------- 
!                                                                       
!  *****************************************************************    
!  *                                                               *    
!  *                         W R O 1 L R                           *    
!  *                                                               *    
!  *  Writes an orbital element record in an orbital element file  *    
!  *   (single-line format, different epochs, keplerian elements)  *    
!  *                     (multi-line format)                       *    
!  *                                                               *    
!  *****************************************************************    
!                                                                       
! OUTPUT:   UNIT      -  Output FORTRAN unit                            
!           NAME      -  Name of planet/asteroid/comet                  
!           ELEM(6)   -  Orbital element vector (keplerian elements)    
!           ELTYPE    -  Type of orbital elements (KEP/EQU/CAR)         
!           T0        -  Epoch of orbital elements (MJD, TDT)           
!           H         -  H absolute magnitude (if <-100, missing)       
!           G         -  G slope parameter                              
!                                                                       
! WARNING: the routine does not write the header of the file: this      
!          must be generated by calling subroutine wro1lh               
!                                                                       
SUBROUTINE wro1lr(unit,name,elem,eltype,t0,h,g,ngr_act) 
  USE fund_const
  IMPLICIT NONE
! Begin interface
  INTEGER,           INTENT(IN)           :: unit 
  DOUBLE PRECISION,  INTENT(IN)           :: elem(6),t0,h,g 
  CHARACTER*(*),     INTENT(IN)           :: name,eltype 
  INTEGER,           INTENT(IN), OPTIONAL :: ngr_act
! End interface
! Expected max length of name                                           
  INTEGER, PARAMETER :: namtl=12 
  DOUBLE PRECISION   :: cnvele(6), princ                                         
  INTEGER            ::ln,nb,i 
  CHARACTER*(namtl)  :: blanks
  INTEGER            ::lench 
  EXTERNAL lench 
  CHARACTER*80       ::elefmt 
  INTEGER            ::ngr_act_loc
!============================================================
! Name                                                                  
  ln=MAX(1,lench(name)) 
  nb=MAX(1,namtl-ln) 
  blanks=' ' 
! Check if there are non-grav parameters
  IF(PRESENT(ngr_act))THEN
     ngr_act_loc=ngr_act
  ELSE
     ngr_act_loc=0
  END IF
! Convert to degrees 
  cnvele=elem 
  IF(eltype .eq. 'CAR')THEN 
     elefmt='('''''''',A,'''''''',A,1x,F13.6,'//                    &
     &        '6(1x,F17.13),'//                                     &
     &        '2F6.2,1X,I2)'                                                  
  ELSEIF(eltype .eq. 'KEP')THEN 
     cnvele(3:6)=cnvele(3:6)*degrad 
     elefmt='('''''''',A,'''''''',A,F13.6,'//                       &
     &        '1p,6(e25.16),'//                                     &
     &        '0p,2F6.2,1X,I2)'                                               
!        elefmt='('''''''',A,'''''''',A,F13.6,'//                       
!    +        'F16.12,1x,F12.10,4(1x,F12.7),'//                         
!    +        '2F6.2)' 
  ELSEIF(eltype .eq. 'COM')THEN 
     cnvele(3:5)=cnvele(3:5)*degrad 
     elefmt='('''''''',A,'''''''',A,F13.6,'//                       &
     &        '1p,5(e25.16),0p,f14.6,'//                            &
     &        '0p,2F6.2,1X,I2)'                             
  ELSEIF(eltype .eq. 'COT')THEN 
     cnvele(3:6)=cnvele(3:6)*degrad 
     elefmt='('''''''',A,'''''''',A,F13.6,'//                       &
     &        '1p,6(e25.16),0p,2F6.2,1X,I2)'
  ELSEIF(eltype .eq. 'EQU')THEN 
     cnvele(6)=princ(cnvele(6))
     cnvele(6)=cnvele(6)*degrad 
     elefmt='('''''''',A,'''''''',A,F13.6,'//                       &
     &        'F16.12,4(1x,f12.9),1x,F12.7,'//                      &
     &        '2F6.2,1X,I2)'                                                  
  ELSEIF(eltype.eq.'ATT')THEN
     WRITE(*,*)' wro1lr: warning, the ATT elelments in  this format have no station code '
     cnvele(1:4)=cnvele(1:4)*degrad
     elefmt='('''''''',A,'''''''',A,F13.6,'//                       &
     &        '1p,6(e25.16),'//                                     &
     &        '0p,2F6.2,1X,I2)'  
  ELSE 
     WRITE(*,*)'*** wro1lr: unsupported coordinate type', eltype
     STOP '*** wro1lr: unsupported coordinate type' 
  ENDIF
! Write                                                                 
  IF(h.GT.-100.d0) THEN 
     WRITE(unit,elefmt) name(1:ln),blanks(1:nb),t0,cnvele,h,g,ngr_act_loc
  ELSE 
     WRITE(unit,elefmt) name(1:ln),blanks(1:nb),t0,cnvele,ngr_act_loc
  ENDIF                                                                      
END SUBROUTINE wro1lr

SUBROUTINE wro1lr_matlab(unit,name,elem,eltype,t0,h,ngr_act,chi,rms,succ) 
  USE fund_const
  IMPLICIT NONE
! Begin interface
  INTEGER,           INTENT(IN)           :: unit 
  DOUBLE PRECISION,  INTENT(IN)           :: elem(6),t0,h
  CHARACTER*(*),     INTENT(IN)           :: name,eltype 
  INTEGER,           INTENT(IN), OPTIONAL :: ngr_act
  DOUBLE PRECISION,  INTENT(IN), OPTIONAL :: chi
  DOUBLE PRECISION,  INTENT(IN), OPTIONAL :: rms
  INTEGER,           INTENT(IN), OPTIONAL :: succ
! End interface
! Expected max length of name                                           
  INTEGER, PARAMETER :: namtl=12 
  DOUBLE PRECISION   :: cnvele(6), princ                                         
  INTEGER            :: ln,nb,i 
  CHARACTER*(namtl)  :: blanks
  INTEGER            :: lench 
  EXTERNAL lench 
  CHARACTER*80       :: elefmt 
  INTEGER            :: ngr_act_loc
  DOUBLE PRECISION   :: chi_loc
  DOUBLE PRECISION   :: rms_loc
  INTEGER            :: succ_loc
!============================================================
! Name                                                                  
  ln=MAX(1,lench(name)) 
  nb=MAX(1,namtl-ln) 
  blanks=' ' 
! Check if there are non-grav parameters
  IF(PRESENT(ngr_act))THEN
     ngr_act_loc=ngr_act
  ELSE
     ngr_act_loc=0
  END IF
  IF(PRESENT(chi))THEN
     chi_loc = chi
  ELSE
     chi_loc = 0.d0
  END IF
  IF(PRESENT(rms))THEN
     rms_loc = rms
  ELSE
     rms_loc = 0.d0
  END IF
  IF(PRESENT(succ))THEN
     succ_loc = succ
  END IF
! Convert to degrees 
  cnvele=elem 
  IF(eltype .eq. 'CAR')THEN 
     elefmt='(A,A,1x,F13.6,'//                    &
     &        '6(1x,F17.13),'//                                     &
     &        '2F6.2,1X,I2,1P,1X,D12.5,1X,D12.5)'                                                  
  ELSEIF(eltype .eq. 'KEP')THEN 
     cnvele(3:6)=cnvele(3:6)*degrad
     IF(PRESENT(succ)) THEN
        elefmt='(A,A,F13.6,'//                       &
             &        '1p,6(e25.16),'//                                     &
             &        '0p,F6.2,1X,I2,1P,1X,D12.5,1X,D12.5,1X,I1)'
     ELSE
        elefmt='(A,A,F13.6,'//                       &
             &        '1p,6(e25.16),'//                                     &
             &        '0p,F6.2,1X,I2,1P,1X,D12.5,1X,D12.5)'
     END IF
!        elefmt='('''''''',A,'''''''',A,F13.6,'//                       
!    +        'F16.12,1x,F12.10,4(1x,F12.7),'//                         
!    +        '2F6.2)' 
  ELSEIF(eltype .eq. 'COM')THEN 
     cnvele(3:5)=cnvele(3:5)*degrad 
     IF(PRESENT(succ))THEN
        elefmt='(A,A,F13.6,'//                       &
             &        '1p,5(e25.16),0p,f14.6,'//                            &
             &        '0p,F6.2,1X,I2,1P,1X,D12.5,1X,D12.5,1X,I1)'
     ELSE
        elefmt='(A,A,F13.6,'//                       &
             &        '1p,5(e25.16),0p,f14.6,'//                            &
             &        '0p,F6.2,1X,I2,1P,1X,D12.5,1X,D12.5,1X,I1)'
     END IF
  ELSEIF(eltype .eq. 'COT')THEN 
     cnvele(3:6)=cnvele(3:6)*degrad 
     elefmt='(A,A,F13.6,'//                       &
     &        '1p,6(e25.16),0p,2F6.2,1X,I2,1P,1X,D12.5,1X,D12.5)'
  ELSEIF(eltype .eq. 'EQU')THEN 
     cnvele(6)=princ(cnvele(6))
     cnvele(6)=cnvele(6)*degrad 
     elefmt='(A,A,F13.6,'//                       &
     &        'F16.12,4(1x,f12.9),1x,F12.7,'//                      &
     &        '2F6.2,1X,I2,1P,1X,D12.5,1X,D12.5)'                                                  
  ELSEIF(eltype.eq.'ATT')THEN
     WRITE(*,*)' wro1lr: warning, the ATT elelments in  this format have no station code '
     cnvele(1:4)=cnvele(1:4)*degrad
     elefmt='(A,A,F13.6,'//                       &
     &        '1p,6(e25.16),'//                                     &
     &        '0p,2F6.2,1X,I,1P,1X,D12.52,1X,D12.5)'  
  ELSE 
     WRITE(*,*)'*** wro1lr_matlab: unsupported coordinate type', eltype
     STOP '*** wro1lr_matlab: unsupported coordinate type' 
  ENDIF
! Write                                                                 
  IF(h.GT.-100.d0) THEN
     IF(PRESENT(succ))THEN 
        WRITE(unit,elefmt) name(1:ln),blanks(1:nb),t0,cnvele,h,ngr_act_loc, chi_loc, rms_loc, succ_loc
     ELSE
        WRITE(unit,elefmt) name(1:ln),blanks(1:nb),t0,cnvele,h,ngr_act_loc, chi_loc, rms_loc
     END IF
  ELSE 
     IF(PRESENT(succ))THEN
        WRITE(unit,elefmt) name(1:ln),blanks(1:nb),t0,cnvele,ngr_act_loc, chi_loc, rms_loc, succ_loc
     ELSE
        WRITE(unit,elefmt) name(1:ln),blanks(1:nb),t0,cnvele,ngr_act_loc, chi_loc, rms_loc
     END IF
  ENDIF                                                                      
END SUBROUTINE wro1lr_matlab

!  *****************************************************************    
!  *                                                               *    
!  *                         W R O 1 L H 2                         *    
!  *                                                               *    
!  *       Writes the header of an orbital element file            *    
!  *   (single-line format, different epochs, elements)            *    
!  *                                                               *    
!  *****************************************************************    
!                                                                       
! INPUT:    UNIT      -  Output FORTRAN unit                            
!           RSYS      -  Reference system type (EQUM/EQUT/ECLM)         
!           EPOCH     -  Reference system epoch (J2000/OFDATE)          
!           ELTYPE    -  Type of orbital elements (EQU/KEP/CAR)         
!           VERS      -  Version (old: OEF1.1 or new: OEF2.0 )
! 
                                                                      
SUBROUTINE wro1lh2(unit,rsys,epoch,eltype,version) 
  IMPLICIT NONE
! Begin interface
  INTEGER,      INTENT(IN)           :: unit 
  CHARACTER*(*),INTENT(IN)           :: rsys,epoch,eltype
  CHARACTER*6,  INTENT(IN), OPTIONAL :: version
! End interface
  INCLUDE 'parcmc.h90'                                                  
  INTEGER l1,l2,l3,nb 
  CHARACTER*100 bl
  CHARACTER*6 :: vers
  INTEGER lench 
  EXTERNAL lench                     
!========================================================
  ! Default is OEF2.0
  IF(PRESENT(version)) THEN
     vers=version
  ELSE
     vers='OEF2.0'
  END IF
  l1=lench(rsys) 
  l2=lench(epoch) 
  l3=lench(eltype) 
  nb=MAX(14-l1-l2,1) 
  bl=' '                                                                
  WRITE(unit,100)     vers(1:6),comcha,comcha,eltype(1:l3),comcha,      &
       &                rsys(1:l1),epoch(1:l2),                         &
       &                bl(1:nb),comcha                                   
100 FORMAT("format  = '",A,"'       ",A,' file format'/                     &
         &       'rectype = ''1L''           ',A,' record type (1L/ML)'/    &
         &       'elem    = ''',A,'''          ',A,                         &
         &                   ' type of orbital elements'/                   &
         &       'refsys  = ',A,1X,A,A,A,' default reference system'/       &
         &       'END_OF_HEADER')                                           
  
  IF(eltype.EQ.'KEP') THEN 
     WRITE(unit,201) comcha 
201  FORMAT(A,' Name, Epoch(MJD), a, e, i, long. node,',' arg. peric., mean anomaly, ' &
          & 'absolute magnitude, slope param, non-grav param') 
  ELSEIF(eltype.EQ.'CAR') THEN 
     WRITE(unit,202) comcha 
202  FORMAT(A,' Name, Epoch(MJD), cartesian position and velocity',' vectors, ' &
          & 'absolute magnitude, slope param, non-grav param')
  ELSEIF(eltype.EQ.'EQU') THEN 
     WRITE(unit,203) comcha 
203  FORMAT(A,' Name, Epoch(MJD), a, e*sin(LP), e*cos(LP),',' tan(i/2)*sin(LN), tan(i/2)*cos(LN), mean long., ' &
          & 'absolute magnitude, slope param, non-grav param')
  ELSEIF(eltype.EQ.'COM') THEN
     WRITE(unit,204) comcha 
204  FORMAT(A,' Name, Epoch(MJD), q, e, i, long. node,',' arg. peric., perihelion time, ' &
          & 'absolute magnitude, slope param, non-grav param')
  ELSEIF(eltype.EQ.'COT') THEN
     WRITE(unit,205) comcha 
205  FORMAT(A,' Name, Epoch(MJD), q, e, i, long. node,',' arg. peric., true anomaly, ' &
          & 'absolute magnitude, slope param, non-grav param')
  ELSEIF(eltype.eq.'ATT')THEN
     WRITE(unit,206) comcha 
206  FORMAT(A,' Name, Epoch(MJD), R.A., DEC, R.A.dot, DECdot, r, rdot, ' &
          & 'absolute magnitude, slope param, non-grav param')
  END IF
END SUBROUTINE wro1lh2

!  *****************************************************************    
!  *                                                               *    
!  *                   W R O 1 L H 2 _ N G R                       *    
!  *                                                               *    
!  *       Writes the header of a file containing non-grav param   *    
!  *              in case of a single line catalog                 *    
!  *                                                               *    
!  *****************************************************************    
!                                                                       
! INPUT:    UNIT      -  Output FORTRAN unit                            
!           RSYS      -  Reference system type (EQUM/EQUT/ECLM)         
!           EPOCH     -  Reference system epoch (J2000/OFDATE)          
!           VERS      -  Version ( OEF2.0 )
! 
                                                                      
SUBROUTINE wro1lh2_ngr(unit,rsys,epoch,vers,n_dp) 
  IMPLICIT NONE
! Begin interface
  INTEGER,      INTENT(IN)            :: unit 
  CHARACTER*(*),INTENT(IN)            :: rsys,epoch
  CHARACTER*6,  INTENT(IN)            :: vers
  INTEGER,      INTENT(IN), OPTIONAL  :: n_dp
! End interface
  INCLUDE 'parcmc.h90'                                                  
  INTEGER       :: l1,l2,nb 
  CHARACTER*100 :: bl
  INTEGER       :: lench 
  EXTERNAL lench                                                          
!======================================================
  l1=lench(rsys) 
  l2=lench(epoch) 
  nb=MAX(14-l1-l2,1) 
  bl=' '                                                                
  WRITE(unit,500)     vers(1:6),comcha,comcha,                              &
       &              rsys(1:l1),epoch(1:l2),                               &
       &              bl(1:nb),comcha                                   
500 FORMAT("format  = '",A,"'       ",A,' file format'/                     &
         &       'rectype = ''1L''           ',A,' record type (1L/ML)'/    &
         &       'refsys  = ',A,1X,A,A,A,' default reference system'/       &
         &       'END_OF_HEADER')                                           
  
  IF(PRESENT(n_dp)) THEN
     IF(n_dp.EQ.2) THEN
        WRITE(unit,501) comcha 
501     FORMAT(A,' Name, dynamical model, actual number of parameters in use, A/M, A2 ') 
     ELSE IF(n_dp.EQ.3) THEN
        WRITE(unit,502) comcha 
502     FORMAT(A,' Name, dynamical model, actual number of parameters in use, A1, A2, A3 ') 
     ELSE IF(n_dp.EQ.4) THEN
        WRITE(unit,503) comcha 
503     FORMAT(A,' Name, dynamical model, actual number of parameters in use, A1, A2, A3, delay ') 
     ELSE IF(n_dp.EQ.1) THEN
        WRITE(unit,504) comcha 
504     FORMAT(A,' Name, dynamical model, actual number of parameters in use, eta ') 
     ELSE IF(n_dp.NE.0) THEN
        WRITE(*,*) 'wro1lh2_ngr: n_dp, wrong value ', n_dp
        STOP
     END IF
  ELSE
     WRITE(unit,505) comcha
505     FORMAT(A,' Name, dynamical model, actual number of parameters in use, A/M, A2' &
             & '            (asteroid case)'/                                          &
             & '!                                                             A1, A2'  &
             & ', A3, delay (comet case)')
  END IF
END SUBROUTINE wro1lh2_ngr

                                      
! --------------------------------------------------------------------- 
!                                                                       
!  *****************************************************************    
!  *                                                               *    
!  *                         W R O M L H 2                           *    
!  *                                                               *    
!  *       Writes the header of an orbital element file            *    
!  *                     (multi-line format)                       *    
!  *                                                               *    
!  *****************************************************************    
!                                                                       
! OUTPUT:   UNIT      -  Output FORTRAN unit                            
!           RSYS      -  Reference system type (EQUM/EQUT/ECLM)         
!           EPOCH     -  Reference system epoch (J2000/OFDATE)          
!           VERS      -  Version (OEF2.0)
!                                                                       
SUBROUTINE wromlh2(unit,rsys,epoch,version) 
  IMPLICIT NONE
! Begin interface
  INTEGER,       INTENT(IN)           :: unit 
  CHARACTER*(*), INTENT(IN)           :: rsys,epoch
  CHARACTER*6,   INTENT(IN), OPTIONAL :: version
! End interface
  INCLUDE 'parcmc.h90'                                                  
  INTEGER       :: l1,l2,nb 
  CHARACTER*100 :: bl 
  CHARACTER*6   :: vers
  INTEGER       :: lench 
  EXTERNAL lench
!============================================================= 
  IF(PRESENT(version)) THEN
     vers=version
  ELSE
     vers='OEF2.0'
  END IF
  l1=lench(rsys) 
  l2=lench(epoch) 
  nb=MAX(14-l1-l2,1) 
  bl=' '                                                                       
  WRITE(unit,100) vers(1:6),comcha,comcha,rsys(1:l1),epoch(1:l2),       &
       &                bl(1:nb),comcha                                   
100 FORMAT("format  = '",A6,"'       ",A,' file format'/                    &
         &       'rectype = ''ML''           ',A,' record type (1L/ML)'/    &
         &       'refsys  = ',A,1X,A,A,A,' default reference system'/       &
         &       'END_OF_HEADER')   
END SUBROUTINE wromlh2


! Copyright (C) 1997-2003 by Mario Carpino (carpino@brera.mi.astro.it), 
!                            Andrea Milani (milani@dm.unipi.it),        
!                            Zoran Knezevic (zoran@aob.aob.bg.ac.yu)    
! Version: 12 August 2003                                            
! --------------------------------------------------------------------- 
!                                                                       
!  *****************************************************************    
!  *                                                               *    
!  *                         R D E L E M                           *    
!  *                                                               *    
!  *          Read orbital elements for a list of objects          *    
!  *                     from a list of files                      *    
!  *                                                               *    
!  *****************************************************************    
!                                                                       
! INPUT:    UNIT      -  Output FORTRAN unit (report file)              
!           OBJNAM    -  Object names                                   
!           NOBJ      -  Number of objects                              
!           INFILES   -  Input files                                    
!           NFIL      -  Number of input files                          
!           DEFORB    -  Is the orbit already defined?                  
!                                                                       
! OUTPUT:   DEFORB    -  Was the orbit found?                           
!           DEFCN     -  Tells whether covariance/normal matrices       
!                            are defined                                
!           ELTYPE    -  Type of orbital elements (EQU/KEP/CAR)         
!           TELEM     -  Epoch of orbital elements (MJD, TDT)           
!           ELEM      -  Orbital elements (ECLM J2000)                  
!           COVE      -  Covariance matrix of orbital elements          
!           NORE      -  Normal matrix of orbital elements              
!           MASS      -  Mass (solar masses)                            
!           HMAG      -  H absolute magnitude (if <-100, missing)       
!           GMAG      -  G slope parameter                              
!           COMELE    -  Comment on orbital elements   
!                                                                       
! WARNING: the routine is designed to allow multiple calls on           
!          different lists of files/objects: for this reason,           
!          orbital elements of objects which have already an            
!          orbit defined (DEFORB=.true.) ARE NOT modified.              
!          For this reason, the user must define suitably DEFORB        
!          BEFORE calling RDELEM                                        
!                                                                       
SUBROUTINE rdelem(unit,objnam,nobj,infiles,nfil,deforb,defcn,      &
     &                  eltype,telem,elem,cove,nore,              &
     &                  mass,h,g,comele)  

  INTEGER nfilx 
  PARAMETER(nfilx=20) 
  INTEGER,          INTENT(IN)    :: unit,nobj,nfil 
  CHARACTER(LEN=*), INTENT(IN)    :: objnam(nobj),infiles(nfil)
  LOGICAL,          INTENT(INOUT) :: deforb(nobj)
  LOGICAL,          INTENT(OUT)   :: defcn(nobj)
  CHARACTER(LEN=*), INTENT(OUT)   :: eltype(nobj),comele(nobj) 
  DOUBLE PRECISION, INTENT(OUT)   :: telem(nobj),elem(6,nobj) 
  DOUBLE PRECISION, INTENT(OUT)   :: cove(ndimx,ndimx,nobj), nore(ndimx,ndimx,nobj)
  DOUBLE PRECISION, INTENT(OUT)   :: mass(nobj),h(nobj),g(nobj) 
!*********************************************************************************
  INTEGER   :: i,lf,is1,lfo,uniin,n 
  CHARACTER :: form*10,infil1*100
  LOGICAL   :: opened 
  INTEGER   :: lench 
  EXTERNAL lench 
                                                                        
! Number of input files                                                 
  IF(nfil.GT.nfilx) STOP '**** rdelem: nfil > nfilx ****' 
                                                                        
! Nothing to do                                                         
  IF(nfil.LE.0) RETURN 
            
! Loop on files                                                         
  DO 2 i=1,nfil 
     lf=lench(infiles(i)) 
     opened=.false. 
                                                                        
! Understand file format, if given in the argument infils(i)
     is1=index(infiles(i)(1:lf),'[') 
     IF(is1.GT.0 .AND. infiles(i)(lf:lf).EQ.']') THEN 
        form=infiles(i)(is1+1:lf-1) 
        lf=is1-1 
        infil1=infiles(i)(1:lf) 
     ELSE 
        infil1=infiles(i)
        form=' ' 
        CALL filopn(uniin,infil1(1:lf),'OLD') 
        opened=.true. 
        CALL oefdet(uniin,infil1(1:lf),form) ! SELF-DETECT FILE FORMAT 
        REWIND(uniin) 
     END IF
     lfo=lench(form) 
     IF(form.NE.' ') WRITE(*,102) infil1(1:lf),form(1:lfo) 
102  FORMAT('Scanning file "',A,'" (format: ',A,')') 
     
! Reading file                                                          
     IF(form.EQ.'OEF') THEN 
        IF(.NOT.opened) CALL filopn(uniin,infil1,'OLD')
        CALL rdoef(uniin,infil1,objnam,nobj,deforb,defcn,eltype,telem,&
             &     elem,cove,nore,mass,h,g,comele)
        CALL clorbf
     ELSEIF(form.EQ.'BA1') THEN 
        WRITE(*,*)' rdelem: obsolete Lowell format pre-1999'
        STOP
!        IF(.NOT.opened) CALL filopn(uniin,infil1,'OLD') 
!        CALL rdast1(uniin,infil1,objnam,nobj,deforb,defcn,            &
!     &                eltype,telem,elem,cove,nore,mass,h,g,comele)      
!        CALL filclo(uniin,' ') 
     ELSEIF(form.EQ.'BA2') THEN 
        IF(.NOT.opened) CALL filopn(uniin,infil1,'OLD') 
        CALL rdast2(uniin,infil1,objnam,nobj,deforb,defcn,            &
             &                eltype,telem,elem,cove,nore,mass,h,g,comele)      
        CALL filclo(uniin,' ') 
     ELSEIF(form.EQ.'MPC-A') THEN 
        IF(.NOT.opened) CALL filopn(uniin,infil1,'OLD') 
        CALL rdmpca(uniin,infil1,objnam,nobj,deforb,defcn,            &
             &                eltype,telem,elem,cove,nore,mass,h,g,comele)      
        CALL filclo(uniin,' ') 
     ELSEIF(form.EQ.' ') THEN 
        WRITE(*,120) infil1(1:lf) 
        STOP '**** rdelem: abnormal end ****' 
     ELSE 
        WRITE(*,121) form(1:lfo),infil1(1:lf) 
        STOP '**** rdelem: abnormal end ****' 
     END IF
120  FORMAT('ERROR: unknown format type for file "',A,'"') 
121  FORMAT('ERROR: unsupported format "',A,'" for file "',A,'"') 
! set default for magnitude data                                        
     DO n=1,nobj 
        IF(deforb(n))THEN
           IF(h(n).lt.-1.d6)THEN 
! leave it                                                              
           ENDIF
           IF(g(n).lt.-1.d6)THEN 
              g(n)=0.15d0 
           ENDIF
        ENDIF
     ENDDO
! End of loop on files                                                  
2 END DO

END SUBROUTINE rdelem

! Copyright (C) 1998 by Mario Carpino (carpino@brera.mi.astro.it),      
! Version: June 19, 1998
! --------------------------------------------------------------------- 
!                                                                       
!  *****************************************************************    
!  *                                                               *    
!  *                          R D O E F                            *    
!  *                                                               *    
!  *          Read orbital elements for a list of objects          *    
!  *         from a file written in internal ORBFIT format         *    
!  *                                                               *    
!  *****************************************************************    
!                                                                       
! INPUT:    UNIIN     -  Input unit (file is already opened)            
!           FILE      -  Input file name (for error messages only)      
!           OBJNAM    -  Object names (without embedded blanks)         
!           NOBJ      -  Number of objects                              
!                                                                       
! OUTPUT:   DEFORB    -  Tells whether orbital elements are defined     
!           DEFCN     -  Tells whether covariance/normal matrices       
!                            are defined                                
!           ELTYPE    -  Type of orbital elements (EQU/KEP/CAR/COM?ATT)         
!           TELEM     -  Epoch of orbital elements (MJD, TDT)           
!           ELEM      -  Orbital elements (ECLM J2000)                  
!           COVE      -  Covariance matrix of orbital elements          
!           NORE      -  Normal matrix of orbital elements              
!           MASS      -  Mass (solar masses)                            
!           H         -  H absolute magnitude (if <-100, missing)       
!           G         -  G slope parameter                              
!           COMELE    -  Comment on orbital elements                    
!                                                                       
! WARNING: the routine assumes that objects having DEFORB=.true.        
!          have already orbital elements defined (possibly from another 
!          input file) and does not overwrite them                      
!                                                                       
SUBROUTINE rdoef(uniin,file,objnam,nobj,deforb,defcn,eltype,telem,&
     &                 elem,cove,nore,mass,h,g,comele)                  
  USE fund_const
  USE reference_systems
  USE name_rules
  IMPLICIT NONE
  ! BEGIN INTERFACE
  INTEGER, INTENT(IN) :: uniin 
  INTEGER, INTENT(IN) :: nobj 
  CHARACTER*(*), INTENT(IN) :: file,objnam(nobj)
  DOUBLE PRECISION, INTENT(OUT) :: telem(nobj),elem(6,nobj) 
  DOUBLE PRECISION, INTENT(OUT) :: cove(ndimx,ndimx,nobj), nore(ndimx,ndimx,nobj)
  DOUBLE PRECISION, INTENT(OUT) :: mass(nobj),h(nobj),g(nobj) 
  CHARACTER*(*), INTENT(OUT) :: eltype(nobj),comele(nobj) 
  LOGICAL, INTENT(INOUT) :: deforb(nobj),defcn(nobj)
  ! END INTERFACE
 
  CHARACTER*(20) filedum 
  INTEGER kr,k,j1,j2,lf,ln,lc,lc1,nrem 
  DOUBLE PRECISION t1,gmsun,gma,gma1,enne,h1,g1,m1 
  DOUBLE PRECISION rot(3,3) 
  TYPE(orbit_elem) :: elem_loc
  TYPE(orb_uncert) :: unc, unc2
  DOUBLE PRECISION elem1(6),xv(6),cove1(ndimx,ndimx),cove2(ndimx,ndimx) 
  DOUBLE PRECISION de(6,6),nore1(ndimx,ndimx),nore2(ndimx,ndimx) 
  CHARACTER eltyp1*3,rsys*10,epoch*10,krc*10 
  LOGICAL defcov,defnor,end,err 
  CHARACTER*(idnamvir_len) name1,nc1
  INTEGER lench 
  INTEGER :: nlsloc,lsloc(ndyx)
  !Dimensions
  INTEGER :: nd
  EXTERNAL lench 

  gmsun=gms 
                                                                        
! Number of remaining object (orbit not yet found)                      
  nrem=0 
  DO 1 k=1,nobj 
     IF(.NOT.deforb(k)) nrem=nrem+1 
1 END DO
  IF(nrem.LE.0) RETURN 
                                                            
  CALL oporbf(file,uniin) 
  lf=lench(file)
                                                                        
3 CONTINUE 
!  CALL rdorb(name1,elem1,eltyp1,t1,cove1,defcov,nore1,defnor,       &
!     &           h1,g1,m1,rsys,epoch,kr,end)                            
  CALL read_elems(elem_loc,name1,end,nlsloc,err,lsloc,UNC=unc,FILE=file(1:lf),UNIT=uniin,MASS=m1,&
       & RSYS=rsys, EPOCH=epoch)
! Save dyn%dp, nls, and ls
  IF(dyn%ndp.GT.0.AND.nlsloc.NE.0.AND..NOT.ngr_opt)THEN
     nls=nlsloc
     ls=lsloc
  ENDIF
! Set dimension
  nd=6+nls
  
!  IF(end) GOTO 20 
! Name match is performed disregarding embedded blanks                  
  nc1=name1 
  CALL rmsp(nc1,lc1) 
  IF(lc1.LE.0) GOTO 3 
  DO 4 k=1,nobj
     IF(deforb(k)) GOTO 4 
     IF(nc1(1:lc1).EQ.objnam(k)) THEN 
        IF(rsys.EQ.'ECLM' .AND. epoch.EQ.'J2000') THEN
           eltype(k)=elem_loc%coo
           DO j1=1,nd
              DO  j2=1,nd 
                 cove(j1,j2,k)=unc%g(j1,j2) 
                 nore(j1,j2,k)=unc%c(j1,j2) 
              ENDDO
           ENDDO
           DO j1=1,6
              elem(j1,k)=elem_loc%coord(j1) 
           END DO
        ELSE 
           gma=gmsun*m1 
           gma1=gma+gmsun 
           IF(unc%succ) THEN 
! Transformation in cartesian coordinates                               
              CALL cooder(elem_loc%coord,elem_loc%coo,gma1,xv,'CAR',enne,de) 
              !CALL convertcovdp(6,cove1,de,cove2)
              !CALL norprsdp(nore1,de,6,nore2,error) 
              CALL convertunc(unc,de,unc2)
              !IF(error) THEN 
              !   ln=lench(name1) 
              !   WRITE(*,120) file(1:lf),name1(1:ln) 
              !   defnor=.false. 
              !END IF
! Transformation of reference system                                    
              CALL rotpn(rot,rsys,epoch,elem_loc%t,'ECLM','J2000',0.d0) 
              CALL prodmv(elem(1,k),rot,xv(1)) 
              CALL prodmv(elem(4,k),rot,xv(4)) 
              DO j1=1,3 
                 DO j2=1,3 
                    de(j1,j2)=rot(j1,j2) 
                    de(j1+3,j2)=0 
                    de(j1,j2+3)=0 
                    de(j1+3,j2+3)=rot(j1,j2) 
                 ENDDO
              ENDDO
               !IF(defcov) CALL convertcovdp(6,cove2,de,cove(1:6,1:6,k))
               !IF(defnor) THEN 
               !   CALL norprsdp(nore2,de,6,nore(1,1,k),error) 
               !   IF(error) THEN 
               !      ln=lench(name1) 
               !      WRITE(*,120) file(1:lf),name1(1:ln) 
               !      defnor=.false. 
               !   END IF
               !END IF
              IF(unc%succ)THEN
                 CALL convertunc(unc2,de,unc)
              END IF
           ELSE 
              CALL coocha(elem_loc,eltyp1,gma1,xv,'CAR',enne) 
              CALL rotpn(rot,rsys,epoch,t1,'ECLM','J2000',0.d0) 
              CALL prodmv(elem(1,k),rot,xv(1)) 
              CALL prodmv(elem(4,k),rot,xv(4)) 
           END IF
           eltype(k)='CAR' 
        END IF
        deforb(k)=.true.
        telem(k)=elem_loc%t 
        IF(unc%succ)THEN
           defcov=.TRUE.
           defnor=.TRUE.
        END IF
        CALL fixcnm(defcov,defnor,defcn(k),cove(1,1,k),nore(1,1,k)) 
        mass(k)=m1 
        h(k)=elem_loc%h_mag
        g(k)=elem_loc%g_mag
!        WRITE(krc,101) kr 
 !       CALL rmsp(krc,lc) 
        comele(k)='read from file "'//file(1:lf)//'"'
! &              '" at record '//krc(1:lc)                           
        nrem=nrem-1 
        IF(nrem.LE.0) GOTO 20 
     END IF
101  FORMAT(I6) 
120  FORMAT(' rdoef: error in transforming normal matrix'/             &
          &       '        (file "',A,'", object "',A,'")')                  
4 END DO
  GOTO 3 
20 CONTINUE 
END SUBROUTINE rdoef


END MODULE orbit_elements

! ===================================================================== 
!                                                                       
!  ***************************************************************      
!  *                                                             *      
!  *                         A P P M A G                         *      
!  *                                                             *      
!  *     Calcolo della magnitudine apparente di un asteroide     *      
!  *                                                             *      
!  ***************************************************************      
!                                                                       
!                                                                       
! INPUT:    H         -  Magnitudine assoluta                           
!           G         -  Parametro di slope                             
!           DS        -  Distanza dal Sole (au)                         
!           DT        -  Distanza dall'osservatore (au)                     
!           BETA      -  Angolo di fase solare (rad): angolo tra il     
!                        Sole e l'osservatore (Terra), visto dall'astero
!                                                                       
DOUBLE PRECISION FUNCTION appmag(h,g,ds,dt,beta) 
  IMPLICIT NONE
  DOUBLE PRECISION h,g,ds,dt,beta 
  DOUBLE PRECISION a1,a2,b1,b2 
  DOUBLE PRECISION tb2,phi1,phi2 
  SAVE a1,a2,b1,b2
! Costanti per il calcolo della magnitudine                             
  DATA a1,a2,b1,b2/3.33d0,1.87d0,0.63d0,1.22d0/ 
  tb2=TAN(beta/2.d0) 
  phi1=EXP(-a1*tb2**b1) 
  phi2=EXP(-a2*tb2**b2) 
  appmag=5.d0*LOG10(ds*dt)+h-2.5d0*LOG10((1.d0-g)*phi1+g*phi2) 
END FUNCTION appmag



