Newer
Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
c
c ELEMENT6.FOR (ErikSoft 5 June 2001)
c
c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
c
c Author: John E. Chambers
c
c Makes output files containing Keplerian orbital elements from data created
c by Mercury6 and higher.
c
c The user specifies the names of the required objects in the file elements.in
c See subroutine M_FORMAT for the identities of each element in the EL array
c e.g. el(1)=a, el(2)=e etc.
c
c------------------------------------------------------------------------------
c
implicit none
include 'mercury.inc'
c
integer itmp,i,j,k,l,iback(NMAX),precision,lenin
integer nmaster,nopen,nwait,nbig,nsml,nbod,nsub,lim(2,100)
integer year,month,timestyle,line_num,lenhead,lmem(NMESS)
integer nchar,algor,centre,allflag,firstflag,ninfile,nel,iel(22)
integer nbod1,nbig1,unit(NMAX),code(NMAX),master_unit(NMAX)
real*8 time,teval,t0,t1,tprevious,rmax,rcen,rfac,rhocgs,temp
real*8 mcen,jcen(3),el(22,NMAX),s(3),is(NMAX),ns(NMAX),a(NMAX)
real*8 mio_c2re, mio_c2fl,fr,theta,phi,fv,vtheta,vphi,gm
real*8 x(3,NMAX),v(3,NMAX),xh(3,NMAX),vh(3,NMAX),m(NMAX)
logical test
character*250 string,fout,header,infile(50)
character*80 mem(NMESS),cc,c(NMAX)
character*8 master_id(NMAX),id(NMAX)
character*5 fin
character*1 check,style,type,c1
character*2 c2
c
c------------------------------------------------------------------------------
c
allflag = 0
tprevious = 0.d0
rhocgs = AU * AU * AU * K2 / MSUN
c
c Read in output messages
inquire (file='message.in', exist=test)
if (.not.test) then
write (*,'(/,2a)') ' ERROR: This file is needed to continue: ',
% ' message.in'
stop
end if
open (14, file='message.in', status='old')
10 continue
read (14,'(i3,1x,i2,1x,a80)',end=20) j,lmem(j),mem(j)
goto 10
20 close (14)
c
c Open file containing parameters for this programme
inquire (file='element.in', exist=test)
if (test) then
open (10, file='element.in', status='old')
else
call mio_err (6,mem(81),lmem(81),mem(88),lmem(88),' ',1,
% 'element.in',9)
end if
c
c Read number of input files
30 read (10,'(a250)') string
if (string(1:1).eq.')') goto 30
call mio_spl (250,string,nsub,lim)
read (string(lim(1,nsub):lim(2,nsub)),*) ninfile
c
c Make sure all the input files exist
do j = 1, ninfile
40 read (10,'(a250)') string
if (string(1:1).eq.')') goto 40
call mio_spl (250,string,nsub,lim)
infile(j)(1:(lim(2,1)-lim(1,1)+1)) = string(lim(1,1):lim(2,1))
inquire (file=infile(j), exist=test)
if (.not.test) call mio_err (6,mem(81),lmem(81),mem(88),
% lmem(88),' ',1,infile(j),80)
end do
c
c What type elements does the user want?
centre = 0
45 read (10,'(a250)') string
if (string(1:1).eq.')') goto 45
call mio_spl (250,string,nsub,lim)
c2 = string(lim(1,nsub):(lim(1,nsub)+1))
if (c2.eq.'ce'.or.c2.eq.'CE'.or.c2.eq.'Ce') then
centre = 0
else if (c2.eq.'ba'.or.c2.eq.'BA'.or.c2.eq.'Ba') then
centre = 1
else if (c2.eq.'ja'.or.c2.eq.'JA'.or.c2.eq.'Ja') then
centre = 2
else
call mio_err (6,mem(81),lmem(81),mem(107),lmem(107),' ',1,
% ' Check element.in',23)
end if
c
c Read parameters used by this programme
timestyle = 1
do j = 1, 4
50 read (10,'(a250)') string
if (string(1:1).eq.')') goto 50
call mio_spl (250,string,nsub,lim)
c1 = string(lim(1,nsub):lim(2,nsub))
if (j.eq.1) read (string(lim(1,nsub):lim(2,nsub)),*) teval
teval = abs(teval) * .999d0
if (j.eq.2.and.(c1.eq.'d'.or.c1.eq.'D')) timestyle = 0
if (j.eq.3.and.(c1.eq.'y'.or.c1.eq.'Y')) timestyle = timestyle+2
if (j.eq.4) call m_format (string,timestyle,nel,iel,fout,header,
% lenhead)
end do
c
c Read in the names of the objects for which orbital elements are required
nopen = 0
nwait = 0
nmaster = 0
60 continue
read (10,'(a250)',end=70) string
call mio_spl (250,string,nsub,lim)
if (string(1:1).eq.')'.or.lim(1,1).eq.-1) goto 60
c
c Either open an aei file for this object or put it on the waiting list
nmaster = nmaster + 1
itmp = min(7,lim(2,1)-lim(1,1))
master_id(nmaster)=' '
master_id(nmaster)(1:itmp+1) = string(lim(1,1):lim(1,1)+itmp)
if (nopen.lt.NFILES) then
nopen = nopen + 1
master_unit(nmaster) = 10 + nopen
call mio_aei (master_id(nmaster),'.aei',master_unit(nmaster),
% header,lenhead,mem,lmem)
else
nwait = nwait + 1
master_unit(nmaster) = -2
end if
goto 60
c
70 continue
c If no objects are listed in ELEMENT.IN assume that all objects are required
if (nopen.eq.0) allflag = 1
close (10)
c
c------------------------------------------------------------------------------
c
c LOOP OVER EACH INPUT FILE CONTAINING INTEGRATION DATA
c
90 continue
firstflag = 0
do i = 1, ninfile
line_num = 0
open (10, file=infile(i), status='old')
c
c Loop over each time slice
100 continue
line_num = line_num + 1
read (10,'(3a1)',end=900,err=666) check,style,type
line_num = line_num - 1
backspace 10
c
c Check if this is an old style input file
if (ichar(check).eq.12.and.(style.eq.'0'.or.style.eq.'1'.or.
% style.eq.'2'.or.style.eq.'3'.or.style.eq.'4')) then
write (*,'(/,2a)') ' ERROR: This is an old style data file',
% ' Try running m_elem5.for instead.'
stop
end if
if (ichar(check).ne.12) goto 666
c
c------------------------------------------------------------------------------
c
c IF SPECIAL INPUT, READ TIME, PARAMETERS, NAMES, MASSES ETC.
c
if (type.eq.'a') then
line_num = line_num + 1
read (10,'(3x,i2,a62,i1)') algor,cc(1:62),precision
c
c Decompress the time, number of objects, central mass and J components etc.
time = mio_c2fl (cc(1:8))
nbig = int(.5d0 + mio_c2re(cc(9:16), 0.d0, 11239424.d0, 3))
nsml = int(.5d0 + mio_c2re(cc(12:19),0.d0, 11239424.d0, 3))
mcen = mio_c2fl (cc(15:22))
jcen(1) = mio_c2fl (cc(23:30))
jcen(2) = mio_c2fl (cc(31:38))
jcen(3) = mio_c2fl (cc(39:46))
rcen = mio_c2fl (cc(47:54))
rmax = mio_c2fl (cc(55:62))
rfac = log10 (rmax / rcen)
c
c Read in strings containing compressed data for each object
do j = 1, nbig + nsml
line_num = line_num + 1
read (10,'(a)',err=666) c(j)(1:51)
end do
c
c Create input format list
if (precision.eq.1) nchar = 2
if (precision.eq.2) nchar = 4
if (precision.eq.3) nchar = 7
Loading full blame...