[SIZE=0]
c Last revision: 9/25/12 S.R.
subroutine io_init_tp_theia(ntp,npl,x_earth,y_earth,z_earth,
& vx_earth,vy_earth,vz_earth,
& xht,yht,zht,vxht,vyht,vzht,istat,rstat,frame)
USE IFPORT
include '../scatr.inc'
include 'io.inc'
c... Input
integer ntp,npl
real*8 x_earth,y_earth,z_earth
real*8 vx_earth,vy_earth,vz_earth
character*80 frame
c... Output
real*8 xht(NTPMAX),yht(NTPMAX),zht(NTPMAX)
real*8 vxht(NTPMAX),vyht(NTPMAX),vzht(NTPMAX)
real*8 rstat(NTPMAX,NSTATR)
integer istat(NTPMAX,NSTAT)
c... Internal
integer iseed
real*8 rran
integer i,j,ierr,ns,ftest
real*8 v_esc,r_earth
real*8 y,phi,vn,y_min,y_max
real*8 nx,ny,nz
integer SCATTERING_ANGLE
c-----
c... Executable code
c Initialize
c 1 AU=149597870691 m
c Initialize random generator
iseed = 35791246
call srand(iseed)
SCATTERING_ANGLE = 20
v_esc = 6.4604555902823027300785018765024e-3 ! AU/day = 11186.0 m/s
r_earth = 4.2563891254523551325458216979349e-5; ! mean radius of earth in AU (=6367467.5 m)
y_min = 0.5 * (1.0 - cos((90.0 - SCATTERING_ANGLE)/PI*180.0))
y_max = 0.5 * (1.0 - cos((90.0 + SCATTERING_ANGLE)/PI*180.0))
if(ntp.gt.NTPMAX) then
write(*,*) ' SWIFT ERROR: in io_init_tp_theia: '
write(*,*) ' The number of test bodies,',ntp,','
write(*,*) ' is too large, it must be less than',NTPMAX
call util_exit(1)
endif
write(*,*) ' '
write(*,*) 'ntp : ',ntp
if(ntp.eq.0) then
write(*,*) ' '
return
endif ! <===== NOTE
c... Determine the number of istat and rstat variables. In what follows,
c... we assume that they are the same.
ns = npl + 2
if(ns.ne.NSTAT) then
write(*,*) 'Warning: The size of istat and rstat arrays is '
write(*,*) ' not NSTAT=',NSTAT,', but is ',ns
endif
c Read in the x's and v's and istat(*,*)
write(*,*) ' '
do i=1,ntp
rran = rand()
y = rran * (y_max - y_min) + y_min ! add random number instead of 0.5d0
rran = rand()
phi = rran * TWOPI ! add random number instead of 0.5d0
rran = rand()
vn = rran * 0.4d0 + 1.0d0 ! add random number instead of 0.5d0
nx = 2.0d0 * sqrt(y - y * y) * cos(phi)
ny = 2.0d0 * sqrt(y - y * y) * sin(phi)
nz = 1.0d0 - 2.0d0 * y
xht(i) = -9.814911942588033e-001 + nx * r_earth
yht(i) = -1.866158863169087e-001 + ny * r_earth
zht(i) = 5.097029871362697e-006 + nz * r_earth
vxht(i) = 2.931839043977104e-003 + vn * v_esc * nx
vyht(i) = -1.697334005001157e-002 + vn * v_esc * ny
vzht(i) = 3.022067795326298e-007 + vn * v_esc * nz
do j=1,ns
istat(i,j) = 0
enddo
do j=ns+1,NSTAT
istat(i,j) = 0
if (j.eq.NSTAT-2.and.
& (frame(1:3).eq.'hel'.or.frame(1:3).eq.'HEL')) then
istat(i,j) = 1
endif
enddo
do j=1,NSTATR
rstat(i,j) = 0.0d0
enddo
enddo
return
end ! io_init_tp_theia.f
c-----------------------------------------------------------------[/SIZE]