! Copyright (c) 2015-2021, the ELSI team.
! All rights reserved.
!
! This file is part of ELSI and is distributed under the BSD 3-clause license,
! which may be found in the LICENSE file in the ELSI root directory.

!>
!! Provide routines to set parameters of ELSI and the solvers.
!!
module ELSI_SET

   use ELSI_CONSTANT, only: MULTI_PROC
   use ELSI_DATATYPE, only: elsi_handle
   use ELSI_MALLOC, only: elsi_allocate
   use ELSI_MPI
   use ELSI_PRECISION, only: r8,i4
   use ELSI_UTIL, only: elsi_check_init

   implicit none

   private

   public :: elsi_set_output
   public :: elsi_set_output_unit
   public :: elsi_set_output_log
   public :: elsi_set_save_ovlp
   public :: elsi_set_unit_ovlp
   public :: elsi_set_zero_def
   public :: elsi_set_sparsity_mask
   public :: elsi_set_illcond_check
   public :: elsi_set_illcond_tol
   public :: elsi_set_spin_degeneracy
   public :: elsi_set_energy_gap
   public :: elsi_set_spectrum_width
   public :: elsi_set_dimensionality
   public :: elsi_set_extrapolation
   public :: elsi_set_elpa_solver
   public :: elsi_set_elpa_n_single
   public :: elsi_set_elpa_gpu
   public :: elsi_set_elpa_autotune
   public :: elsi_set_occ_non_aufbau
   public :: elsi_set_omm_flavor
   public :: elsi_set_omm_n_elpa
   public :: elsi_set_omm_tol
   public :: elsi_set_pexsi_method
   public :: elsi_set_pexsi_n_mu
   public :: elsi_set_pexsi_n_pole
   public :: elsi_set_pexsi_np_per_pole
   public :: elsi_set_pexsi_np_symbo
   public :: elsi_set_pexsi_temp
   public :: elsi_set_pexsi_mu_min
   public :: elsi_set_pexsi_mu_max
   public :: elsi_set_pexsi_inertia_tol
   public :: elsi_set_eigenexa_method
   public :: elsi_set_sips_n_elpa
   public :: elsi_set_sips_n_slice
   public :: elsi_set_sips_inertia_tol
   public :: elsi_set_sips_ev_min
   public :: elsi_set_sips_ev_max
   public :: elsi_set_ntpoly_method
   public :: elsi_set_ntpoly_isr
   public :: elsi_set_ntpoly_tol
   public :: elsi_set_ntpoly_filter
   public :: elsi_set_ntpoly_max_iter
   public :: elsi_set_magma_solver
   public :: elsi_set_mu_broaden_scheme
   public :: elsi_set_mu_broaden_width
   public :: elsi_set_mu_tol
   public :: elsi_set_mu_mp_order
   public :: elsi_set_n_frozen
   public :: elsi_set_frozen_idx
   public :: elsi_set_frozen_method

   ! Deprecated
   public :: elsi_set_write_unit
   public :: elsi_set_sing_check
   public :: elsi_set_sing_tol
   public :: elsi_set_elpa_gpu_kernels
   public :: elsi_set_pexsi_gap
   public :: elsi_set_pexsi_delta_e
   public :: elsi_set_sips_interval
   public :: elsi_set_mu_spin_degen
   public :: elsi_set_uuid

   interface elsi_set_write_unit
      module procedure elsi_set_output_unit
   end interface

   interface elsi_set_sing_check
      module procedure elsi_set_illcond_check
   end interface

   interface elsi_set_sing_tol
      module procedure elsi_set_illcond_tol
   end interface

   interface elsi_set_mu_spin_degen
      module procedure elsi_set_spin_degeneracy
   end interface

   interface elsi_set_pexsi_gap
      module procedure elsi_set_energy_gap
   end interface

   interface elsi_set_pexsi_delta_e
      module procedure elsi_set_spectrum_width
   end interface

contains

!>
!! Set the output level of ELSI.
!!
subroutine elsi_set_output(eh,output)

   implicit none

   type(elsi_handle), intent(inout) :: eh !< Handle
   integer(kind=i4), intent(in) :: output !< Output level

   character(len=*), parameter :: caller = "elsi_set_output"

   call elsi_check_init(eh%bh,eh%handle_init,caller)

   if(output < 0) then
      eh%bh%print_info = 0
   else if(output > 3) then
      eh%bh%print_info = 3
   else
      eh%bh%print_info = output
   end if

end subroutine

!>
!! Set the unit to be used by ELSI output.
!!
subroutine elsi_set_output_unit(eh,print_unit)

   implicit none

   type(elsi_handle), intent(inout) :: eh !< Handle
   integer(kind=i4), intent(in) :: print_unit !< Unit

   character(len=*), parameter :: caller = "elsi_set_output_unit"

   call elsi_check_init(eh%bh,eh%handle_init,caller)

   eh%bh%print_unit = print_unit

end subroutine

!>
!! Set whether a JSON log file should be output.
!!
subroutine elsi_set_output_log(eh,output_log)

   implicit none

   type(elsi_handle), intent(inout) :: eh !< Handle
   integer(kind=i4), intent(in) :: output_log !< Output log

   character(len=*), parameter :: caller = "elsi_set_output_log"

   call elsi_check_init(eh%bh,eh%handle_init,caller)

   eh%bh%print_json = output_log

end subroutine

!>
!! Set whether to save the overlap matrix for density matrix extrapolation.
!!
subroutine elsi_set_save_ovlp(eh,save_ovlp)

   implicit none

   type(elsi_handle), intent(inout) :: eh !< Handle
   integer(kind=i4), intent(in) :: save_ovlp !< Save overlap?

   character(len=*), parameter :: caller = "elsi_set_save_ovlp"

   call elsi_check_init(eh%bh,eh%handle_init,caller)

   if(save_ovlp == 0) then
      eh%ph%save_ovlp = .false.
   else
      eh%ph%save_ovlp = .true.
   end if

end subroutine

!>
!! Set the overlap matrix to be identity.
!!
subroutine elsi_set_unit_ovlp(eh,unit_ovlp)

   implicit none

   type(elsi_handle), intent(inout) :: eh !< Handle
   integer(kind=i4), intent(in) :: unit_ovlp !< Overlap is identity?

   character(len=*), parameter :: caller = "elsi_set_unit_ovlp"

   call elsi_check_init(eh%bh,eh%handle_init,caller)

   if(unit_ovlp == 0) then
      eh%ph%unit_ovlp = .false.
   else
      eh%ph%unit_ovlp = .true.
   end if

end subroutine

!>
!! Set the threshold to define "zero".
!!
subroutine elsi_set_zero_def(eh,zero_def)

   implicit none

   type(elsi_handle), intent(inout) :: eh !< Handle
   real(kind=r8), intent(in) :: zero_def !< Zero threshold

   character(len=*), parameter :: caller = "elsi_set_zero_def"

   call elsi_check_init(eh%bh,eh%handle_init,caller)

   eh%bh%def0 = zero_def

end subroutine

!>
!! Set the sparsity pattern used in dense to sparse matrix conversion.
!!
subroutine elsi_set_sparsity_mask(eh,sparsity_mask)

   implicit none

   type(elsi_handle), intent(inout) :: eh !< Handle
   integer(kind=i4), intent(in) :: sparsity_mask !< Sparsity pattern mask

   character(len=200) :: msg

   character(len=*), parameter :: caller = "elsi_set_sparsity_mask"

   call elsi_check_init(eh%bh,eh%handle_init,caller)

   if(sparsity_mask < 0 .or. sparsity_mask > 2) then
      write(msg,"(A)") "Input value should be 0, 1, or 2"
      call elsi_stop(eh%bh,msg,caller)
   end if

   eh%ph%sparsity_mask = sparsity_mask

end subroutine

!>
!! Set whether ill-conditioning should be checked.
!!
subroutine elsi_set_illcond_check(eh,illcond_check)

   implicit none

   type(elsi_handle), intent(inout) :: eh !< Handle
   integer(kind=i4), intent(in) :: illcond_check !< Check singularity?

   character(len=*), parameter :: caller = "elsi_set_illcond_check"

   call elsi_check_init(eh%bh,eh%handle_init,caller)

   if(illcond_check == 0) then
      eh%ph%ill_check = .false.
      eh%ph%ill_ovlp = .false.
      eh%ph%n_good = eh%ph%n_basis
   else
      eh%ph%ill_check = .true.
   end if

end subroutine

!>
!! Set the tolerance of ill-conditioning.
!!
subroutine elsi_set_illcond_tol(eh,illcond_tol)

   implicit none

   type(elsi_handle), intent(inout) :: eh !< Handle
   real(kind=r8), intent(in) :: illcond_tol !< Singularity tolerance

   character(len=*), parameter :: caller = "elsi_set_illcond_tol"

   call elsi_check_init(eh%bh,eh%handle_init,caller)

   eh%ph%ill_tol = illcond_tol

end subroutine

!>
!! Set the spin degeneracy that controls the maximum number of electrons on a
!! state.
!!
subroutine elsi_set_spin_degeneracy(eh,spin_degeneracy)

   implicit none

   type(elsi_handle), intent(inout) :: eh !< Handle
   real(kind=r8), intent(in) :: spin_degeneracy !< Spin degeneracy

   character(len=*), parameter :: caller = "elsi_set_spin_degeneracy"

   call elsi_check_init(eh%bh,eh%handle_init,caller)

   eh%ph%spin_degen = spin_degeneracy
   eh%ph%spin_is_set = .true.

end subroutine

!>
!! Set an estimate of the band (or HOMO/LUMO) gap.
!!
subroutine elsi_set_energy_gap(eh,gap)

   implicit none

   type(elsi_handle), intent(inout) :: eh !< Handle
   real(kind=r8), intent(in) :: gap !< Energy gap

   character(len=200) :: msg

   character(len=*), parameter :: caller = "elsi_set_energy_gap"

   call elsi_check_init(eh%bh,eh%handle_init,caller)

   if(gap < 0.0_r8) then
      write(msg,"(A)") "Input value cannot be negative"
      call elsi_stop(eh%bh,msg,caller)
   end if

   eh%ph%energy_gap = gap
   eh%ph%pexsi_options%gap = gap

end subroutine

!>
!! Set an estimate of the eigenspectrum width.
!!
subroutine elsi_set_spectrum_width(eh,width)

   implicit none

   type(elsi_handle), intent(inout) :: eh !< Handle
   real(kind=r8), intent(in) :: width !< Spectrum width

   character(len=200) :: msg

   character(len=*), parameter :: caller = "elsi_set_spectrum_width"

   call elsi_check_init(eh%bh,eh%handle_init,caller)

   if(width < 0.0_r8) then
      write(msg,"(A)") "Input value cannot be negative"
      call elsi_stop(eh%bh,msg,caller)
   end if

   eh%ph%spectrum_width = width
   eh%ph%pexsi_options%deltaE = width

end subroutine

!>
!! Set the dimensionality of the simulating system.
!!
subroutine elsi_set_dimensionality(eh,dimensionality)

   implicit none

   type(elsi_handle), intent(inout) :: eh !< Handle
   integer(kind=i4), intent(in) :: dimensionality !< Dimensionality

   character(len=200) :: msg

   character(len=*), parameter :: caller = "elsi_set_dimensionality"

   call elsi_check_init(eh%bh,eh%handle_init,caller)

   if(dimensionality < 1 .or. dimensionality > 3) then
      write(msg,"(A)") "Input value should be 1, 2, or 3"
      call elsi_stop(eh%bh,msg,caller)
   end if

   eh%ph%dimensionality = dimensionality

end subroutine

!>
!! Set the method for density matrix extrapolation.
!!
subroutine elsi_set_extrapolation(eh,extrapolation)

   implicit none

   type(elsi_handle), intent(inout) :: eh !< Handle
   integer(kind=i4), intent(in) :: extrapolation !< extrapolation

   character(len=200) :: msg

   character(len=*), parameter :: caller = "elsi_set_extrapolation"

   call elsi_check_init(eh%bh,eh%handle_init,caller)

   if(extrapolation < 0 .or. extrapolation > 1) then
      write(msg,"(A)") "Input value should be 0 or 1"
      call elsi_stop(eh%bh,msg,caller)
   end if

   eh%ph%extrapolation = extrapolation

end subroutine

!>
!! Set the ELPA solver.
!!
subroutine elsi_set_elpa_solver(eh,solver)

   implicit none

   type(elsi_handle), intent(inout) :: eh !< Handle
   integer(kind=i4), intent(in) :: solver !< ELPA solver

   character(len=200) :: msg

   character(len=*), parameter :: caller = "elsi_set_elpa_solver"

   call elsi_check_init(eh%bh,eh%handle_init,caller)

   if(solver < 1 .or. solver > 2) then
      write(msg,"(A)") "Input value should be 1 or 2"
      call elsi_stop(eh%bh,msg,caller)
   end if

   eh%ph%elpa_solver = solver

end subroutine

!>
!! Set the number of single precision steps with ELPA.
!!
subroutine elsi_set_elpa_n_single(eh,n_single)

   implicit none

   type(elsi_handle), intent(inout) :: eh !< Handle
   integer(kind=i4), intent(in) :: n_single !< Single-precision steps

   character(len=*), parameter :: caller = "elsi_set_elpa_n_single"

   call elsi_check_init(eh%bh,eh%handle_init,caller)

   if(n_single < 0) then
      eh%ph%elpa_n_single = 0
   else
      eh%ph%elpa_n_single = n_single
   end if

end subroutine

!>
!! Set whether GPU acceleration eigenvectors should be enabled in ELPA. No
!! effect if no GPU acceleration available.
!!
subroutine elsi_set_elpa_gpu(eh,gpu)

   implicit none

   type(elsi_handle), intent(inout) :: eh !< Handle
   integer(kind=i4), intent(in) :: gpu !< Use GPU acceleration?

   character(len=*), parameter :: caller = "elsi_set_elpa_gpu"

   call elsi_check_init(eh%bh,eh%handle_init,caller)

   if(gpu == 0) then
      eh%ph%elpa_gpu = 0
   else
      eh%ph%elpa_gpu = 1
   end if

end subroutine

!>
!! Set whether auto-tuning should be enabled in ELPA. No effect if no
!! auto-tuning available.
!!
subroutine elsi_set_elpa_autotune(eh,autotune)

   implicit none

   type(elsi_handle), intent(inout) :: eh !< Handle
   integer(kind=i4), intent(in) :: autotune !< Use auto-tuning?

   character(len=*), parameter :: caller = "elsi_set_elpa_autotune"

   call elsi_check_init(eh%bh,eh%handle_init,caller)

   if(autotune == 0) then
      eh%ph%elpa_autotune = 0
   else
      eh%ph%elpa_autotune = 1
   end if

end subroutine

!>
!! Set the flavor of libOMM.
!!
subroutine elsi_set_omm_flavor(eh,flavor)

   implicit none

   type(elsi_handle), intent(inout) :: eh !< Handle
   integer(kind=i4), intent(in) :: flavor !< libOMM flavor

   character(len=200) :: msg

   character(len=*), parameter :: caller = "elsi_set_omm_flavor"

   call elsi_check_init(eh%bh,eh%handle_init,caller)

   if(flavor /= 0 .and. flavor /= 2) then
      write(msg,"(A)") "Input value should be 0 or 2"
      call elsi_stop(eh%bh,msg,caller)
   end if

   eh%ph%omm_flavor = flavor

end subroutine

!>
!! Set the number of ELPA steps when using libOMM.
!!
subroutine elsi_set_omm_n_elpa(eh,n_elpa)

   implicit none

   type(elsi_handle), intent(inout) :: eh !< Handle
   integer(kind=i4), intent(in) :: n_elpa !< ELPA steps

   character(len=*), parameter :: caller = "elsi_set_omm_n_elpa"

   call elsi_check_init(eh%bh,eh%handle_init,caller)

   if(n_elpa < 1) then
      eh%ph%omm_n_elpa = 1
   else
      eh%ph%omm_n_elpa = n_elpa
   end if

end subroutine

!>
!! Set the convergence tolerance of OMM minimization.
!!
subroutine elsi_set_omm_tol(eh,tol)

   implicit none

   type(elsi_handle), intent(inout) :: eh !< Handle
   real(kind=r8), intent(in) :: tol !< Tolerance of OMM minimization

   character(len=200) :: msg

   character(len=*), parameter :: caller = "elsi_set_omm_tol"

   call elsi_check_init(eh%bh,eh%handle_init,caller)

   if(tol <= 0.0_r8) then
      write(msg,"(A)") "Input value should be positive"
      call elsi_stop(eh%bh,msg,caller)
   end if

   eh%ph%omm_tol = tol

end subroutine

!>
!! Set the pole expansion method in PEXSI.
!!
subroutine elsi_set_pexsi_method(eh,method)

   implicit none

   type(elsi_handle), intent(inout) :: eh !< Handle
   integer(kind=i4), intent(in) :: method !< Pole expansion method

   character(len=200) :: msg

   character(len=*), parameter :: caller = "elsi_set_pexsi_method"

   call elsi_check_init(eh%bh,eh%handle_init,caller)

   if(method < 1 .or. method > 3) then
      write(msg,"(A)") "Input value should be 1, 2, or 3"
      call elsi_stop(eh%bh,msg,caller)
   end if

   eh%ph%pexsi_options%method = method

end subroutine

!>
!! Set the number of chemical potential points in PEXSI.
!!
subroutine elsi_set_pexsi_n_mu(eh,n_mu)

   implicit none

   type(elsi_handle), intent(inout) :: eh !< Handle
   integer(kind=i4), intent(in) :: n_mu !< Number of mu points

   character(len=200) :: msg

   character(len=*), parameter :: caller = "elsi_set_pexsi_n_mu"

   call elsi_check_init(eh%bh,eh%handle_init,caller)

   if(n_mu < 1) then
      write(msg,"(A)") "Input value cannot be smaller than 1"
      call elsi_stop(eh%bh,msg,caller)
   end if

   eh%ph%pexsi_options%nPoints = n_mu

end subroutine

!>
!! Set the number of poles in the pole expansion.
!!
subroutine elsi_set_pexsi_n_pole(eh,n_pole)

   implicit none

   type(elsi_handle), intent(inout) :: eh !< Handle
   integer(kind=i4), intent(in) :: n_pole !< Number of poles

   character(len=200) :: msg

   character(len=*), parameter :: caller = "elsi_set_pexsi_n_pole"

   call elsi_check_init(eh%bh,eh%handle_init,caller)

   if(n_pole < 1) then
      write(msg,"(A)") "Input value cannot be smaller than 1"
      call elsi_stop(eh%bh,msg,caller)
   end if

   eh%ph%pexsi_options%numPole = n_pole

end subroutine

!>
!! Set the number of MPI tasks assigned for one pole.
!!
subroutine elsi_set_pexsi_np_per_pole(eh,np_per_pole)

   implicit none

   type(elsi_handle), intent(inout) :: eh !< Handle
   integer(kind=i4), intent(in) :: np_per_pole !< Number of tasks per pole

   character(len=200) :: msg

   character(len=*), parameter :: caller = "elsi_set_pexsi_np_per_pole"

   call elsi_check_init(eh%bh,eh%handle_init,caller)

   if(np_per_pole < 1) then
      write(msg,"(A)") "Input value cannot be smaller than 1"
      call elsi_stop(eh%bh,msg,caller)
   end if

   eh%ph%pexsi_np_per_pole = np_per_pole

end subroutine

!>
!! Set the number of MPI tasks for the symbolic factorization.
!!
subroutine elsi_set_pexsi_np_symbo(eh,np_symbo)

   implicit none

   type(elsi_handle), intent(inout) :: eh !< Handle
   integer(kind=i4), intent(in) :: np_symbo !< Number of tasks for symbolic

   character(len=200) :: msg

   character(len=*), parameter :: caller = "elsi_set_pexsi_np_symbo"

   call elsi_check_init(eh%bh,eh%handle_init,caller)

   if(np_symbo < 1) then
      write(msg,"(A)") "Input value cannot be smaller than 1"
      call elsi_stop(eh%bh,msg,caller)
   end if

   eh%ph%pexsi_options%npSymbFact = np_symbo

end subroutine

!>
!! Set the electronic temperature in PEXSI.
!!
subroutine elsi_set_pexsi_temp(eh,temp)

   implicit none

   type(elsi_handle), intent(inout) :: eh !< Handle
   real(kind=r8), intent(in) :: temp !< Temperature

   character(len=200) :: msg

   character(len=*), parameter :: caller = "elsi_set_pexsi_temp"

   call elsi_check_init(eh%bh,eh%handle_init,caller)

   if(temp <= 0.0_r8) then
      write(msg,"(A)") "Input value should be positive"
      call elsi_stop(eh%bh,msg,caller)
   end if

   eh%ph%pexsi_options%temperature = temp

end subroutine

!>
!! Set the lower bound of the chemical potential in PEXSI.
!!
subroutine elsi_set_pexsi_mu_min(eh,mu_min)

   implicit none

   type(elsi_handle), intent(inout) :: eh !< Handle
   real(kind=r8), intent(in) :: mu_min !< Lower bound of mu

   character(len=*), parameter :: caller = "elsi_set_pexsi_mu_min"

   call elsi_check_init(eh%bh,eh%handle_init,caller)

   eh%ph%pexsi_options%muMin0 = mu_min

end subroutine

!>
!! Set the upper bound of the chemical potential in PEXSI.
!!
subroutine elsi_set_pexsi_mu_max(eh,mu_max)

   implicit none

   type(elsi_handle), intent(inout) :: eh !< Handle
   real(kind=r8), intent(in) :: mu_max !< Upper bound of mu

   character(len=*), parameter :: caller = "elsi_set_pexsi_mu_max"

   call elsi_check_init(eh%bh,eh%handle_init,caller)

   eh%ph%pexsi_options%muMax0 = mu_max

end subroutine

!>
!! Set the tolerance of the estimation of the chemical potential in the inertia
!! counting procedure.
!!
subroutine elsi_set_pexsi_inertia_tol(eh,inertia_tol)

   implicit none

   type(elsi_handle), intent(inout) :: eh !< Handle
   real(kind=r8), intent(in) :: inertia_tol !< Inertia counting tolerance

   character(len=200) :: msg

   character(len=*), parameter :: caller = "elsi_set_pexsi_inertia_tol"

   call elsi_check_init(eh%bh,eh%handle_init,caller)

   if(inertia_tol <= 0.0_r8) then
      write(msg,"(A)") "Input value should be positive"
      call elsi_stop(eh%bh,msg,caller)
   end if

   eh%ph%pexsi_options%muInertiaTolerance = inertia_tol

end subroutine

!>
!! Set the EigenExa algorithm.
!!
subroutine elsi_set_eigenexa_method(eh,method)

   implicit none

   type(elsi_handle), intent(inout) :: eh !< Handle
   integer(kind=i4), intent(in) :: method !< Method

   character(len=200) :: msg

   character(len=*), parameter :: caller = "elsi_set_eigenexa_method"

   call elsi_check_init(eh%bh,eh%handle_init,caller)

   if(method < 1 .or. method > 2) then
      write(msg,"(A)") "Input value should be 1 or 2"
      call elsi_stop(eh%bh,msg,caller)
   end if

   eh%ph%exa_method = method

end subroutine

!>
!! Set the number of ELPA steps when using SLEPc-SIPs.
!!
subroutine elsi_set_sips_n_elpa(eh,n_elpa)

   implicit none

   type(elsi_handle), intent(inout) :: eh !< Handle
   integer(kind=i4), intent(in) :: n_elpa !< ELPA steps

   character(len=*), parameter :: caller = "elsi_set_sips_n_elpa"

   call elsi_check_init(eh%bh,eh%handle_init,caller)

   if(n_elpa < 1) then
      eh%ph%sips_n_elpa = 1
   else
      eh%ph%sips_n_elpa = n_elpa
   end if

end subroutine

!>
!! Set the number of slices in SLEPc-SIPs.
!!
subroutine elsi_set_sips_n_slice(eh,n_slice)

   implicit none

   type(elsi_handle), intent(inout) :: eh !< Handle
   integer(kind=i4), intent(in) :: n_slice !< Number of slices

   character(len=200) :: msg

   character(len=*), parameter :: caller = "elsi_set_sips_n_slice"

   call elsi_check_init(eh%bh,eh%handle_init,caller)

   if(n_slice < 1) then
      write(msg,"(A)") "Input value cannot be smaller than 1"
      call elsi_stop(eh%bh,msg,caller)
   end if

   eh%ph%sips_n_slices = n_slice

end subroutine

!>
!! Set the tolerance to stop inertia counting in SLEPc-SIPs.
!!
subroutine elsi_set_sips_inertia_tol(eh,inertia_tol)

   implicit none

   type(elsi_handle), intent(inout) :: eh !< Handle
   real(kind=r8), intent(in) :: inertia_tol !< Stopping criterion

   character(len=200) :: msg

   character(len=*), parameter :: caller = "elsi_set_sips_inertia_tol"

   call elsi_check_init(eh%bh,eh%handle_init,caller)

   if(inertia_tol <= 0.0_r8) then
      write(msg,"(A)") "Input value should be positive"
      call elsi_stop(eh%bh,msg,caller)
   end if

   eh%ph%sips_inertia_tol = inertia_tol

end subroutine

!>
!! Set the lower bound of the interval to be solved by SLEPc-SIPs.
!!
subroutine elsi_set_sips_ev_min(eh,ev_min)

   implicit none

   type(elsi_handle), intent(inout) :: eh !< Handle
   real(kind=r8), intent(in) :: ev_min !< Lower bound

   character(len=*), parameter :: caller = "elsi_set_sips_ev_min"

   call elsi_check_init(eh%bh,eh%handle_init,caller)

   eh%ph%sips_interval(1) = ev_min

end subroutine

!>
!! Set the upper bound of the interval to be solved by SLEPc-SIPs.
!!
subroutine elsi_set_sips_ev_max(eh,ev_max)

   implicit none

   type(elsi_handle), intent(inout) :: eh !< Handle
   real(kind=r8), intent(in) :: ev_max !< Upper bound

   character(len=*), parameter :: caller = "elsi_set_sips_ev_max"

   call elsi_check_init(eh%bh,eh%handle_init,caller)

   eh%ph%sips_interval(2) = ev_max

end subroutine

!>
!! Set the global interval to be solved by SLEPc-SIPs.
!!
subroutine elsi_set_sips_interval(eh,lower,upper)

   implicit none

   type(elsi_handle), intent(inout) :: eh !< Handle
   real(kind=r8), intent(in) :: lower !< Lower bound
   real(kind=r8), intent(in) :: upper !< Upper bound

   character(len=200) :: msg

   character(len=*), parameter :: caller = "elsi_set_sips_interval"

   call elsi_check_init(eh%bh,eh%handle_init,caller)

   if(lower >= upper) then
      write(msg,"(A)") "Lower bound should be smaller than upper bound"
      call elsi_stop(eh%bh,msg,caller)
   end if

   eh%ph%sips_interval(1) = lower
   eh%ph%sips_interval(2) = upper

end subroutine

!>
!! Set the density matrix purification method in NTPoly.
!!
subroutine elsi_set_ntpoly_method(eh,method)

   implicit none

   type(elsi_handle), intent(inout) :: eh !< Handle
   integer(kind=i4), intent(in) :: method !< Purification method

   character(len=200) :: msg

   character(len=*), parameter :: caller = "elsi_set_ntpoly_method"

   call elsi_check_init(eh%bh,eh%handle_init,caller)

   if(method < 0 .or. method > 3) then
      write(msg,"(A)") "Input value should be 0, 1, 2, or 3"
      call elsi_stop(eh%bh,msg,caller)
   end if

   eh%ph%nt_method = method

end subroutine

!>
!! Set the inverse square root method in NTPoly.
!!
subroutine elsi_set_ntpoly_isr(eh,isr)

   implicit none

   type(elsi_handle), intent(inout) :: eh !< Handle
   integer(kind=i4), intent(in) :: isr !< Inverse square root method

   character(len=200) :: msg

   character(len=*), parameter :: caller = "elsi_set_ntpoly_isr"

   call elsi_check_init(eh%bh,eh%handle_init,caller)

   if(isr /= 2 .and. isr /= 3 .and. isr /= 5) then
      write(msg,"(A)") "Input value should be 2, 3, or 5"
      call elsi_stop(eh%bh,msg,caller)
   end if

   eh%ph%nt_isr = isr

end subroutine

!>
!! Set the tolerance of the density matrix purification.
!!
subroutine elsi_set_ntpoly_tol(eh,tol)

   implicit none

   type(elsi_handle), intent(inout) :: eh !< Handle
   real(kind=r8), intent(in) :: tol !< Tolerance

   character(len=200) :: msg

   character(len=*), parameter :: caller = "elsi_set_ntpoly_tol"

   call elsi_check_init(eh%bh,eh%handle_init,caller)

   if(tol <= 0.0_r8) then
      write(msg,"(A)") "Input value should be positive"
      call elsi_stop(eh%bh,msg,caller)
   end if

   eh%ph%nt_tol = tol

end subroutine

!>
!! Set the threshold to filter intermediate results in NTPoly.
!!
subroutine elsi_set_ntpoly_filter(eh,filter)

   implicit none

   type(elsi_handle), intent(inout) :: eh !< Handle
   real(kind=r8), intent(in) :: filter !< Filter

   character(len=*), parameter :: caller = "elsi_set_ntpoly_filter"

   call elsi_check_init(eh%bh,eh%handle_init,caller)

   eh%ph%nt_filter = filter

end subroutine

!>
!! Set the maximum number of density matrix purification steps.
!!
subroutine elsi_set_ntpoly_max_iter(eh,max_iter)

   implicit none

   type(elsi_handle), intent(inout) :: eh !< Handle
   integer(kind=i4), intent(in) :: max_iter !< Number of iterations

   character(len=200) :: msg

   character(len=*), parameter :: caller = "elsi_set_ntpoly_max_iter"

   call elsi_check_init(eh%bh,eh%handle_init,caller)

   if(max_iter < 1) then
      write(msg,"(A)") "Input value cannot be smaller than 1"
      call elsi_stop(eh%bh,msg,caller)
   end if

   eh%ph%nt_max_iter = max_iter

end subroutine

!>
!! Set the solver of MAGMA.
!!
subroutine elsi_set_magma_solver(eh,solver)

   implicit none

   type(elsi_handle), intent(inout) :: eh !< Handle
   integer(kind=i4), intent(in) :: solver !< MAGMA solver

   character(len=200) :: msg

   character(len=*), parameter :: caller = "elsi_set_magma_solver"

   call elsi_check_init(eh%bh,eh%handle_init,caller)

   if(solver < 1 .or. solver > 2) then
      write(msg,"(A)") "Input value should be 1 or 2"
      call elsi_stop(eh%bh,msg,caller)
   end if

   eh%ph%magma_solver = solver

end subroutine

!>
!! Set the broadening scheme to determine the chemical potential and the
!! occupation numbers.
!!
subroutine elsi_set_mu_broaden_scheme(eh,broaden_scheme)

   implicit none

   type(elsi_handle), intent(inout) :: eh !< Handle
   integer(kind=i4), intent(in) :: broaden_scheme !< Broadening method

   character(len=200) :: msg

   character(len=*), parameter :: caller = "elsi_set_mu_broaden_method"

   call elsi_check_init(eh%bh,eh%handle_init,caller)

   if(broaden_scheme < 0 .or. broaden_scheme > 4) then
      write(msg,"(A)") "Input value should be 0, 1, 2, 3, or 4"
      call elsi_stop(eh%bh,msg,caller)
   end if

   eh%ph%mu_scheme = broaden_scheme

end subroutine

!>
!! Set the broadening width to determine the chemical potential and the
!! occupation numbers.
!!
subroutine elsi_set_mu_broaden_width(eh,broaden_width)

   implicit none

   type(elsi_handle), intent(inout) :: eh !< Handle
   real(kind=r8), intent(in) :: broaden_width !< Broadening width

   character(len=200) :: msg

   character(len=*), parameter :: caller = "elsi_set_mu_broaden_width"

   call elsi_check_init(eh%bh,eh%handle_init,caller)

   if(broaden_width <= 0.0_r8) then
      write(msg,"(A)") "Input value should be positive"
      call elsi_stop(eh%bh,msg,caller)
   end if

   eh%ph%mu_width = broaden_width

end subroutine

!>
!! Set the desired accuracy of the determination of the chemical potential and
!! the occupation numbers.
!!
subroutine elsi_set_mu_tol(eh,tol)

   implicit none

   type(elsi_handle), intent(inout) :: eh !< Handle
   real(kind=r8), intent(in) :: tol !< Accuracy of mu

   character(len=200) :: msg

   character(len=*), parameter :: caller = "elsi_set_mu_tol"

   call elsi_check_init(eh%bh,eh%handle_init,caller)

   if(tol < 0.0_r8) then
      write(msg,"(A)") "Input value cannot be negative"
      call elsi_stop(eh%bh,msg,caller)
   end if

   eh%ph%mu_tol = tol

end subroutine
!>
!! Set the non Aufbau occupation
!!
subroutine elsi_set_occ_non_aufbau(eh,occ_non_aufbau)

   implicit none

   type(elsi_handle), intent(inout) :: eh !< Handle
   logical, intent(in) :: occ_non_aufbau !< whether to use non Aufbau occupation

   character(len=200) :: msg

   character(len=*), parameter :: caller = "elsi_set_occ_non_aufbau"

   call elsi_check_init(eh%bh,eh%handle_init,caller)

   eh%ph%occ_non_aufbau = occ_non_aufbau

end subroutine

!>
!! Set the order of the Methfessel-Paxton broadening scheme.
!!
subroutine elsi_set_mu_mp_order(eh,mp_order)

   implicit none

   type(elsi_handle), intent(inout) :: eh !< Handle
   integer(kind=i4), intent(in) :: mp_order !< Order

   character(len=200) :: msg

   character(len=*), parameter :: caller = "elsi_set_mu_mp_order"

   call elsi_check_init(eh%bh,eh%handle_init,caller)

   if(mp_order < 0) then
      write(msg,"(A)") "Input value cannot be negative"
      call elsi_stop(eh%bh,msg,caller)
   end if

   eh%ph%mu_mp_order = mp_order

end subroutine

!>
!! Set the number of frozen core orbitals.
!!
subroutine elsi_set_n_frozen(eh,n_frozen)

   implicit none

   type(elsi_handle), intent(inout) :: eh !< Handle
   integer(kind=i4), intent(in) :: n_frozen !< Number of frozen orbitals

   integer(kind=i4) :: ierr
   character(len=200) :: msg

   integer(kind=i4), external :: numroc

   character(len=*), parameter :: caller = "elsi_set_n_frozen"

   call elsi_check_init(eh%bh,eh%handle_init,caller)

   if(n_frozen < 0) then
      write(msg,"(A)") "Input value cannot be negative"
      call elsi_stop(eh%bh,msg,caller)
   end if

   eh%ph%n_basis_c = n_frozen
   eh%ph%n_basis_v = eh%ph%n_basis-n_frozen

   if(eh%ph%parallel_mode == MULTI_PROC) then
      eh%ph%n_lrow_v = numroc(eh%ph%n_basis_v,eh%bh%blk,eh%bh%my_prow,0,&
         eh%bh%n_prow)
      eh%ph%n_lcol_v = numroc(eh%ph%n_basis_v,eh%bh%blk,eh%bh%my_pcol,0,&
         eh%bh%n_pcol)

      call descinit(eh%ph%desc_v,eh%ph%n_basis_v,eh%ph%n_basis_v,eh%bh%blk,&
           eh%bh%blk,0,0,eh%bh%blacs_ctxt,max(1,eh%ph%n_lrow_v),ierr)
   end if

   call elsi_allocate(eh%bh,eh%perm_fc,eh%ph%n_basis,"perm_fc",caller)

end subroutine

!>
!! Set the indices of frozen core orbitals.
!!
subroutine elsi_set_frozen_idx(eh,frozen_idx)

   implicit none

   type(elsi_handle), intent(inout) :: eh !< Handle
   integer(kind=i4), intent(in) :: frozen_idx(eh%ph%n_basis_c) !< Indices

   integer :: i
   integer :: j
   integer :: tmp
   character(len=200) :: msg

   character(len=*), parameter :: caller = "elsi_set_frozen_idx"

   call elsi_check_init(eh%bh,eh%handle_init,caller)

   if(eh%ph%n_basis_c > 0) then
      if(any(frozen_idx <= 0)) then
         write(msg,"(A)") "Input value should be positive"
         call elsi_stop(eh%bh,msg,caller)
      end if

      if(any(frozen_idx > eh%ph%n_basis)) then
         write(msg,"(A)") "Input value cannot be larger than number of basis"
         call elsi_stop(eh%bh,msg,caller)
      end if

      if(any(frozen_idx > eh%ph%n_basis_c)) then
         eh%ph%fc_perm = .true.

         do i = 1,eh%ph%n_basis
            eh%perm_fc(i) = i
         end do

         do i = 1,eh%ph%n_basis_c
            j = frozen_idx(i)
            tmp = eh%perm_fc(i)
            eh%perm_fc(i) = eh%perm_fc(j)
            eh%perm_fc(j) = tmp
         end do
      end if
   end if

end subroutine

!>
!! Set the frozen core approximation method.
!!
subroutine elsi_set_frozen_method(eh,frozen_method)

   implicit none

   type(elsi_handle), intent(inout) :: eh !< Handle
   integer(kind=i4), intent(in) :: frozen_method !< Frozen core method

   character(len=200) :: msg

   character(len=*), parameter :: caller = "elsi_set_frozen_method"

   call elsi_check_init(eh%bh,eh%handle_init,caller)

   if(frozen_method < 0 .or. frozen_method > 2) then
      write(msg,"(A)") "Input value should be 0, 1, or 2"
      call elsi_stop(eh%bh,msg,caller)
   end if

   eh%ph%fc_method = frozen_method

end subroutine

!>
!! Set a UUID.
!! (Deprecated)
!!
subroutine elsi_set_uuid(eh,uuid)

   implicit none

   type(elsi_handle), intent(inout) :: eh !< Handle
   character(len=*), intent(in) :: uuid !< UUID

   character(len=36) :: tmp

   character(len=*), parameter :: caller = "elsi_set_uuid"

   call elsi_check_init(eh%bh,eh%handle_init,caller)

   tmp = uuid

end subroutine

!>
!! Set whether GPU kernels for back-transforming eigenvectors should be enabled
!! in ELPA.
!! (Deprecated)
!!
subroutine elsi_set_elpa_gpu_kernels(eh,gpu_kernels)

   implicit none

   type(elsi_handle), intent(inout) :: eh !< Handle
   integer(kind=i4), intent(in) :: gpu_kernels !< Use GPU kernels?

   integer(kind=i4) :: tmp

   character(len=*), parameter :: caller = "elsi_set_elpa_gpu_kernels"

   call elsi_check_init(eh%bh,eh%handle_init,caller)

   tmp = gpu_kernels

end subroutine

end module ELSI_SET


text:
      module m
      implicit none
      type vector
      integer::x
      integer::y
      end type vector
      interface operator (.plus.)
      module procedure addvec
      end interface operator (.plus.)
      interface assignment (=)
      module procedure assignvec
      endinterface assignment (=)
      contains
      function addvec(a,b)
      type(vector),intent(in)::a,b
      type(vector)::addvec
      addvec%x=a%x+b%x
      addvec%y=a%y+b%y
      end function addvec
      subroutine assignvec (b,a)
      type(vector),intent(in)::a
      real,intent(out)::b
      b = sqrt((a%x**2.0) + (a%y**2.0))
      endsubroutine assignvec

      end module m

      program t
      use m
      implicit none
      type(vector)::f,g
      integer::h,i
      real::magnitude
      f%x=7
      f%y=8
      g%x=9
      g%y=10
      h=7
      i=8
      print *,f.plus.g
      print *,h+i
      magnitude =f
      print *,magnitude
      end program t

ext: f

max_stars_repo_name: maddenp/ppp

max_issues_repo_name: maddenp/ppp

max_forks_repo_name: maddenp/ppp


text:
    program test
    USE, INTRINSIC :: ISO_C_BINDING
    USE OPENST
    use omp_lib
    REAL (C_DOUBLE), TARGET, ALLOCATABLE :: U(:,:,:), V(:,:,:)
    INTEGER (C_SIZE_T) :: NI , NJ, NK
    REAL (C_DOUBLE) :: HI, HJ, HK, SRCI, SRCJ, SRCK, RCVI, RCVJ, RCVK
    TYPE (C_PTR) :: U_PTR, V_PTR
    REAL (C_DOUBLE) :: di, dj, dk, dist
    INTEGER (C_INT) :: errcode
    INTEGER (C_SIZE_T):: i,j,k
    INTEGER (C_SIZE_T) :: MAX_SEG
    double precision :: start, finish

    REAL (C_DOUBLE) :: TSTEP
    REAL (C_DOUBLE), POINTER, DIMENSION(:,:) :: RAY
    TYPE (C_PTR), TARGET :: RAY_PTR, RAY_PTR_PTR
    INTEGER (C_SIZE_T), TARGET :: RAY_NI, RAY_NJ
    TYPE (C_PTR) :: RAY_NI_PTR, RAY_NJ_PTR

    NI = 50
    NJ = 50
    NK = 50

    RCVI = 1.0
    RCVJ = 0.5
    RCVK = 0.5

    SRCI = 0.5
    SRCJ = 0.5
    SRCK = 0.5

    HI = 1.0 / (NI - 1)
    HJ = 1.0 / (NJ - 1)
    HK = 1.0 / (NK - 1)

    TSTEP = MIN(HI,HJ,HK) / 1.0 * 0.999
    MAX_SEG = (NI * NJ * NK);

    print '("H: ",E15.6,E15.6,E15.6)', HI,HJ,HK
    print '("TSTEP: ",E15.6)', TSTEP
    print '("MAX_SEG: ",I22)', MAX_SEG

    ALLOCATE ( U (NI, NJ, NK) )
    ALLOCATE ( V (NI, NJ, NK) )

    U_PTR = C_LOC(U)
    V_PTR = C_LOC(V)

    do k = 1,NK
        do j = 1,NJ
            do i = 1,NI
                V(i,j,k) = 1.0
                di = SRCI - DBLE(i - 1) * HI
                dj = SRCJ - DBLE(j - 1) * HJ
                dk = SRCK - DBLE(k - 1) * HK
                dist = SQRT(di * di + dj * dj + dk * dk)
                U(i,j,k) = dist
            end do
        end do
    end do

    RAY_PTR_PTR = C_LOC(RAY_PTR)
    RAY_NI_PTR = C_LOC(RAY_NI)
    RAY_NJ_PTR = C_LOC(RAY_NJ)

    start = omp_get_wtime()
    errcode = OpenST_BRT3D_Trace(U_PTR,V_PTR,NI,NJ,NK,HI,HJ,HK,TSTEP,RCVI,RCVJ,RCVK, &
        SRCI,SRCJ,SRCK,MAX_SEG,RAY_PTR_PTR,RAY_NI_PTR,RAY_NJ_PTR)
    finish = omp_get_wtime()
    print '("BRT3D: ",E15.6," seconds.")',finish-start

    if (errcode.NE.0) then
        print '("OPENST_ERROR: code ",i2)',errcode
        STOP 1
    end if

    CALL C_F_POINTER(RAY_PTR, RAY, [RAY_NI, RAY_NJ])
    RAY = RESHAPE( RAY, [RAY_NI,RAY_NJ], order=[2,1])

    do i = 1,RAY_NI
        print '(I3, E15.6, E15.6, E15.6)', i, RAY(i,1), RAY(i,2), RAY(i,3)
    end do

    CALL OpenST_INTEROP_FreePointer(RAY_PTR)

    end program

max_stars_repo_name: aanikitin/seis-eikonal

max_issues_repo_name: aanikitin/seis-eikonal

max_forks_repo_name: aanikitin/seis-eikonal


text:
!	Find estimated reliabilities of eaps of latent vectors.
!	Input:

!	covtheta is the array of conditional covariance matrices.
!	meantheta is the array of conditional means.

!	Output:

!	obsweight provides observation weights.
!	coveap is the covariance matrix of the eaps.
!	covartheta is the covariance matrix of the latent vector.
!	meancovtheta is the average conditional covariance of theta given the observations.
!	meaneap is the average eap.
!	releap is the reliability for each element of the latent vector.

subroutine reltheta(covtheta,meantheta,obsweight,coveap,covartheta,meancovtheta,meaneap,releap)
implicit none
interface
!	Multinomial mean.
	function mmeanth(density,quadpoint)
		implicit none
		real(kind=8),intent(in)::density(:),quadpoint(:,:)
		real(kind=8)::mmeanth(size(quadpoint,1))
	end function mmeanth
!	Multinomial covariance matrix.
	function mcovth(density,meanth,quadpoint)
		implicit none
		real(kind=8),intent(in)::density(:),meanth(:),quadpoint(:,:)
		real(kind=8)::mcovth(size(meanth),size(meanth))
	end function mcovth

end interface
real(kind=8),intent(in)::covtheta(:,:,:),meantheta(:,:),obsweight(:)
real(kind=8),intent(out)::coveap(:,:),covartheta(:,:),meancovtheta(:,:),meaneap(:),releap(:)
!	row is row number.
integer::row
!	Average conditional variance, mean eap, variance of eaps, total weight
real(kind=8)::condvar,mean,var,total
total=sum(obsweight)
releap=0.0_8
meaneap=mmeanth(obsweight,meantheta)/total
coveap=mcovth(obsweight,meaneap,meantheta)/total
do row=1,size(meaneap)
	meancovtheta(:,row)=mmeanth(obsweight,covtheta(:,row,:))/total
	covartheta(:,row)=meancovtheta(:,row)+coveap(:,row)
	if(covartheta(row,row)>0.0_8)releap(row)=coveap(row,row)/covartheta(row,row)
end do

return
end subroutine reltheta

max_stars_repo_name: EducationalTestingService/MIRT

max_issues_repo_name: EducationalTestingService/MIRT

max_forks_repo_name: EducationalTestingService/MIRT


text:
subroutine intpar(neqn,nout,nsteps,t0,tf,abserr,relerr)
!
!  Subroutine intpar sets the parameters to control the
!  integration of the linear PDE problem
!
!  Double precision coding is used
implicit double precision(a-h,o-z)
!
!  Number of ODEs
neqn=21
!
!  Number of output points
nout=6
!
!  Maximum number of steps in the interval t0 to tf
nsteps=250
!
!  Initial, final values of the independent variable
t0=0.0d0
tf=0.2d0
!
!  Error tolerances
abserr=1.0d-05
relerr=1.0d-05
return
!
!  End of intpar
end

max_stars_repo_name: cunyizju/Runge-Kutta-ODE-FORTRAN

max_issues_repo_name: cunyizju/Runge-Kutta-ODE-FORTRAN

max_forks_repo_name: cunyizju/Runge-Kutta-ODE-FORTRAN


text:
C Copyright(C) 2009-2017 National Technology & Engineering Solutions
C of Sandia, LLC (NTESS).  Under the terms of Contract DE-NA0003525 with
C NTESS, the U.S. Government retains certain rights in this software.
C
C Redistribution and use in source and binary forms, with or without
C modification, are permitted provided that the following conditions are
C met:
C
C     * Redistributions of source code must retain the above copyright
C       notice, this list of conditions and the following disclaimer.
C
C     * Redistributions in binary form must reproduce the above
C       copyright notice, this list of conditions and the following
C       disclaimer in the documentation and/or other materials provided
C       with the distribution.
C
C     * Neither the name of NTESS nor the names of its
C       contributors may be used to endorse or promote products derived
C       from this software without specific prior written permission.
C
C THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
C "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
C LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
C A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
C OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
C SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
C LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
C DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
C THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
C (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
C OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

C=======================================================================
      SUBROUTINE SQZSTR (STRING, LSTR)
C=======================================================================

C   --*** SQZSTR *** (STRLIB) Remove extra blanks from string
C   --   Written by Amy Gilkey - revised 06/02/87
C   --
C   --SQZSTR deletes leading and extra blanks within a string.
C   --To prevent problems, an empty string is returned with a length of 1.
C   --
C   --Parameters:
C   --   STRING - IN/OUT - the string to be compressed, returned, may be
C   --      up to 1024 characters long
C   --   LSTR - OUT - the new string length

C   --Routines Called:
C   --   LENSTR - (STRLIB) Find string length

      CHARACTER*(*) STRING
      INTEGER LSTR

      PARAMETER (MXSTR = 1024)
      CHARACTER*(MXSTR) TMPSTR

      LSTR = 1
      IF (STRING .EQ. ' ') RETURN

      LSTR = LENSTR (STRING)
      if (lstr .gt. MXSTR) then
        call prterr ('PROGRAM', 'String is too long in SQZSTR')
        return
      end if

      IBLK = INDEX (STRING, '  ')
  100 CONTINUE
      IF ((IBLK .GT. 0) .AND. (IBLK .LE. LSTR)) THEN

         INOBLK = IBLK + 2
  110    CONTINUE
         IF (STRING(INOBLK:INOBLK) .EQ. ' ') THEN
            INOBLK = INOBLK + 1
            GOTO 110
         END IF

         TMPSTR = STRING(INOBLK:LSTR)
         STRING(IBLK+1:LSTR) = TMPSTR
         LSTR = LSTR - (INOBLK-IBLK-1)
         IBLK = INDEX (STRING, '  ')
         GOTO 100
      END IF

      IF (STRING(1:1) .EQ. ' ') THEN
         TMPSTR = STRING(2:LSTR)
         STRING(1:LSTR) = TMPSTR
         LSTR = LSTR - 1
      END IF

      RETURN
      END

max_stars_repo_name: mathstuf/seacas

max_stars_repo_licenses: BSD-3-Clause NetCDF Zlib MIT

max_issues_repo_name: mathstuf/seacas

max_issues_repo_licenses: BSD-3-Clause NetCDF Zlib MIT

max_forks_repo_name: mathstuf/seacas

max_forks_repo_licenses: BSD-3-Clause NetCDF Zlib MIT


MODULE xc_ggapbe96_2_module

  use gradient_module
  use grid_module, only: grid
  use xc_variables, only: xcpot, xcene
  use fd_module, only: fd, construct_nabla_fd
  use lattice_module, only: lattice, get_aa_lattice, get_reciprocal_lattice
  use parallel_module
  use basic_type_factory

  implicit none

  PRIVATE
  PUBLIC :: calc_GGAPBE96_2

  integer,parameter :: DP=kind(0.0d0)
!#ifdef _NO_QPRECISION_
  integer,parameter :: QP=kind(0.0d0)
!#else
!  integer,parameter :: QP=kind(0.0q0)
!#endif

  real(QP),parameter :: zero_density = 1.e-10_QP
  real(QP),allocatable :: nab(:)
  real(QP),allocatable :: vx(:,:),vc(:,:)
  real(QP) :: Ex,Ec
  real(QP) :: b(3,3)
  integer :: Md
  integer :: SYStype=0

  real(8) :: mu=0.21951d0
  real(8) :: Kp=0.804d0
  real(8) :: beta=0.066725d0

CONTAINS

  SUBROUTINE calc_GGAPBE96_2( rgrid, rho, ene, pot, mu_in, Kp_in )

    implicit none
    type( grid ),intent(IN) :: rgrid
    type( GSArray ),intent(IN) :: rho
    type( xcene ) :: ene
    type( xcpot ) :: pot
    real(8),optional,intent(IN) :: mu_in, Kp_in
    type(gradient16) :: grad16
    type(fd) :: nabla
    type(lattice) :: aa, bb
    integer :: m1,m2,n1,n2,i,i1,i2,i3
    real(QP) :: Pi
    real(DP) :: sb(2),rb(2)

! ---

    if ( present(mu_in) ) mu=mu_in
    if ( present(Kp_in) ) Kp=Kp_in
    beta = mu*3.0d0/acos(-1.0d0)**2

    call construct_gradient16( rgrid, rho, grad16 )

    call construct_nabla_fd( nabla )

    Md = nabla%md

    if ( .not.allocated(nab) ) allocate( nab(-Md:Md) )
    nab(:) = nabla%coef(:)

    call get_aa_lattice( aa )
    call get_reciprocal_lattice( aa, bb )

    Pi=acos(-1.0_QP)
    b(1:3,1)=aa%Length(1)*bb%LatticeVector(1:3,1)/( 2*Pi*rgrid%spacing(1) )
    b(1:3,2)=aa%Length(2)*bb%LatticeVector(1:3,2)/( 2*Pi*rgrid%spacing(2) )
    b(1:3,3)=aa%Length(3)*bb%LatticeVector(1:3,3)/( 2*Pi*rgrid%spacing(3) )

! ---

    m1 = pot%xc%g_range%head
    m2 = pot%xc%g_range%tail
    n1 = rho%s_range%head
    n2 = rho%s_range%tail

    allocate( vx(m1:m2,n1:n2) ) ; vx=0.0_QP

    call calc_PBE_x( rgrid, rho, grad16 )

    allocate( vc(m1:m2,n1:n2) ) ; vc=0.0_QP

    call calc_PBE_c( rgrid, rho, grad16 )

! ---

    sb(1)=Ex*rgrid%VolumeElement
    sb(2)=Ec*rgrid%VolumeElement
    call MPI_ALLREDUCE( sb, rb, 2, MPI_REAL8, MPI_SUM, comm_grid, i )

    ene%Ex  = rb(1)
    ene%Ec  = rb(2)
    ene%Exc = rb(1)+rb(2)

    n1 = pot%xc%s_range%head
    n2 = pot%xc%s_range%tail

    pot%xc%val(:,:)  = vx(:,n1:n2) + vc(:,n1:n2)
    if ( allocated(pot%x%val) ) pot%x%val(:,:) = vx(:,n1:n2)
    if ( allocated(pot%c%val) ) pot%c%val(:,:) = vc(:,n1:n2)

! ---

    deallocate( vc )
    deallocate( vx )

    call destruct_gradient16( grad16 )

  END SUBROUTINE calc_GGAPBE96_2

  SUBROUTINE calc_PBE_x( rgrid, rho, grad16 )

    implicit none
    type( grid ) :: rgrid
    type( GSArray ) :: rho
    type( gradient16 ) :: grad16
!   real(8),parameter :: mu=0.21951d0, Kp=0.804d0
    integer :: i,ispin,m,ierr
    real(QP) :: trho,cm
    real(QP),allocatable :: rrrr(:,:),rtmp(:)
    real(QP) :: kf, vx_lda, ex_lda, Fx, Pi, g2, factor
    real(QP) :: onethr,const1,const2
    integer :: i1,i2,i3,j,j1,j2,j3,k1,k2,k3
    integer :: mm,m1,m2
    real(8),allocatable :: f(:), Gf(:)

    Pi = acos(-1.0_QP)

    factor = 1.0_QP
    if ( rho%s_range%size_global == 2 ) factor = 2.0_QP

    onethr = 1.0_QP/3.0_QP
    const1 = 3.0_QP*Pi*Pi
    const2 = 3.0_QP/(4.0_QP*Pi)

    m1 = rho%g_range%head
    m2 = rho%g_range%tail
    allocate( rtmp(m1:m2) ) ; rtmp=0.0_QP

    Ex = 0.0_QP

    do ispin=rho%s_range%head,rho%s_range%tail

       do i=rho%g_range%head,rho%g_range%tail

          trho = factor*rho%val(i,ispin)

          if ( trho <= zero_density ) cycle

          kf = (const1*trho)**onethr

          ex_lda = -const2*kf
          vx_lda = -kf/Pi

          g2 = grad16%gg(i)

          Fx = 1.0_QP + Kp - 4.0_QP*Kp*Kp*(trho*kf)**2 &
                            /( 4.0_QP*Kp*(trho*kf)**2 + mu*g2 )

          Ex = Ex + trho*ex_lda*Fx

          vx(i,ispin) = vx(i,ispin) &
               + Fx*vx_lda + ( 24.0_QP*Pi*Kp*Kp*mu*trho**3*g2 ) &
                            /( 4.0_QP*Kp*(trho*kf)**2 + mu*g2 )**2

          rtmp(i) = -18.0_QP*Pi*Kp*Kp*mu*trho**4 &
                    /( 4.0_QP*Kp*(trho*kf)**2 + mu*g2 )**2

       end do ! i

       allocate( rrrr(m1:m2,3)  ) ; rrrr=0.0_QP
       rrrr(m1:m2,1) = rtmp(m1:m2)*grad16%gx(m1:m2)
       rrrr(m1:m2,2) = rtmp(m1:m2)*grad16%gy(m1:m2)
       rrrr(m1:m2,3) = rtmp(m1:m2)*grad16%gz(m1:m2)

       allocate( f(m1:m2)  ) ; f=0.0d0
       allocate( Gf(m1:m2) ) ; Gf=0.0d0

       do j=1,3
          f(:) = rrrr(m1:m2,1)*b(1,j) + rrrr(m1:m2,2)*b(2,j) + rrrr(m1:m2,3)*b(3,j)
          call calc_abc_gradient( j, rgrid, f, Gf )
          vx(:,ispin) = vx(:,ispin) - Gf(:)
       end do

       deallocate( Gf )
       deallocate( f  )
       deallocate( rrrr )

    end do ! ispin

    Ex = Ex/factor

    deallocate( rtmp )

  END SUBROUTINE calc_PBE_x

  SUBROUTINE calc_PBE_c( rgrid, rho, grad16 )

    implicit none
    type( grid ) :: rgrid
    type( GSArray ) :: rho
    type( gradient16 ) :: grad16
    real(8),parameter :: A00  =0.031091d0,A01  =0.015545d0,A02  =0.016887d0
    real(8),parameter :: alp10=0.21370d0 ,alp11=0.20548d0 ,alp12=0.11125d0
    real(8),parameter :: bt10 =7.5957d0  ,bt11 =1.41189d1 ,bt12 =1.0357d1
    real(8),parameter :: bt20 =3.5876d0  ,bt21 =6.1977d0  ,bt22 =3.6231d0
    real(8),parameter :: bt30 =1.6382d0  ,bt31 =3.3662d0  ,bt32 =0.88026d0
    real(8),parameter :: bt40 =0.49294d0 ,bt41 =0.62517d0 ,bt42 =0.49671d0
!   real(8),parameter :: C1=2.14611945579107d0,C2=0.031091d0
    real(QP) :: C1,C2
    real(QP) :: const1, const2, factor
    integer :: i,j,ispin,m,i1,i2,i3,j1,j2,j3,k1,k2,k3,ierr
    integer :: m1,m2
    real(QP) :: trho, rhoa, rhob
    real(QP) :: kf, rs, ec_U, ec_P, ec_lda, phi, g2
    real(QP) :: dac_dn, dfz_dz, deU_dn, deP_dn, H, A, T, alpc, fz
    real(QP) :: dec_dz, dphi_dz, dH_dA, dH_dT, tmp, dH_dphi, Ai
    real(QP) :: dA_dn, dec_dn, Hs, zeta, dT_dphi
    real(QP) :: dz_dn(2), cm, rssq
    real(QP) :: dac_drs, deP_drs, deU_drs, drs_dn
    real(QP),allocatable :: rrrr(:,:), rtmp(:)
    real(QP) :: Pi, one, two, fouthr, onethr, zero
    real(QP) :: sevthr, twothr, thrtwo, ThrFouPi
    real(8),allocatable :: f(:),Gf(:)

    const1 = 2.0_QP**(4.0_QP/3.0_QP)-2.0_QP
    const2 = 9.0_QP*(2.0_QP**(1.0_QP/3.0_QP)-1.0_QP)/4.0_QP

    factor = 1.0_QP
    if ( rho%s_range%size_global == 1 ) factor = 0.5_QP

    Pi       = acos(-1.0_QP)
    zero     = 0.0_QP
    one      = 1.0_QP
    two      = 2.0_QP
    fouthr   = 4.0_QP/3.0_QP
    onethr   = 1.0_QP/3.0_QP
    ThrFouPi = 3.0_QP/(4.0_QP*Pi)
    thrtwo   = 3.0_QP/2.0_QP
    twothr   = 2.0_QP/3.0_QP
    sevthr   = 7.0_QP/3.0_QP

    C2 = ( 1.0_QP-log(2.0_QP) )/Pi**2 ! "gamma" in PBE paper
    C1 = beta/C2

    Ec = 0.0_QP

    m1 = rho%g_range%head
    m2 = rho%g_range%tail
    allocate( rtmp(m1:m2) ) ; rtmp=0.0_QP

    do i=rho%g_range%head,rho%g_range%tail

       rhoa = rho%val(i,1)*factor
       rhob = rho%val(i,rho%s_range%size_global)*factor
       trho = rhoa + rhob

       if ( trho <= zero_density ) cycle

       zeta = ( rhoa - rhob )/trho
       if ( abs(zeta) > one ) zeta=sign(one,zeta)

       fz = ( (one+zeta)**fouthr + (one-zeta)**fouthr - two )*const1

       kf = ( 3.0_QP*Pi*Pi*trho )**onethr

       rs = ( ThrFouPi/trho )**onethr

       rssq = sqrt(rs)

       ec_U = -two*A00*( one + alp10*rs ) &
            *log( one + one/( two*A00 &
            *(bt10*rssq + bt20*rs + bt30*rs*rssq + bt40*rs*rs) ) )
       ec_P = -two*A01*( one + alp11*rs ) &
            *log( one + one/( two*A01 &
            *(bt11*rssq + bt21*rs + bt31*rs*rssq + bt41*rs*rs) ) )
       alpc = -two*A02*( one + alp12*rs ) &
            *log( one + one/( two*A02 &
            *(bt12*rssq + bt22*rs + bt32*rs*rssq + bt42*rs*rs) ) )

       ec_lda = ec_U - alpc*fz*const2*(one-zeta**4) &
            + ( ec_P - ec_U )*fz*zeta**4

       phi = 0.5_QP*( (one+zeta)**twothr + (one-zeta)**twothr )

       g2 = grad16%gg(i)

       T = Pi*g2/( 16.0_QP*phi**2*kf*trho**2 )

       A = C1/( exp( -ec_lda/(C2*phi**3) ) - one )

       Hs = C2 * phi**3 * log( one + C1*T*(one+A*T)/(one+A*T+(A*T)**2) )

       Ec = Ec + trho*( ec_lda + Hs )

       drs_dn = -4.0_QP*Pi*rs**4/9.0_QP

       tmp = bt10*rssq + bt20*rs + bt30*rs*rssq + bt40*rs*rs
       deU_drs = alp10*ec_U/( one + alp10*rs ) &
            +A00*( one + alp10*rs )/rssq &
            *( bt10 + two*bt20*rssq + 3.0_QP*bt30*rs + 4.0_QP*bt40*rs*rssq ) &
            /( two*A00*tmp*tmp+tmp )

       tmp = bt11*rssq + bt21*rs + bt31*rs*rssq + bt41*rs*rs
       deP_drs = alp11*ec_P/( one + alp11*rs ) &
            +A01*( one + alp11*rs )/rssq &
            *( bt11 + two*bt21*rssq + 3.0_QP*bt31*rs + 4.0_QP*bt41*rs*rssq ) &
            /( two*A01*tmp*tmp+tmp )

       tmp = bt12*rssq + bt22*rs + bt32*rs*rssq + bt42*rs*rs
       dac_drs = alp12*alpc/( one + alp12*rs ) &
            +A02*( one + alp12*rs )/rssq &
            *( bt12 + two*bt22*rssq + 3.0_QP*bt32*rs + 4.0_QP*bt42*rs*rssq ) &
            /( two*A02*tmp*tmp+tmp )

       deU_dn = deU_drs * drs_dn
       deP_dn = deP_drs * drs_dn
       dac_dn = dac_drs * drs_dn

       dfz_dz = fouthr*( (one+zeta)**onethr - (one-zeta)**onethr )*const1

       dec_dz = -alpc*dfz_dz*const2*(one-zeta**4) &
               + 4.0_QP*alpc*fz*const2*zeta**3 &
               +(ec_P-ec_U)*( dfz_dz*zeta**4 + fz*4.0_QP*zeta**3 )

       if ( abs(zeta) == one ) then
          dphi_dz = zero
       else
          dphi_dz = ( (one+zeta)**(-onethr)-(one-zeta)**(-onethr) )/3.0_QP
       end if

       tmp = one + A*T + (A*T)**2

       dH_dA = -phi**3*C1*C2*A*T**3*(two+A*T)/(tmp**2+C1*T*(one+A*T)*tmp)

       dH_dT =  phi**3*C1*C2*(one+two*A*T)/(tmp**2+C1*T*(one+A*T)*tmp)

       dH_dphi = 3.0_QP*Hs/phi

       dz_dn(1)   = 2.0_QP*rhob/trho**2
       dz_dn(rho%s_range%size_global) =-2.0_QP*rhoa/trho**2

       do ispin=rho%s_range%head,rho%s_range%tail

          dec_dn = deU_dn - dac_dn*fz*const2*(one-zeta**4) &
               +(deP_dn-deU_dn)*fz*zeta**4 + dec_dz*dz_dn(ispin)

          dA_dn = A*(C1+A)/(C1*C2*phi**3) &
               *( dec_dn - 3.0_QP*ec_lda/phi*dphi_dz*dz_dn(ispin) )

          vc(i,ispin) = vc(i,ispin) + ec_lda + Hs + trho*dec_dn &
               + trho*dH_dA*dA_dn &
               + dH_dT*(-sevthr*T - trho*two*T/phi*dphi_dz*dz_dn(ispin) ) &
               + trho*dH_dphi*dphi_dz*dz_dn(ispin)

       end do ! ispin

       rtmp(i) = dH_dT*Pi/(8.0_QP*kf*trho*phi**2)

    end do ! i

    allocate( rrrr(m1:m2,3)  ) ; rrrr=0.0_QP
    rrrr(m1:m2,1) = rtmp(m1:m2)*grad16%gx(m1:m2)
    rrrr(m1:m2,2) = rtmp(m1:m2)*grad16%gy(m1:m2)
    rrrr(m1:m2,3) = rtmp(m1:m2)*grad16%gz(m1:m2)

    allocate( f(m1:m2)  ) ; f=0.0d0
    allocate( Gf(m1:m2) ) ; Gf=0.0d0

    do j=1,3
       f(:) = rrrr(m1:m2,1)*b(1,j) + rrrr(m1:m2,2)*b(2,j) + rrrr(m1:m2,3)*b(3,j)
       call calc_abc_gradient( j, rgrid, f, Gf )
       do ispin=rho%s_range%head,rho%s_range%tail
          vc(:,ispin) = vc(:,ispin) - Gf(:)
       end do
    end do

    deallocate( Gf )
    deallocate( f  )
    deallocate( rrrr )
    deallocate( rtmp )

    return
  END SUBROUTINE calc_PBE_c

END MODULE xc_ggapbe96_2_module


text: The Stonegate Tennis Academy provides an opportunity for people of all ages to experience the sport of tennis in personal and group settings. They offer private one on one lessons, clinics for individuals and small groups, camps throughout the year and tout one of Californias top ranked http://www.usta.com/ USTA junior teams. Their goal: to help you strive for excellence in every aspect of the sport, whether it be as an individual or incorporating you into one of their teams.

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
ccc   By Trifon Trifonov trifon@hku.hk
ccc   You can modify if as you want, but please
ccc   do not distribute without permision.
ccc   This is not a final version!!!
ccc   The final version will be available in the Python RVMod lib
ccc   Trifonov et al. (in prep).

      implicit none
      real*8 PI, twopi
      parameter (PI=3.14159265358979d0)
      integer npl,ndset,idset,ndata,ma,mfit,i,j,NDSMAX,NPLMAX,MMAX
      integer writeflag_best_par,hkl,gr_flag
      integer writeflag_RV,writeflag_fit, amoebastarts
      parameter (NDSMAX=20, NPLMAX=20, MMAX=200)
      integer idsmax(NDSMAX),ia(MMAX),nt, ts(20000),ii, iter
      real*8 x(20000),y(20000),sig(20000),y_in(20000)
      real*8 a(MMAX),covar(MMAX,MMAX),alpha(MMAX,MMAX)
      real*8 rms,mstar, mass(NPLMAX),ap(NPLMAX)
      real*8 swift_mass(NPLMAX),s_mass(NPLMAX),j_mass(NPLMAX)
      real*8 chisq,alamda,ochisq,dchisq, epsil, deltat
      real*8 sigscale,x0,xmax, incl(NPLMAX),cap0m(NPLMAX)
      real*8 t0,t1,t2,dt,offset,t_max,loglik,dy,sig2i
      real*8 st_mass,sini,m1,a1,m2,a2,epoch,ftol
      real*8 ymod(20000),dyda(MMAX), p(MMAX+1,MMAX),yamoeba(MMAX+1)
      real*8 loglikk, ologlikk, dloglikk,best_w,best_we
      external rvkep, compute_abs_loglik
      character*80 infile
      character*80 version_input, version

      real*4 t_stop,when_to_kill, model_max,model_min

      common /DSBLK/ npl,ndset,idsmax,idset,gr_flag

      version = "0.07"

      CALL getarg(1, version_input)
      if(version_input.eq.'-version') then
          write(*,*) version
          goto 222
      endif

      twopi=2.d0*PI
      ftol=0.000001d0

c     first two just for consistency with dynamical input, not really used
      read (*,*) epsil,deltat, amoebastarts,
     &          when_to_kill, nt, model_max, model_min ,gr_flag

c      write(*,*) 'Stellar mass'
      read (*,*) st_mass, writeflag_best_par, writeflag_RV,
     & writeflag_fit

      call io_read_data (ndata,x,ts,y,sig,epoch,
     &               x0,t_max,a,ia,ma,incl,cap0m,hkl)

      mfit = 0
      do j = 1,ma
          if (ia(j).ne.0) mfit = mfit + 1
      enddo

c      call prepare_for_amoeba(p,MMAX+1,MMAX,yamoeba,a,ia,ma,mfit,
c     & compute_abs_loglik,ndata,x,y,ymod,dyda,ts,sig)

      i = 0
 500  continue

         if (i.eq.amoebastarts) then
             i = 0
             goto 502
         endif

         i = i + 1
         ologlikk = loglikk

         call prepare_for_amoeba(p,MMAX+1,MMAX,yamoeba,a,ia,ma,mfit,
     & compute_abs_loglik,ndata,x,y,ymod,dyda,ts,sig, i,hkl)
         call amoeba(p,yamoeba,MMAX+1,MMAX,mfit,ftol,compute_abs_loglik,
     & iter,ndata,x,y,ymod,dyda,ma,ts,sig,a,ia,loglikk,hkl)

         CALL SECOND(t_stop)
         if (t_stop.ge.when_to_kill) then
            write(*,*) 'Max. time=',when_to_kill, 'sec ',
     &                 'exceeded t_stop =', t_stop, 'sec '
            goto 502
         endif

         loglikk = yamoeba(1)
         dloglikk = ologlikk - loglikk

         j=0
         do ii=1,ma
           if (ia(ii).ne.0) then
              j=j+1
              a(ii)=p(1,j)
           endif
         enddo

      if (dabs(dloglikk).ge.0.000001d0) goto 500

502   idset = 1
      chisq=0.d0
      loglik=0.d0

      if (hkl.eq.0) then

          do i = 1,npl
             j = 6*(i-1)

             if (a(j+2).lt.0.d0) then  ! if P<0, set P>0
                a(j+2) = abs(a(j+2))
             endif

             if (a(j+1).lt.0.d0) then  ! if K<0, set K>0 and w = w+PI
                a(j+4) = a(j+4) + PI
                a(j+1) = abs(a(j+1))
                if (a(j+4).gt.2.d0*PI) a(j+4) = a(j+4)-2.d0*PI
             endif
             if (a(j+3).lt.0.d0) then  ! if e<0, set e>0 and w=w+PI, M0=M0-PI
                a(j+3) = abs(a(j+3))
                a(j+4) = a(j+4) +  PI
                if (a(j+4).gt.2.d0*PI) a(j+4) = a(j+4)-2.d0*PI
                a(j+5) = a(j+5) - PI
                if (a(j+5).lt.0.d0) a(j+5) = a(j+5)+2.d0*PI
             endif
             if (a(j+4).lt.0.d0) a(j+4) = dmod(a(j+4)+2.d0*PI, 2.d0*PI)
             if (a(j+5).lt.0.d0) a(j+5) = dmod(a(j+5)+2.d0*PI, 2.d0*PI)
c             if (a(j+6).lt.0.d0) a(j+6) = dmod(a(j+6)+2.d0*PI, 2.d0*PI)

             if (a(j+4).gt.2.d0*PI) a(j+4) = dmod(a(j+4), 2.d0*PI )
             if (a(j+5).gt.2.d0*PI) a(j+5) = dmod(a(j+5), 2.d0*PI )
c             if (a(j+6).gt.2.d0*PI) a(j+6) = dmod(a(j+6), 2.d0*PI )

          enddo

      else

          do i = 1,npl
             j = 6*(i-1)
             if (a(j+1).lt.0.d0) then  ! if K<0, set K>0 and w = w+PI
                a(j+4) = -1.d0*a(j+4)       !     which is h = -h, k = -k
                a(j+3) = -1.d0*a(j+3)
                a(j+1) = abs(a(j+1))
             endif

             if (a(j+5).lt.0.d0) a(j+5) = dmod(a(j+5)+2.d0*PI, 2.d0*PI)
             if (a(j+5).gt.2.d0*PI) a(j+5) = dmod(a(j+5), 2.d0*PI )
             if (a(j+6).lt.0.d0) a(j+6) = dmod(a(j+6)+2.d0*PI, 2.d0*PI)
             if (a(j+6).gt.2.d0*PI) a(j+6) = dmod(a(j+6), 2.d0*PI )

c             write(*,*) a(j+4),a(j+4),ecc(i) ,omega(i) ,capmm(i)
          enddo

      endif

      do i = 1,ndata

          idset = ts(i)
          call RVKEP (x(i),a,ymod(i),dyda,ma,idset,hkl)

          y_in(i) = y(i) - a(6*npl+idset) - a(6*npl+2*ndset+1)*x(i) -
     &      a(6*npl +2*ndset + 2)*x(i)**2
	      ymod(i) = ymod(i) - a(6*npl+idset)
     &    - a(6*npl +2*ndset + 1)*x(i) -
     &      a(6*npl +2*ndset + 2)*x(i)**2

          dy = y_in(i) - ymod(i)

          if (writeflag_RV.gt.0) then
              write(*,*) x0 + x(i),
     &        ymod(i), y_in(i) + a(6*npl+2*ndset+1)*x(i) +
     &        a(6*npl +2*ndset + 2)*x(i)**2,
     &        dy, sig(i), idset

          endif

          sig2i = 1.d0/(sig(i)**2 + a(6*npl+ndset+idset)**2)

 	      chisq  = chisq + dy*dy*sig2i

	      loglik =  loglik - 0.5*dy*dy*sig2i -
     &               0.5*dlog(twopi*(sig(i)**2
     &                + a(6*npl+ndset+idset)**2))
          rms = rms + dy**2
      enddo

      rms = dsqrt(rms/dble(ndata))

      write(*,*) 'loglik, reduced chi^2, chi^2, rms:'
      write(*,*) loglik, chisq/dble(ndata-mfit),chisq,rms

51    format(f10.3,f10.3,f10.3,f10.3,f10.3,f10.3,f10.3,f10.3)
52    format(a,f14.3)
53    format(a,i4,a,i4,a,f7.3,a,f7.3,a,f12.3)

      call MA_J (a,ma,npl,st_mass,sini,mass,ap,hkl,gr_flag)

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

      if (writeflag_best_par.gt.0) then

          do j = 1,npl+1
	      j_mass(j) = mass(j)/1.2667d17
c             s_mass(j) = mass(j)/1.32712497d20
c              swift_mass(j) = (mass(j)/1.32712497d20)*((4.d0*PI*PI)
c     &           /(365.25*365.25))
          enddo

          write (*,*) 'Best-fit K [m/s], P [days], e, w [deg],
     & M0 [deg], i[deg], cap0m[deg], w dot [deg/yr], and their errors'

          do j = 1,npl
              i = 6*(j-1)

              if (hkl.eq.0) then
                  best_w = a(i+4)*180.d0/PI
                  best_we = dsqrt(covar(i+4,i+4))*180.d0/PI
              else
                  best_w = a(i+4)
                  best_we = dsqrt(covar(i+4,i+4))
              endif

              write (*,*) a(i+1),a(i+2),a(i+3),best_w,
     &        a(i+5)*180.d0/PI,incl(j),cap0m(j),a(i+6)*180.d0/PI
              write (*,*) dsqrt(covar(i+1,i+1)),dsqrt(covar(i+2,i+2)),
     &                 dsqrt(covar(i+3,i+3)),
     &                 best_we,
     &                 dsqrt(covar(i+5,i+5))*180.d0/PI,0.d0, 0.d0,
     &                 dsqrt(covar(i+6,i+6))*180.d0/PI

          enddo

          write (*,*) 'Best-fit V0 [m/s] and their error bars:'

          do j = 1,ndset
              i = 6*npl + j
              write (*,*) a(i)
              write (*,*) dsqrt(covar(i,i))
          enddo

          write (*,*) 'Jitters for each data set:'
          do j = 1,ndset
              write (*,*) a(6*npl+ndset+j)
              write (*,*) '0'
          enddo

          write (*,*) 'linear trend [m/s per day]:'
          write (*,*) a(6*npl + 2*ndset + 1)
          write (*,*) dsqrt(covar(6*npl + ndset + 1,6*npl + ndset + 1))

          write (*,*) 'quad. trend [m/s per day]:'
          write (*,*) a(6*npl + 2*ndset + 2)
          write (*,*) dsqrt(covar(6*npl + ndset + 2,6*npl + ndset + 2))

          write (*,*) ' ndata =',ndata
          write (*,*) ' mfit =',mfit
          write (*,*) ' RMS =',rms
          write (*,*) ' Chi^2 =',chisq/dble(ndata-mfit)
          write (*,*) ' epoch = ', x0

          write (*,*) 'Jupiter mass'
          write (*,*) (j_mass(i+1),i=1,npl)

          write(*,*) 'semi-major axes in Jacobi'
          write(*,*)  (ap(i)/1.49597892d11,i=1,npl)

      endif

      if(writeflag_fit.gt.0) then
          dt = (x(ndata)+model_max+model_min)/dble(nt - 1)
          do i = 1,nt
	      x(i) = ((i-1)*dt)-model_min

              do j = 1,ndset
                  a(6*npl + j) = 0.0
              enddo
              call RVKEP (x(i),a,ymod(i),dyda,ma,1,hkl)
              write(*,*) x0 + x(i), ymod(i)
          enddo
      endif

c      stop
222   end

      subroutine compute_abs_loglik(ndata,x,y,a2,ymod,dyda,ma,mfit,ts,
     & sig,loglik, num,a,ia,hkl)
      implicit none

      integer MMAX,NDSMAX,npl,ndset,idset,num, mfit,gr_flag
      parameter (MMAX=200, NDSMAX=20)
      real*8 loglik, PI, TWOPI
      parameter (PI=3.14159265358979d0)
      parameter (TWOPI=2.0*PI)
      integer ndata, i, j, ma, ts(20000), ia(MMAX), idsmax(NDSMAX),hkl
      real*8 dy, sig(20000), dyda(MMAX), x(20000), y(20000)
      real*8 ymod(20000),a(MMAX),a2(mfit),a3(MMAX),sig2i,y_in(20000)
     & , y2(20000)

      common /DSBLK/ npl,ndset,idsmax,idset,gr_flag

      loglik=0.d0
      j=0
      do i=1,ma
         if (ia(i).ne.0) then
             j=j+1
             a3(i)=a2(j)
         else
             a3(i)=a(i)
         endif
      enddo

        do i = 1,ndata
              idset = ts(i)
              call RVKEP (x(i),a3,y2(i),dyda,ma,idset,hkl)
              y_in(i) = y(i) - a3(6*npl+idset)-
     &                 a3(6*npl+2*ndset+1)*x(i)
     &               - a3(6*npl+2*ndset+2)*x(i)**2

	          y2(i) = y2(i) - a3(6*npl+idset) -
     &         a3(6*npl+2*ndset+1)*x(i)
     &       - a3(6*npl+2*ndset+2)*x(i)**2

          dy = y_in(i) - y2(i)

	      sig2i = 1.d0/(sig(i)**2 + a3(6*npl+ndset+idset)**2)

	      loglik =  loglik + 0.5*dy*dy*sig2i +
     &               dlog(dsqrt(TWOPI*(sig(i)**2 +
     &               a3(6*npl+ndset+idset)**2)))
     &               - dlog(dsqrt(TWOPI))
c        write(*,*) loglik
        enddo

      return
      end

      subroutine io_read_data(ndata,t,ts,ys,sigs,epoch,t0,t_max,
     &          ar,iar,ma,incl,cap0m,hkl)

      implicit none
      integer ndset,idset,ndata,NDSMAX,NPLMAX,MMAX,npl,ma

      real*8 t(20000),y(20000),sig(20000),ys(20000),sigs(20000),PI
      parameter (NDSMAX=20,NPLMAX=20,MMAX=200)
      parameter(PI=3.14159265358979d0)
      real*8 ar(MMAX),incl(NPLMAX),cap0m(NPLMAX)
      integer iar(MMAX),u_off(NDSMAX),u_jit(NDSMAX),hkl
      integer idsmax(NDSMAX),ts(20000), u_incl, u_cap0m
      real*8 jitt(NDSMAX),sigscale,t0,t_max, epoch
      real*8 off(NDSMAX),loglik
      integer i,k,j,gr_flag
      character*80 infile

      common /DSBLK/ npl,ndset,idsmax,idset,gr_flag

c      write (*,*)  ' Number of Data Sets: '
      read (*,*) ndset
      if (ndset.gt.NDSMAX) stop ' KEPFIT: ndset > NDSMAX.'

      ndata = 1
      do i = 1,ndset
          read (*,50) infile
 50       format (a)
          read (*,*) off(i)
          read (*,*) u_off(i)
          read (*,*) jitt(i)
          read (*,*) u_jit(i)
          open (unit=10,file=infile)
 100      continue
             read (10,*,err=200,end=200) t0,y(ndata),sig(ndata)

c              sig(ndata) = dsqrt(sig(ndata)**2 + jitt(i)**2)
              if (ndata.eq.1) then              ! make sequences for datapoints
                 t(ndata) = t0
                 ts(ndata) = i
                 ys(ndata) = y(ndata)
                 sigs(ndata) = sig(ndata)
              else
                 j = 1
 1002            if (t0.lt.t(j)) then
                    do k = ndata-1,j,-1
                       t(k+1) = t(k)
                       ts(k+1) = ts(k)
                       ys(k+1) = ys(k)
                       sigs(k+1) = sigs(k)
                    enddo
                    t(j) = t0
                    ts(j) = i
                    ys(j) = y(ndata)
                    sigs(j) = sig(ndata)
                    goto 1001
                 else
                    j = j + 1
                    if (j.eq.ndata) then
                       t(j) = t0
                       ts(j) = i
                       ys(j) = y(ndata)
                       sigs(j) = sig(ndata)
                       goto 1001
                    endif
                    goto 1002
                 endif
              endif

 1001         ndata = ndata + 1
          goto 100
 200      idsmax(i) = ndata - 1
          close (unit=10)

          enddo
          read (*,*) npl
      if (npl.gt.NPLMAX) stop ' KEPFIT: npl > NPLMAX.'

      do i = 1,ndset
          ar(6*npl+i)=off(i)
          iar(6*npl+i)=u_off(i)
          ar(6*npl+ndset+i)=jitt(i)
          iar(6*npl+ndset+i)=u_jit(i)
      enddo

      ma = 6*npl + 2*ndset + 2
      do j = 1,npl
          i = 6*(j-1)
          read (*,*) ar(i+1),ar(i+2),ar(i+3),ar(i+4),ar(i+5),incl(j),
     &    cap0m(j),ar(i+6)
          read (*,*) iar(i+1),iar(i+2),iar(i+3),iar(i+4),iar(i+5),
     &    u_incl, u_cap0m,iar(i+6)

c         inclinations and cap0m are always ignored in the fit, just for consistency with dynamical input and output

      enddo

      read (*,*) ar(6*npl+ 2*ndset+1)
      read (*,*) iar(6*npl+2*ndset+1)

      read (*,*) ar(6*npl+ 2*ndset+2)
      read (*,*) iar(6*npl+2*ndset+2)

      ndata = ndata - 1

c      write(*,*) 'for epoch :'
      read (*,*) epoch

      t_max = t(ndata)

      if (epoch.eq.0) then
         t0 = t(1)
      else
         t0 = epoch
      endif

      read (*,*) hkl

      do j = 1,npl
          i = 6*(j-1)
          if (hkl.eq.0) then
              ar(i+4) = ar(i+4)*PI/180.d0
          endif
          ar(i+5) = ar(i+5)*PI/180.d0
          ar(i+6) = ar(i+6)*PI/180.d0
      enddo

      do i = 1,ndata
         t(i) = (t(i) - t0)             ! time unit is day
      enddo

      return
      end

      subroutine RVKEP (x,a,y,dyda,ma,ts,hkl)
      implicit none
      real*8 PI,TWOPI
      parameter (PI=3.14159265358979d0)
      parameter (TWOPI=2.0d0*PI)
      integer npl,ndset,idset,ma,i,j,NDSMAX,ts,hkl,gr_flag
      parameter (NDSMAX=20)
      integer idsmax(NDSMAX)
      real*8 x,y,a(ma),a2(ma),dyda(ma),mass(10),ap(10)
      real*8 cosw,sinw,capm,cape,cose,sine,cosf,sinf,fac1,fac2,fac3
      real*8 orbel_ehybrid, f, coswf,omega(10),capmm(10),ecc(10)
      real*8 ecc2,wm,sinwm,coswm,sin2wm,cos2wm,sin3wm,cos3wm,omegad(10)

      common /DSBLK/ npl,ndset,idsmax,idset,gr_flag

      y = 0.d0

      do i = 1,ma
          a2(i)=a(i)
      enddo

      if (hkl.eq.0) then

          do i = 1,npl
             j = 6*(i-1)

             if (a2(j+2).lt.0.d0) then  ! if P<0, set P>0
                a2(j+2) = dabs(a2(j+2))
             endif

             if (a2(j+1).lt.0.d0) then  ! if K<0, set K>0 and w = w+PI
                a2(j+4) = a2(j+4) + PI
                a2(j+1) = dabs(a2(j+1))
                if (a2(j+4).gt.2.d0*PI) a2(j+4) = a2(j+4)-2.d0*PI
             endif
             if (a2(j+3).lt.0.d0) then  ! if e<0, set e>0 and w=w+PI, M0=M0-PI
                a2(j+3) = dabs(a2(j+3))
                a2(j+4) = a2(j+4) +  PI
                if (a2(j+4).gt.2.d0*PI) a2(j+4) = a2(j+4)-2.d0*PI
                a2(j+5) = a2(j+5) - PI
                if (a2(j+5).lt.0.d0) a2(j+5) = a2(j+5)+2.d0*PI
             endif
             if (a2(j+4).lt.0.d0) a2(j+4)=dmod(a2(j+4)+2.d0*PI,2.d0*PI)
             if (a2(j+5).lt.0.d0) a2(j+5)=dmod(a2(j+5)+2.d0*PI,2.d0*PI)
c             if (a2(j+6).lt.0.d0) a2(j+6)=dmod(a2(j+6)+2.d0*PI,2.d0*PI)

             if (a2(j+4).gt.2.d0*PI) a2(j+4)=dmod(a2(j+4), 2.d0*PI)
             if (a2(j+5).gt.2.d0*PI) a2(j+5)=dmod(a2(j+5), 2.d0*PI)
c             if (a2(j+6).gt.2.d0*PI) a2(j+6)=dmod(a2(j+6), 2.d0*PI)

             ecc(i) = a2(j+3)
             omega(i) = a2(j+4)
             capmm(i) = a2(j+5)
c             omegad(i) = a2(j+6)

             if(gr_flag.ne.0) call MA_J (a,ma,npl,1.0d0,1.0d0,
     & 	           mass,ap,hkl,gr_flag)

             omegad(i) = a(j+6)

c             write(*,*) ecc(i) ,omega(i) ,capmm(i),omegad(i)

          enddo

      else

          do i = 1,npl
             j = 6*(i-1)
             if (a2(j+1).lt.0.d0) then  ! if K<0, set K>0 and w = w+PI
                a2(j+4) = -1.d0*a2(j+4)       !     which is h = -h, k = -k
                a2(j+3) = -1.d0*a2(j+3)
                a2(j+1) = dabs(a2(j+1))
             endif

             ecc(i) = dsqrt(a2(j+3)**2 + a2(j+4)**2)
             omega(i) = datan2(a2(j+3),a2(j+4))

             if(omega(i).lt.0.d0)omega(i)=dmod(omega(i)+2.d0*PI,2.d0*PI)
             if(omega(i).gt.0.d0)omega(i)=dmod(omega(i),        2.d0*PI)
             if (a2(j+5).lt.0.d0) a2(j+5)=dmod(a2(j+5)+2.d0*PI, 2.d0*PI)
             if (a2(j+5).gt.2.d0*PI) a2(j+5) = dmod(a2(j+5), 2.d0*PI)

             capmm(i) = a2(j+5) - omega(i)

             if(capmm(i).lt.0.d0)capmm(i)=dmod(capmm(i)+2.d0*PI,2.d0*PI)
             if(capmm(i).gt.0.d0)capmm(i)=dmod(capmm(i),        2.d0*PI)

c             write(*,*) a2(j+4),a2(j+4),ecc(i) ,omega(i) ,capmm(i)
          enddo

      endif

      if (hkl.eq.0) then
      do j = 1,npl

          i = 6*(j-1)
          cosw = dcos(omega(j)+omegad(j)*x/365.25d0)
          sinw = dsin(omega(j)+omegad(j)*x/365.25d0)

          capm = TWOPI*x/a2(2+i) + capmm(j)
          capm = dmod(capm,  2.d0*PI )

          cape = ORBEL_EHYBRID (ecc(j),capm)
          cose = dcos(cape)
          sine = dsin(cape)

          cosf = (cose - ecc(j))/(1.d0 - ecc(j)*cose)
          sinf = (dsqrt(1.d0 - ecc(j)**2)*sine)/(1.d0 - ecc(j)*cose)
c          f = 2.0d0*datan2( dsqrt(1.d0 - ecc(j))*dcos(cape/2.d0),
c     &                 dsqrt(1.d0 + ecc(j))*dsin(cape/2.d0))

c          coswf = dcos(omega(j)+f)
c          fac1 = coswf + ecc(j)*cosw

          fac1 = cosw*cosf - sinw*sinf + ecc(j)*cosw

          fac2 = (cosw*sinf + sinw*cosf)/(1.d0 - ecc(j)*cose)**2
          fac3 = -a2(1+i)*dsqrt(1.d0 - ecc(j)**2)*fac2

          y = y + a2(1+i)*fac1
          dyda(1+i) = fac1
          dyda(2+i) = -TWOPI*fac3*x/a2(2+i)**2
          dyda(3+i) = -a2(1+i)*sine*(2.d0 - ecc(j)**2 - ecc(j)*cose)*
     &                 fac2/dsqrt(1.d0 - ecc(j)**2)
          dyda(4+i) = -a2(1+i)*(sinw*cosf + cosw*sinf + ecc(j)*sinw)
          dyda(5+i) = fac3
          dyda(6+i) = -a2(1+i)*(sinw*cosf + cosw*sinf
     &                 + ecc(j)*sinw)*x/365.25d0
      enddo

      else

      do j = 1,npl

          i = 6*(j-1)
c          ecc2 = dsqrt(a2(3+i)**2 + a2(4+i)**2)

          if (ecc(j).gt.1.d-2) then
             cosw = dcos(omega(j)+omegad(i)*x/365.25d0)
             sinw = dsin(omega(j)+omegad(i)*x/365.25d0)

             capm = TWOPI*x/a2(2+i) + capmm(j)
             capm = dmod(capm,  2.d0*PI )
c             write(*,*) capm
             cape = ORBEL_EHYBRID (ecc(j),capm)
             cose = dcos(cape)
             sine = dsin(cape)
             cosf = (cose - ecc(j))/(1.d0 - ecc(j)*cose)
             sinf = (dsqrt(1.d0 - ecc(j)**2)*sine)/(1.d0 - ecc(j)*cose)

             fac1 = cosw*cosf - sinw*sinf + ecc(j)*cosw
             fac2 = cosw*sinf + sinw*cosf
             fac3 = -a2(1+i)*dsqrt(1.d0 - ecc(j)**2)*fac2/
     &              (1.d0 - ecc(j)*cose)**2

             y = y + a2(1+i)*fac1
             dyda(1+i) = fac1
             dyda(2+i) = -TWOPI*fac3*x/a2(2+i)**2
             dyda(3+i) = -a2(1+i)*fac2*((2.d0-ecc(j)**2-ecc(j)*cose)*
     &                   sinw*sine/dsqrt(1.d0-ecc(j)**2) -
     &                   dsqrt(1.d0 - ecc(j)*2)*cosw/ecc(j))/
     &                   (1.d0 - ecc(j)*cose)**2 -
     &                   a2(1+i)*fac2*cosw/ecc(j)
             dyda(4+i) = -a2(1+i)*fac2*((2.d0-ecc(j)**2-ecc(j)*cose)*
     &                   cosw*sine/dsqrt(1.d0 - ecc(j)**2) +
     &                   dsqrt(1.d0 - ecc(j)*2)*sinw/ecc(j))/
     &                   (1.d0 - ecc(j)*cose)**2 +
     &                   a2(1+i)*fac2*sinw/ecc(j)
             dyda(5+i) = fac3
             dyda(6+i) = dyda(4+i)*x/365.25d0

          else
             wm = TWOPI*x/a2(2+i) + a2(5+i)
             wm = dmod(wm,  2.d0*PI )

             coswm = dcos(wm)
             sinwm = dsin(wm)
             cos2wm = dcos(2.d0*wm)
             sin2wm = dsin(2.d0*wm)
             cos3wm = dcos(3.d0*wm)
             sin3wm = dsin(3.d0*wm)

             fac1 = coswm + a2(3+i)*sin2wm - a2(4+i)*(1.d0 - cos2wm) -
     &              a2(3+i)**2*(0.875d0*coswm + 1.125d0*cos3wm) -
     &              a2(3+i)*a2(4+i)*(0.25d0*sinwm - 2.25d0*sin3wm) -
     &              a2(4+i)**2*1.125d0*(coswm - cos3wm)
             fac3 = -sinwm + 2.d0*a2(3*i)*cos2wm - 2.d0*a2(4+i)*sin2wm +
     &              a2(3+i)**2*(0.875d0*sinwm - 3.375d0*sin3wm) -
     &              a2(3+i)*a2(4+i)*(0.25d0*coswm - 6.75d0*cos3wm) +
     &              a2(4+i)**2*(1.125d0*coswm - 3.375*sin3wm)

             y = y + a2(1+i)*fac1
             dyda(1+i) = fac1
             dyda(2+i) = -a2(1+i)*TWOPI*fac3*x/a2(2+i)**2
             dyda(3+i) = a2(1+i)*(sin2wm -
     &                 a2(3+i)*(1.75d0*coswm + 2.25d0*cos3wm) -
     &                 a2(4+i)*(0.25d0*sinwm - 2.25d0*sin3wm))
             dyda(4+i) = a2(1+i)*(-1.d0 + a2(4+i)*cos2wm -
     &              a2(3+i)**(0.25d0*sinwm - 2.25d0*sin3wm) -
     &              a2(4+i)*2.25d0*(coswm - cos3wm))
             dyda(5+i) = a2(1+i)*fac3
             dyda(6+i) = dyda(4+i)*x/365.25d0
          endif

      enddo

      endif

c      do i = 1,idset
      y = y + a2(6*npl+ts)
      dyda(6*npl+ts) = 1.d0
c      enddo

c      do i = 1,idset
c          y = y + a(5*npl+i)
c          dyda(5*npl+i) = 1.d0
c      enddo

      y = y + a2(6*npl +2*ndset + 1)*x
      y = y + a2(6*npl +2*ndset + 2) + a2(6*npl+2*ndset + 2)*x**2
c      write(*,*) a2(6*npl +2*ndset + 1), a2(6*npl+ndset + 2)
c      dyda(6*npl + ndset + 1) = x
c      dyda(6*npl + ndset + 2) = x**2

      do i = ts+1,ndset
          dyda(6*npl+i) = 0.d0
      enddo

      return
      end

      subroutine prepare_for_amoeba(p,mp,np,y,a,ia,ma,mfit,funk,ndata,
     & x,z,ymod,dyda,ts,sig, it,hkl)
      integer MMAX,NDSMAX,ma,ts(20000), ndata,mp,np,mfit,it
      parameter(MMAX=200,NDSMAX=20)
      REAL*8 ftol,p(mp,np),y(mp),a(MMAX), a2(mfit),fr,frjitt
      real*8 x(20000),z(20000),ymod(20000)
      real*8 dyda(MMAX), sig(20000), loglik
      parameter(fr=0.05, frjitt=0.05)
      INTEGER i,j,k, ia(MMAX), idsmax(NDSMAX),hkl,gr_flag
      external funk

      common /DSBLK/ npl,ndset,idsmax,idset,gr_flag

      k=0
      do j=1,ma
          if(ia(j).ne.0) then
          k=k+1
          p(1,k)=a(j)
          do i=2,mfit+1
              if (k.eq.(i-1)) then
                  if (j.gt.(6*npl+ndset)) then
                  p(i,k)=(1+frjitt)*(p(1,k)+0.1)
                  else
                  if (mod(j,5).eq.2) then
                     p(i,k)=(1+fr)*(p(1,k) + 0.1)
                  else if (mod(j,5).eq.3) then
                     p(i,k)=(1+frjitt)*(p(1,k)+0.1)
                  else
                     p(i,k)=(1+fr)*(p(1,k)+0.1)
                  endif
                  endif
              else
                  p(i,k)=p(1,k)
              endif
          enddo
          endif
      enddo
      do i=1,mfit+1
          do j=1,mfit
              a2(j)=p(i,j)
          enddo
          call funk(ndata,x,z,a2,ymod,dyda,ma,mfit,ts,sig,loglik,i,
     &              a,ia,hkl)
          y(i)=loglik
c          write(*,*) a2(1),a2(2),a2(3),a2(4),a2(5),a2(6),a2(7)
c          write(*,*) a2(8),a2(9),a2(10),a2(11),a2(12),a2(13),a2(14)
c          write(*,*) a2(15),a2(16),a2(17),a2(18),a2(19),a2(20),a2(21)
      enddo
      return
      end

      SUBROUTINE amoeba(p,y,mp,np,ndim,ftol,funk,iter,ndata,x,z,ymod,
     & dyda,ma,ts,sig,a,ia,ytry,hkl)
      implicit none
      INTEGER iter,mp,ndim,np,NMAX,ITMAX, MMAX,ma,ts(20000), ndata
      REAL*8 ftol,p(mp,np),y(mp),x(20000),z(20000),ymod(20000)
      PARAMETER (NMAX=20,ITMAX=200000,MMAX=200)
      real*8 dyda(MMAX), sig(20000), loglik, a(MMAX)
      EXTERNAL funk
      INTEGER i,ihi,ilo,inhi,j,m,n, ia(MMAX),hkl
      REAL*8 rtol,summ,swap,ysave,ytry,psum(ndim),amotry
      iter=0
1     do 12 n=1,ndim
        summ=0.d0
        do 11 m=1,ndim+1
          summ=summ+p(m,n)
11      continue
        psum(n)=summ
12    continue
2     ilo=1
      if (y(1).gt.y(2)) then
        ihi=1
        inhi=2
      else
        ihi=2
        inhi=1
      endif
      do 13 i=1,ndim+1
        if(y(i).le.y(ilo)) ilo=i
        if(y(i).gt.y(ihi)) then
          inhi=ihi
          ihi=i
        else if(y(i).gt.y(inhi)) then
          if(i.ne.ihi) inhi=i
        endif
13    continue
      rtol=2.d0*abs(y(ihi)-y(ilo))/(abs(y(ihi))+abs(y(ilo)))
      if (rtol.lt.ftol) then
        swap=y(1)
        y(1)=y(ilo)
        y(ilo)=swap
        do 14 n=1,ndim
          swap=p(1,n)
          p(1,n)=p(ilo,n)
          p(ilo,n)=swap
14      continue
        return
      endif
      if (iter.ge.ITMAX) then
          write (*,*) 'ITMAX exceeded in amoeba'
          return
      endif
      iter=iter+2
      ytry=amotry(p,y,psum,mp,np,ndim,funk,ihi,-1.0d0,ndata,x,z,ymod,
     & dyda,ma,ts,sig,a,ia,hkl)
      if (ytry.le.y(ilo)) then
        ytry=amotry(p,y,psum,mp,np,ndim,funk,ihi,2.0d0,ndata,x,z,ymod,
     & dyda,ma,ts,sig,a,ia,hkl)
      else if (ytry.ge.y(inhi)) then
        ysave=y(ihi)
        ytry=amotry(p,y,psum,mp,np,ndim,funk,ihi,0.5d0,ndata,x,z,ymod,
     & dyda,ma,ts,sig,a,ia,hkl)

        if (ytry.ge.ysave) then
          do 16 i=1,ndim+1
            if(i.ne.ilo)then
              do 15 j=1,ndim
                psum(j)=0.5d0*(p(i,j)+p(ilo,j))
                p(i,j)=psum(j)
15            continue
              call funk(ndata,x,z,psum,ymod,dyda,ma,ndim,ts,sig,loglik,
     &                 i,a,ia,hkl)
              y(i)=loglik
             endif
16        continue
          iter=iter+ndim
          goto 1
        endif
      else
        iter=iter-1
      endif
      goto 2

      END
C  (C) Copr. 1986-92 Numerical Recipes Software 0=M,173+9.

      FUNCTION amotry(p,y,psum,mp,np,ndim,funk,ihi,fac,ndata,x,z,ymod,
     & dyda,ma,ts,sig,a,ia,hkl)
      implicit none
      INTEGER ihi,mp,ndim,np,NMAX, MMAX, ma, ts(20000),ndata
      PARAMETER (NMAX=20, MMAX=200)
      REAL*8 amotry,fac,p(mp,np),psum(np),y(mp),x(20000),z(20000),
     & ymod(20000)
      real*8 dyda(MMAX), sig(20000),loglik
      EXTERNAL funk
      INTEGER j, ia(MMAX),hkl
      REAL*8 fac1,fac2,ytry,ptry(ndim), a(MMAX)
      fac1=(1.0d0-fac)/ndim
      fac2=fac1-fac
      do 11 j=1,ndim
        ptry(j)=psum(j)*fac1-p(ihi,j)*fac2
11    continue
      call funk(ndata,x,z,ptry,ymod,dyda,ma,ndim,ts,sig,loglik,ihi,
     & a,ia,hkl)

      ytry=loglik
C      WRITE(*,*) loglik
      if (ytry.lt.y(ihi)) then
        y(ihi)=ytry
        do 12 j=1,ndim
          psum(j)=psum(j)-p(ihi,j)+ptry(j)
          p(ihi,j)=ptry(j)
12      continue
      endif
      amotry=ytry
      return
      END
C  (C) Copr. 1986-92 Numerical Recipes Software 0=M,173+9.

	    subroutine MA_J (a,ma,npl,m0,sini,mass,ap,hkl,gr_flag)

	    implicit none
	    real*8 m0,PI,TWOPI,THIRD,GMSUN,dm,MSUN
        integer npl,ma,i,j,NPLMAX,hkl,gr_flag
        parameter (NPLMAX=7)
        real*8 sini,mm(NPLMAX),ecc,corr
        real*8 a(ma),mass(NPLMAX),ap(NPLMAX),mpold(NPLMAX),mtotal
	    parameter (THIRD=1.d0/3.d0)
        parameter (PI=3.14159265358979d0,TWOPI=2.d0*PI)
	    parameter (GMSUN=1.32712497d20,MSUN=1.32712497d20)

c*******G is set to be unit, and s, m, kg as unit of time, length and mass
c*******expectively.

        do j = 1,npl
           i = 6*(j-1)

           mm(j) = 2.d0*PI/(a(i+2)*8.64d4)

        enddo

        do i = 0,npl-1

           if (hkl.eq.0) then
               ecc = a(6*i+3)
           else
               ecc = dsqrt(a(6*i+3)**2+a(6*i+4)**2)    !! only for h, k
           endif

           mass(1) = m0
	       mpold(i+1) = 0.d0
 101       continue
           if (i.eq.0) then
           mtotal = m0
	       mass(i+2) = a(6*i+1)*(TWOPI/mm(i+1)*(m0 + mpold(i+1))**2/
     &               (TWOPI*GMSUN))**THIRD*
     &               dsqrt(1.d0 - ecc**2)
           else
              mtotal = m0
              do j = 0, i-1
                 mtotal = mtotal + mass(j+2)
              enddo
              mass(i+2) = a(6*i+1)*(TWOPI/mm(i+1)*(mtotal
     &                  +mpold(i+1))**2/(TWOPI*GMSUN))**THIRD*
     &                  dsqrt(1.d0 - ecc**2)
           endif

	    dm = dabs(mass(i+2)-mpold(i+1))/mass(i+2)
	    mpold(i+1) = mass(i+2)
           if (dm.gt.0) goto 101

	    ap(i+1) = (GMSUN*(mtotal + mass(i+2))*(1.d0/mm(i+1))
     &               **2)**THIRD

        enddo

        do i = 1,npl+1
           mass(i) = mass(i)*MSUN
        enddo

        if(gr_flag.ne.0) then
           do i = 1,npl
             j = 6*(i-1)
             call gr_corr(ap(i),a(j+3),mass(1),mass(i+1),i,corr,1.0d0)
             a(j+6) = corr*365.25
           enddo
        endif

	    return
	    end

        subroutine gr_corr(a,e,gmi,mass,n,corr,dt)
        implicit none

        real*8  a ,e, T, PI, c, GMSUN, AU, st_mass
        integer n,NPLMAX
      	real*8 mass, corr,THIRD,gmi,dt,prec_frac
        parameter (PI=3.14159265358979d0)
        parameter (c = 299792458.0d0)
	    parameter (GMSUN=1.32712497d20, AU=1.49597892d13)
        parameter (THIRD=1.d0/3.d0)

        T = 2.0d0*PI * sqrt((a**3.0d0)/(gmi ) )

         corr = (24.0d0 * (PI**3.0d0) * (a**2.0d0)) / ( (T**2.0d0) *
     &    (c**2.0d0)*(1.0d0-e**2.0d0) )

        corr = corr/T

        return
        end

***********************************************************************
c                    ORBEL_EHYBRID.F
***********************************************************************
*     PURPOSE:  Solves Kepler's eqn.   e is ecc.   m is mean anomaly.
*
*             Input:
*                           e ==> eccentricity anomaly. (real scalar)
*                           m ==> mean anomaly. (real scalar)
*             Returns:
*              orbel_ehybrid ==>  eccentric anomaly. (real scalar)
*
*     ALGORITHM: For e < 0.18 uses fast routine ESOLMD
*	         For larger e but less than 0.8, uses EGET
*	         For e > 0.8 uses EHIE
*     REMARKS: Only EHIE brings M and E into range (0,TWOPI)
*     AUTHOR: M. Duncan
*     DATE WRITTEN: May 25,1992.
*     REVISIONS: 2/26/93 hfl
***********************************************************************

	real*8 function orbel_ehybrid(e,m)

    include 'swift_loglik_Jakub.inc'

c...  Inputs Only:
	real*8 e,m

c...  Internals:
	real*8 orbel_esolmd,orbel_eget,orbel_ehie

c----
c...  Executable code

	if(e .lt. 0.18d0) then
	  orbel_ehybrid = orbel_esolmd(e,m)
	else
	  if( e .le. 0.8d0) then
	     orbel_ehybrid = orbel_eget(e,m)
	  else
	     orbel_ehybrid = orbel_ehie(e,m)
	  endif
	endif

	return
	end     ! orbel_ehybrid

c--------------------------------------------------------------------

***********************************************************************
c                    ORBEL_EHIE.F
***********************************************************************
*     PURPOSE:  Solves Kepler's eqn.   e is ecc.   m is mean anomaly.
*
*             Input:
*                           e ==> eccentricity anomaly. (real scalar)
*                           m ==> mean anomaly. (real scalar)
*             Returns:
*              orbel_ehybrid ==>  eccentric anomaly. (real scalar)
*
*     ALGORITHM: Use Danby's quartic for 3 iterations.
*                Eqn. is f(x) = x - e*sin(x+M). Note  that
*	         E = x + M. First guess is very good for e near 1.
*	         Need to first get M between 0. and PI and use
*		 symmetry to return right answer if M between PI and 2PI
*     REMARKS: Modifies M so that both E and M are in range (0,TWOPI)
*     AUTHOR: M. Duncan
*     DATE WRITTEN: May 25,1992.
*     REVISIONS:
***********************************************************************

      real*8 function orbel_ehie(e,m)

    include 'swift_loglik_Jakub.inc'

c...  Inputs Only:
	real*8 e,m

c...  Internals:
      integer iflag,nper,niter,NMAX
      real*8 dx,x,sa,ca,esa,eca,f,fp

      parameter (NMAX = 3)

c----
c...  Executable code

c In this section, bring M into the range (0,TWOPI) and if
c the result is greater than PI, solve for (TWOPI - M).
	iflag = 0
	nper = m/TWOPI
	m = m - nper*TWOPI
	if (m .lt. 0.d0) m = m + TWOPI

	if (m.gt.PI) then
	   m = TWOPI - m
	   iflag = 1
	endif

c Make a first guess that works well for e near 1.
	x = (6.d0*m)**(1.d0/3.d0) - m
	niter =0

c Iteration loop
	do niter =1,NMAX
	    call orbel_scget(x + m,sa,ca)
	    esa = e*sa
	    eca = e*ca
	    f = x - esa
	    fp = 1.d0 -eca
	    dx = -f/fp
	    dx = -f/(fp + 0.5d0*dx*esa)
	    dx = -f/(fp + 0.5d0*dx*(esa+0.3333333333333333d0*eca*dx))
	    x = x + dx
	enddo

	orbel_ehie = m + x

	if (iflag.eq.1) then
	  orbel_ehie = TWOPI - orbel_ehie
	  m = TWOPI - m
	endif

	return
	end         !orbel_ehie
c------------------------------------------------------------------

***********************************************************************
c	                  ORBEL_SCGET.F
***********************************************************************
*     PURPOSE:  Given an angle, efficiently compute sin and cos.
*
*        Input:
*             angle ==> angle in radians (real scalar)
*
*        Output:
*             sx    ==>  sin(angle)  (real scalar)
*             cx    ==>  cos(angle)  (real scalar)
*
*     ALGORITHM: Obvious from the code
*     REMARKS: The HP 700 series won't return correct answers for sin
*       and cos if the angle is bigger than 3e7. We first reduce it
*       to the range [0,2pi) and use the sqrt rather than cos (it's faster)
*       BE SURE THE ANGLE IS IN RADIANS - NOT DEGREES!
*     AUTHOR:  M. Duncan.
*     DATE WRITTEN:  May 6, 1992.
*     REVISIONS:
***********************************************************************

	subroutine orbel_scget(angle,sx,cx)

    include 'swift_loglik_Jakub.inc'

c...  Inputs Only:
        real*8 angle

c...  Output:
	real*8 sx,cx

c... Internals:
	integer nper
	real*8 x
	real*8 PI3BY2
	parameter(PI3BY2 = 1.5d0*PI)

c----
c...  Executable code

        nper = angle/TWOPI
	x = angle - nper*TWOPI
	if(x.lt.0.d0) then
           x = x + TWOPI
        endif
	sx = sin(x)
	cx= sqrt(1.d0 - sx*sx)
	if( (x .gt. PIBY2) .and. (x .lt.PI3BY2)) then
           cx = -cx
        endif

	return
	end   ! orbel_scget
c-------------------------------------------------------------------

***********************************************************************
c                    ORBEL_EGET.F
***********************************************************************
*     PURPOSE:  Solves Kepler's eqn.   e is ecc.   m is mean anomaly.
*
*             Input:
*                           e ==> eccentricity anomaly. (real scalar)
*                           m ==> mean anomaly. (real scalar)
*             Returns:
*                  orbel_eget ==>  eccentric anomaly. (real scalar)
*
*     ALGORITHM: Quartic convergence from Danby
*     REMARKS: For results very near roundoff, give it M between
*           0 and 2*pi. One can condition M before calling EGET
*           by calling my double precision function MOD2PI(M).
*           This is not done within the routine to speed it up
*           and because it works fine even for large M.
*     AUTHOR: M. Duncan
*     DATE WRITTEN: May 7, 1992.
*     REVISIONS: May 21, 1992.  Now have it go through EXACTLY two iterations
*                with the premise that it will only be called if
*	         we have an ellipse with e between 0.15 and 0.8
***********************************************************************

	real*8 function orbel_eget(e,m)

    include 'swift_loglik_Jakub.inc'

c...  Inputs Only:
	real*8 e,m

c...  Internals:
	real*8 x,sm,cm,sx,cx
	real*8 es,ec,f,fp,fpp,fppp,dx

c----
c...  Executable code

c Function to solve Kepler's eqn for E (here called
c x) for given e and M. returns value of x.
c MAY 21 : FOR e < 0.18 use ESOLMD for speed and sufficient accuracy
c MAY 21 : FOR e > 0.8 use EHIE - this one may not converge fast enough.

	  call orbel_scget(m,sm,cm)

c  begin with a guess accurate to order ecc**3
	  x = m + e*sm*( 1.d0 + e*( cm + e*( 1.d0 -1.5d0*sm*sm)))

c  Go through one iteration for improved estimate
	  call orbel_scget(x,sx,cx)
	  es = e*sx
	  ec = e*cx
	  f = x - es  - m
	  fp = 1.d0 - ec
	  fpp = es
	  fppp = ec
	  dx = -f/fp
	  dx = -f/(fp + dx*fpp/2.d0)
	  dx = -f/(fp + dx*fpp/2.d0 + dx*dx*fppp/6.d0)
	  orbel_eget = x + dx

c Do another iteration.
c For m between 0 and 2*pi this seems to be enough to
c get near roundoff error for eccentricities between 0 and 0.8

	  x = orbel_eget
	  call orbel_scget(x,sx,cx)
	  es = e*sx
	  ec = e*cx
	  f = x - es  - m
	  fp = 1.d0 - ec
	  fpp = es
	  fppp = ec
	  dx = -f/fp
	  dx = -f/(fp + dx*fpp/2.d0)
	  dx = -f/(fp + dx*fpp/2.d0 + dx*dx*fppp/6.d0)

	  orbel_eget = x + dx

	return
	end  ! orbel_eget
c---------------------------------------------------------------------

***********************************************************************
c                    ORBEL_ESOLMD.F
***********************************************************************
*     PURPOSE:  Solves Kepler's eqn.   e is ecc.   m is mean anomaly.
*
*             Input:
*                           e ==> eccentricity anomaly. (real scalar)
*                           m ==> mean anomaly. (real scalar)
*             Returns:
*                orbel_esolmd ==>  eccentric anomaly. (real scalar)
*
*     ALGORITHM: Some sort of quartic convergence from Wisdom.
*     REMARKS: ONLY GOOD FOR SMALL ECCENTRICITY SINCE IT ONLY
*         ITERATES ONCE. (GOOD FOR PLANET CALCS.)
*      	  ALSO DOES NOT PUT M OR E BETWEEN 0. AND 2*PI
*     INCLUDES: needs SCGET.F
*     AUTHOR: M. Duncan
*     DATE WRITTEN: May 7, 1992.
*     REVISIONS: 2/26/93 hfl
***********************************************************************

	real*8 function orbel_esolmd(e,m)

    include 'swift_loglik_Jakub.inc'

c...  Inputs Only:
      real*8 e,m

c...  Internals:
	real*8 x,sm,cm,sx,cx
	real*8 es,ec,f,fp,fpp,fppp,dx

c----
c...  Executable code

c...    Function to solve Kepler's eqn for E (here called
c...    x) for given e and M. returns value of x.

	  call orbel_scget(m,sm,cm)
	  x = m + e*sm*( 1.d0 + e*( cm + e*( 1.d0 -1.5d0*sm*sm)))

	  call orbel_scget(x,sx,cx)
	  es = e*sx
	  ec = e*cx
	  f = x - es  - m
	  fp = 1.d0 - ec
	  fpp = es
	  fppp = ec
	  dx = -f/fp
	  dx = -f/(fp + dx*fpp/2.d0)
	  dx = -f/(fp + dx*fpp/2.d0 + dx*dx*fppp/6.d0)

	  orbel_esolmd = x + dx

	return   ! orbel_esolmd
	end
c--------------------------------------------------------------------

max_stars_repo_name: Simske/exostriker


text:
      SUBROUTINE CONVRT
C
C        'CONVRT' CONVERTS INPUT DATA INTO THE FORM REQUIRED
C        FOR INTERNAL USE
C
      IMPLICIT REAL*8 (A-H,O-Z)
C
      REAL*8 I3,I2OVI3
C
      COMMON/CCMBNZ/ CMRK(10),Z01(4),Z02(4),Z03(4),Z04(4),Z11(12),
     .               Z12(12),Z13(12),Z14(12),Z15(12),Z16(12),Z21(36),
     .               Z22(36),Z23(36),Z24(36),Z25(36),Z26(36),Z27(36),
     .               Z28(36),Z31(108),Z32(108),Z33(108),Z34(108),
     .               Z35(108),Z41(324),Z42(324),Z43(324),ZK21(18),
     .               ZK22(18),ZK23(18),ZK31(54),ZK32(54),ZK33(54),
     .               ZK34(54),ZK35(54),ZK36(54),ZK41(162),ZK42(162),
     .               ZK43(162),ZK44(162),ZK45(162),ZK46(162),ZK47(162),
     .               ZK48(162),ZS01(2),ZS02(2),ZS11(6),ZS12(6),ZS13(6),
     .               ZS21(18),ZS22(18),ZS23(18),ZS24(18),ZS31(54),
     .               ZS32(54),ZS41(162),ZT21(108),ZT22(108),STAO(18),
     .               STA1(18),STA2(18),STBO(18),STB1(18),STB2(18),
     .               Z2S01(2),Z2S12(6),Z2S23(18)
C
      COMMON/CCNVRT/ BDYMI(3,3),DPRMI(3,3),EMODLS(10),RTUBE(10),
     .               HTUBE(10),THERMC(10),TIPMS(10),C(10)
C
      COMMON/CFNALP/ SCO,CIZO(3,3),CIYO(3,3),DONA(3,10),DONB(3,10)
     .              ,CDAMP(3,10)
C
      COMMON/CONSTS/ PI,TWOPI,RADIAN
C
      COMMON/CSOLAR/ SAO(10),SKA(9),SKB(9),SKOA(10,3),SKOB(10,3),
     .               STMK(10),SKAA(10,9),SKBB(10,9)
C
      COMMON/CSTVAL/ TSTART,ZL0(10),ZL1(10),ZL2(10),ZLA(10)
C
      COMMON/EISUBK/ EI(10)
C
      COMMON/ELKDMP/ OMKDMP(3,10),IOMKDM(10)
C
      COMMON/IPOOL1/ IGRAV,IDAMP,IK,K1,ITIM,IAB,IAPS,IBB,IBPS,NK(10),
     .               LK(10),LLK(10)
C
C
      COMMON/RPOOL1/ RHOK(10),TIME,SA(3,3),FM1(3,3),ZLK(10),OMEG(3),
     .               ZLKP(10),ZLKDP(10),CMAT(3,3),GBAR(3,3),YBCM(3),
     .               ZBZK(3,10),FCM(3,3),DTO,PHID,PHI
C
C
      COMMON/RNEWR/ZA(10),I3(10),I2OVI3(10),
     *             ZDQ(10),ZJ(10)
C
      COMMON/RTDIST/ TDIS(10)
C
      COMMON/SATLSH/ TAUK(10),OCULTK(10),RADSH
C
C
      DIMENSION STAUK(10),SOCULT(10)
      DIMENSION RHO(10),SAOK(10),ZLOK(10),ZL1K(10),ZL2K(10),EM(10),
     .          HT(10),RT(10),TI(10),TH(10),ST(10),DONB1(3,10),PROD(10),
     .          IL(10),ILL(10),SKOA1(10,3),SKOB1(10,3),TDIST(10)
C
C
C
      REAL*8 INERTA
C
C
C
      DO 5 I=1,3
      DO 5 J=1,3
      CIYO(I,J)=0.D0
    5 CIZO(I,J)=0.D0
C
C
      CIYO(1,1)=(BDYMI(3,3) + BDYMI(2,2) - BDYMI(1,1))/2.D0
      CIYO(2,2)=(BDYMI(3,3) + BDYMI(1,1) - BDYMI(2,2))/2.D0
      CIYO(3,3)=(BDYMI(1,1) + BDYMI(2,2) - BDYMI(3,3))/2.D0
C
      CIZO(1,1)=(DPRMI(3,3) + DPRMI(2,2) - DPRMI(1,1))/2.D0
      CIZO(2,2)=(DPRMI(3,3) + DPRMI(1,1) - DPRMI(2,2))/2.D0
      CIZO(3,3)=(DPRMI(1,1) + DPRMI(2,2) - DPRMI(3,3))/2.D0
C
      CIYO(1,2)=BDYMI(1,2)
      CIYO(1,3)=BDYMI(1,3)
      CIYO(2,3)=BDYMI(2,3)
C
      CIZO(1,2)=DPRMI(1,2)
      CIZO(1,3)=DPRMI(1,3)
      CIZO(2,3)=DPRMI(2,3)
C
      DO 10 I=1,3
      DO 10 J=I,3
      CIYO(J,I)=CIYO(I,J)
   10 CIZO(J,I)=CIZO(I,J)
C
      DO 20 K=1,IK
      STAUK(K)=TAUK(K)
      SOCULT(K)=OCULTK(K)
      TDIST(K)=TDIS(K)
      RHO(K)=RHOK(K)
      SAOK(K)=SAO(K)
      ZLOK(K)=ZL0(K)
      ZL1K(K)=ZL1(K)
      ZL2K(K)=ZLA(K)
      EM(K)=EMODLS(K)
      HT(K)=HTUBE(K)
      RT(K)=RTUBE(K)
      TI(K)=TIPMS(K)
      TH(K)=THERMC(K)
      ST(K)=STMK(K)
      IL(K)=LK(K)
      ILL(K)=LLK(K)
C
      DO 20 I=1,3
      SKOB1(K,I)=SKOB(K,I)
      SKOA1(K,I)=SKOA(K,I)
   20 DONB1(I,K)=CDAMP(I,K)
C
C
C
      DO 60 K=1,IK
      L=IK-K1+K
      IF((K-K1).GT.0) L=K-K1
C
      TAUK(K)=STAUK(L)
      OCULTK(K)=SOCULT(L)
      TDIS(K)=TDIST(L)
      RHOK(K)=RHO(L)
      SAO(K)=SAOK(L)
      ZL0(K)=ZLOK(L)
      ZL1(K)=ZL1K(L)
      ZLA(K)=ZL2K(L)
      EMODLS(K)=EM(L)
      HTUBE(K)=HT(L)
      RTUBE(K)=RT(L)
      TIPMS(K)=TI(L)
      THERMC(K)=TH(L)
      STMK(K)=ST(L)
      LK(K)=IL(L)
      LLK(K)=ILL(L)
C
      DO 30 I=1,3
      SKOA(K,I)=SKOA1(L,I)
      SKOB(K,I)=SKOB1(L,I)
   30 CDAMP(I,K)=DONB1(I,L)
C
   60 CONTINUE
C
C
      DO 70 I=1,IK
      INERTA=PI*RTUBE(I)**3*HTUBE(I)
      PROD(I)=EMODLS(I)*INERTA/144.D0
      EI(I)=6.D0*PROD(I)*THERMC(I)/RTUBE(I)*TDIS(I)
      CMRK(I)=TIPMS(I)/RHOK(I)
      DO 70 J=1,9
      SKAA(I,J)=PROD(I)*(1.0D0+I2OVI3(I))
   70 SKBB(I,J)=PROD(I)*(1.0D0+I2OVI3(I))
C
C
      DO 80 I=1,9
      SKA(I)=0.D0
   80 SKB(I)=0.D0
C
      DO 90 I=1,IK
      SUM=DSQRT(PROD(I)*RHOK(I))
      DO 90 J=1,3
      GO TO (101,102,103),J
  101 N1=1
      N2=19
      M1=1
      M2=10
      GO TO 104
  102 N1=5
      N2=23
      M1=5
      M2=14
      GO TO 104
  103 N1=9
      N2=27
      M1=9
      M2=18
  104 N=N1
      IF (LK(I).EQ.2) N=N2
      M=M1
      IF (LLK(I).EQ.2) M=M2
      SUM1=DSQRT(2.0D0*Z21(N)*ZK21(M))
      SUM2=DSQRT(2.0D0*Z21(N)*ZK23(M))
      DONA(J,I)=SUM*SUM1*2.0D0*CDAMP(J,I)
      DONB(J,I)=SUM*SUM2*2.0D0*CDAMP(J,I)
      IF(IOMKDM(I).EQ.0) GO TO 90
      DONA(J,I)=2.0D0*CDAMP(J,I)*OMKDMP(J,I)
      DONB(J,I)=2.0D0*CDAMP(J,I)*OMKDMP(J,I)
   90 CONTINUE
C
      RETURN
      END

ext: for