aboutsummaryrefslogtreecommitdiff
path: root/Binary.f
blob: 889f8bac223ca3fd1c4be0231744f96ccfd487c2 (plain) (blame)
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

      subroutine binary
c******************************************************************************
c     This program synthesizes a section of a spectroscopic binary star
c     spectrum, using one model atmosphere for the primary and one model
c     for the secondary, and then compares it to an observed spectrum.
c******************************************************************************

      implicit real*8 (a-h,o-z)
      include 'Atmos.com'
      include 'Factor.com'
      include 'Mol.com'
      include 'Linex.com'
      include 'Pstuff.com'
      include 'Equivs.com'
      include 'Multistar.com'
      integer ncall


c*****examine the parameter file
      begin = 0
      ncall = 1
1     call params
      linprintopt = linprintalt
      if (begin == 0) then
         if (numpecatom > 0) then
            do i=3,95
               binpec(syncount,i) = pec(i)
               do j=1,numatomsyn
                  binpecabund(syncount,i,j) = pecabund(i,j)
               enddo
            enddo
         endif
      else
         if (numpecatom > 0) then
            do i=3,95
               pec(i) = binpec(syncount,i)
               do j=1,numatomsyn
                  pecabund(i,j) = binpecabund(syncount,i,j)
               enddo
            enddo
         endif
      endif


c*****open the files for: standard output, raw spectrum depths, smoothed 
c     spectra, and (if desired) IRAF-style smoothed spectra
      istat = ivcleof(5,1)
      nf1out = 20     
      lscreen = 4
      array = 'STANDARD OUTPUT'
      nchars = 15
      call infile ('output ',nf1out,'formatted  ',0,nchars,
     .             f1out,lscreen)
      nf2out = 21               
      lscreen = lscreen + 2
      array = 'RAW SYNTHESIS OUTPUT'
      nchars = 20
      call infile ('output ',nf2out,'formatted  ',0,nchars,
     .             f2out,lscreen)
      if (syncount == 1) then
         f7out = f2out
      else
         f8out = f2out
      endif
      if (iraf /= 0) then
         nf4out = 23               
         lscreen = lscreen + 2
         array = 'IRAF ("rtext") OUTPUT'
         nchars = 24
         call infile ('output ',nf4out,'formatted  ',0,nchars,
     .                f4out,lscreen)
      endif


c*****open and read the model atmosphere file
      nfmodel = 30
      lscreen = lscreen + 2
      array = 'THE MODEL ATMOSPHERE'
      nchars = 20
      call infile ('input  ',nfmodel,'formatted  ',0,nchars,
     .             fmodel,lscreen)
      call inmodel


c*****open the line list file and the strong line list file
      nflines = 31
      lscreen = lscreen + 2
      array = 'THE LINE LIST'
      nchars = 13
      call infile ('input  ',nflines,'formatted  ',0,nchars,
     .              flines,lscreen)
      if (dostrong > 0) then
         nfslines = 32
         lscreen = lscreen + 2
         array = 'THE STRONG LINE LIST'
         nchars = 20
         call infile ('input  ',nfslines,'formatted  ',0,nchars,
     .                 fslines,lscreen)
      endif
      

c*****do the syntheses
      if (numpecatom == 0 .or. numatomsyn == 0) then
         isorun = 1
         nlines = 0
         mode = 3
         call inlines (1)
         call eqlib
         call nearly (1)
         call synspec
      else
         do n=1,numatomsyn
            isynth = n
            isorun = 1
            start = oldstart
            sstop = oldstop
            mode = 3
            call inlines (1)
            call eqlib
            call nearly (1)
            call synspec
            linprintopt = 0
         enddo
      endif
      if (syncount == 1) then
         fluxprimary = flux
      else
         fluxsecondary = flux
      endif
         

c*****finish the syntheses
      call finish (1)
      istat = ivcleof(4,1)
      if (control /= 'gridend') go to 1


c*****combine the synthetic spectra for plotting
      call binplotprep


c*****now plot the spectrum, maybe iterating abundances, and end the program
      if (plotopt==2 .and. specfileopt>0) then
         nfobs = 33               
         lscreen = lscreen + 2
         array = 'THE OBSERVED SPECTRUM'
         nchars = 21
         if (specfileopt==1 .or. specfileopt==3) then
            call infile ('input  ',nfobs,'unformatted',2880,nchars,
     .                   fobs,lscreen)
         else
            call infile ('input  ',nfobs,'formatted  ',0,nchars,
     .                   fobs,lscreen)
         endif
      endif
      if (plotopt /= 0) then
         control = 'binary '
         nf2out = nf9out
         nf3out = nf10out
         call pltspec (lscreen,ncall)
         if (choice == 'n') then
            do syncount=1,2
               if (numpecatom > 0) then
                  do i=3,95
                     pec(i) = binpec(syncount,i)
                     do j=1,numatomsyn
                        pecabund(i,j) = binpecabund(syncount,i,j)
                     enddo
                  enddo
               endif
               call chabund
               do i=3,95
                  binpec(syncount,i) = pec(i)
                  do j=1,numatomsyn
                     binpecabund(syncount,i,j) = pecabund(i,j)
                  enddo
               enddo
            enddo
            call finish (1)
            linecount = 0
            oldcount = 0
            begin = begin + 1
            choice = '1'
            ncall = 2
            go to 1
         endif
      endif

      call finish (0)

      end