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
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
|
!Copyright (C) 2013- Roger Wesson
!Free under the terms of the GNU General Public License v3
module mod_commandline
use mod_functions
use mod_globals
contains
subroutine readcommandline(redshiftguess_initial,resolutionguess_initial,normalisation)
implicit none
character(len=512), dimension(:), allocatable :: options
integer :: Narg,nargused,i,exclusioncount
real :: excludewavelength,redshiftguess_initial,resolutionguess_initial,normalisation
#ifdef CO
print *,"subroutine: readcommandline"
#endif
c=299792.458 !km/s
spectrumfile=""
narg = 0
nargused = 0 !to count options specified
narg = IARGC() !count input arguments
exclusioncount = 0 !to count lines excluded
subtractcontinuum = .true.
if (narg .eq. 0) then
print *,"Usage: alfa [options] [file]"
print *," [file] is an ascii file with columns for wavelength and flux"
print *," or a FITS file with 1, 2 or 3 dimensions, containing spectra."
print *," see the man page or online documentation for details of the options"
call exit(0)
endif
call get_command(commandline)
allocate (options(Narg))
options=""
print *,gettime(),"command line: ",trim(commandline)
! read command line options into array, counting how many times the exclude line option is present
do i=1,Narg
call get_command_argument(i,options(i))
if (trim(options(i)).eq."-el" .or. trim(options(i)).eq."--exclude-line") then
exclusioncount = exclusioncount + 1
endif
enddo
allocate(exclusions(exclusioncount))
if (exclusioncount .gt. 0) then
exclusioncount = 1 !now repurposing this variable to be an index for the array
else
deallocate(exclusions)
endif
! process the options
do i=1,narg
if ((trim(options(i))=="-n" .or. trim(options(i))=="--normalise")) then
if ((i+1) .le. Narg) then
read (options(i+1),*) normalisation
normalise=.true.
options(i:i+1)=""
else
print *,gettime(),"[100] no value specified for ",trim(options(i))
call exit(100)
endif
endif
if ((trim(options(i))=="-vg" .or. trim(options(i))=="--velocity-guess")) then
if ((i+1) .le. Narg) then
read (options(i+1),*) redshiftguess_initial
options(i:i+1)=""
else
print *,gettime(),"[100] no value specified for ",trim(options(i))
call exit(100)
endif
endif
if ((trim(options(i))=="-rg" .or. trim(options(i))=="--resolution-guess")) then
if ((i+1) .le. Narg) then
read (options(i+1),*) resolutionguess_initial
resolution_estimated=.true.
options(i:i+1)=""
if (resolutionguess_initial .lt. 0.) then
print *,gettime(),"[100] invalid value given for resolution guess"
call exit(100)
endif
else
print *,gettime(),"[100] no value specified for ",trim(options(i))
call exit(100)
endif
endif
if ((trim(options(i))=="-vtol1" .or. trim(options(i))=="--velocity-tolerance-1")) then
if ((i+1) .le. Narg) then
read (options(i+1),*) vtol1
vtol1 = vtol1/c
options(i:i+1)=""
if (vtol1 .lt. 0.) then
print *,gettime(),"[100] invalid value given for vtol1"
call exit(100)
endif
else
print *,gettime(),"[100] no value specified for ",trim(options(i))
call exit(100)
endif
endif
if ((trim(options(i))=="-vtol2" .or. trim(options(i))=="--velocity-tolerance-2")) then
if ((i+1) .le. Narg) then
read (options(i+1),*) vtol2
vtol2 = vtol2/c
options(i:i+1)=""
if (vtol2 .lt. 0.) then
print *,gettime(),"[100] invalid value given for vtol2"
call exit(100)
endif
else
print *,gettime(),"[100] no value specified for ",trim(options(i))
call exit(100)
endif
endif
if ((trim(options(i))=="-rtol1" .or. trim(options(i))=="--resolution-tolerance-1")) then
if ((i+1) .le. Narg) then
read (options(i+1),*) rtol1
options(i:i+1)=""
if (rtol1 .lt. 0.) then
print *,gettime(),"[100] invalid value given for rtol1"
call exit(100)
endif
else
print *,gettime(),"[100] no value specified for ",trim(options(i))
call exit(100)
endif
endif
if ((trim(options(i))=="-rtol2" .or. trim(options(i))=="--resolution-tolerance-2")) then
if ((i+1) .le. Narg) then
read (options(i+1),*) rtol2
options(i:i+1)=""
if (rtol1 .lt. 0.) then
print *,gettime(),"[100] invalid value given for rtol1"
call exit(100)
endif
else
print *,gettime(),"[100] no value specified for ",trim(options(i))
call exit(100)
endif
endif
if (trim(options(i))=="-ss" .or. trim(options(i))=="--subtract-sky") then
subtractsky=.true.
options(i)=""
endif
if (trim(options(i))=="-b" .or. trim(options(i))=="--bad-data") then
if ((i+1) .le. Narg) then
read (options(i+1),*) baddata
options(i:i+1)=""
else
print *,gettime(),"[100] no value specified for ",trim(options(i))
call exit(100)
endif
endif
if ((trim(options(i))=="-o" .or. trim(options(i))=="--output-dir")) then
if ((i+1) .le. Narg) then
read (options(i+1),"(A)") outputdirectory
outputdirectory=trim(outputdirectory)//"/"
inquire(file=trim(outputdirectory), exist=file_exists) ! trailing slash ensures it's looking for a directory
else
print *,gettime(),"[100] no value specified for ",trim(options(i))
call exit(100)
endif
if (.not. file_exists) then
print *,gettime(),"[102] output directory does not exist"
call exit(102)
endif
options(i:i+1)=""
endif
if (trim(options(i))=="--sky-catalogue" .or. trim(options(i))=="-skyc" .or. trim(options(i))=="-skycat") then
if ((i+1) .le. Narg) then
read (options(i+1),"(A)") skylinelistfile
options(i:i+1)=""
else
print *,gettime(),"[100] no value specified for ",trim(options(i))
call exit(100)
endif
endif
if (trim(options(i))=="--strong-catalogue" .or. trim(options(i))=="-sc" .or. trim(options(i))=="-strongcat") then
if ((i+1) .le. Narg) then
read (options(i+1),"(A)") stronglinelistfile
options(i:i+1)=""
else
print *,gettime(),"[100] no value specified for ",trim(options(i))
call exit(100)
endif
endif
if (trim(options(i))=="--deep-catalogue" .or. trim(options(i))=="-dc" .or. trim(options(i))=="-deepcat") then
if ((i+1) .le. Narg) then
read (options(i+1),"(A)") deeplinelistfile
options(i:i+1)=""
else
print *,gettime(),"[100] no value specified for ",trim(options(i))
call exit(100)
endif
endif
if ((trim(options(i))=="-g" .or. trim(options(i))=="--generations")) then
if ((i+1) .le. Narg) then
read (options(i+1),*) generations
options(i:i+1)=""
if (generations .lt. 1) then
print *,gettime(),"[100] invalid value given for generations"
call exit(100)
endif
else
print *,gettime(),"[100] no value specified for ",trim(options(i))
call exit(100)
endif
endif
if ((trim(options(i))=="-ps" .or. trim(options(i))=="--populationsize")) then
if ((i+1) .le. Narg) then
read (options(i+1),*) popsize
options(i:i+1)=""
if (popsize .lt. 1) then
print *,gettime(),"[100] invalid value given for popsize"
call exit(100)
endif
else
print *,gettime(),"[100] no value specified for ",trim(options(i))
call exit(100)
endif
endif
if ((trim(options(i))=="-pr" .or. trim(options(i))=="--pressure")) then
if ((i+1) .le. Narg) then
read (options(i+1),*) pressure
options(i:i+1)=""
if (pressure .lt. 0.d0 .or. pressure .gt. 1.d0) then
print *,"[100] pressure must be between 0 and 1"
call exit(100)
endif
else
print *,gettime(),"[100] no value specified for ",trim(options(i))
call exit(100)
endif
endif
if ((trim(options(i))=="-ul" .or. trim(options(i))=="--upper-limits")) then
upperlimits=.true.
options(i)=""
endif
if ((trim(options(i))=="--collapse")) then
collapse=.true.
options(i)=""
endif
if (trim(options(i))=="--citation") then
print *
print *,"ALFA was described in Wesson, 2016, MNRAS, 456, 3774"
print *,"https://ui.adsabs.harvard.edu/abs/2016MNRAS.456.3774W"
print *,"The bibtex data for the paper is:"
print *
print *,"@ARTICLE{2016MNRAS.456.3774W,"
print *," author = {{Wesson}, R.},"
print *,' title = "{ALFA: an automated line fitting algorithm}",'
print *," journal = {\mnras},"
print *," keywords = {line: identification, methods: data analysis, H II regions,"
print *,"planetary nebulae: general, Astrophysics - Solar and Stellar Astrophysics},"
print *," year = 2016,"
print *," month = mar,"
print *," volume = {456},"
print *," number = {4},"
print *," pages = {3774-3781},"
print *," doi = {10.1093/mnras/stv2946},"
print *,"archivePrefix = {arXiv},"
print *," eprint = {1512.04539},"
print *," primaryClass = {astro-ph.SR},"
print *," adsurl = {https://ui.adsabs.harvard.edu/abs/2016MNRAS.456.3774W},"
print *," adsnote = {Provided by the SAO/NASA Astrophysics Data System}"
print *,"}"
call exit(0)
endif
if ((trim(options(i))=="-ws" .or. trim(options(i))=="--wavelength-scaling")) then
if ((i+1) .le. Narg) then
read (options(i+1),*) wavelengthscaling
options(i:i+1)=""
if (wavelengthscaling .lt. 0.d0) then
print *,gettime(),"[100] invalid value given for wavelengthscaling"
call exit(100)
endif
else
print *,gettime(),"[100] no value specified for ",trim(options(i))
call exit(100)
endif
endif
if ((trim(options(i))=="-el" .or. trim(options(i))=="--exclude-line")) then
if ((i+1) .le. Narg) then
read (options(i+1),*) excludewavelength
options(i:i+1)=""
exclusions(exclusioncount) = excludewavelength
exclusioncount = exclusioncount + 1
else
print *,gettime(),"[100] no value specified for ",trim(options(i))
call exit(100)
endif
endif
if ((trim(options(i))=="-dl" .or. trim(options(i))=="--detection-limit")) then
if ((i+1) .le. Narg) then
read (options(i+1),*) detectionlimit
options(i:i+1)=""
if (detectionlimit .lt. 0) then
detectionlimit = 0.d0
print *,gettime(),"warning: negative sigma detection limit specified - has been reset to zero"
endif
else
print *,gettime(),"[100] no value specified for ",trim(options(i))
call exit(100)
endif
endif
if ((trim(options(i))=="-rb" .or. trim(options(i))=="--rebin")) then
if ((i+1) .le. Narg) then
read (options(i+1),*) rebinfactor
options(i:i+1)=""
if (rebinfactor<1) then
print *,gettime(),"[100] impossible rebin factor specified: ",rebinfactor
call exit(100)
endif
else
print *,gettime(),"[100] no value specified for ",trim(options(i))
endif
endif
if ((trim(options(i))=="-nc" .or. trim(options(i))=="--no-continuum")) then
subtractcontinuum=.false.
options(i)=""
endif
if ((trim(options(i))=="-cw" .or. trim(options(i))=="--continuum-window")) then
if ((i+1) .le. Narg) then
read (options(i+1),*) continuumwindow
options(i:i+1)=""
if (continuumwindow .lt. 1) then
print *,gettime(),"[100] invalid value given for continuum window"
call exit(100)
endif
if (mod(continuumwindow,2).eq.1) then
continuumwindow=continuumwindow+1
print *,gettime(),"warning: continuum window has to be an odd number. incremented by one so it's now ",continuumwindow
endif
else
print *,gettime(),"[100] no value specified for ",trim(options(i))
call exit(100)
endif
endif
! to implement:
! continuum percentile
if ((trim(options(i))=="-wc" .or. trim(options(i))=="--wavelength-column")) then
if ((i+1) .le. Narg) then
read (options(i+1),*) tablewavelengthcolumn
options(i:i+1)=""
if (tablewavelengthcolumn .lt. 1) then
print *,gettime(),"[100] invalid value given for table wavelength column"
call exit(100)
endif
else
print *,gettime(),"[100] no value specified for ",trim(options(i))
endif
endif
if ((trim(options(i))=="-fc" .or. trim(options(i))=="--flux-column")) then
if ((i+1) .le. Narg) then
read (options(i+1),*) tablefluxcolumn
options(i:i+1)=""
if (tablefluxcolumn .lt. 1) then
print *,gettime(),"[100] invalid value given for table flux column"
call exit(100)
endif
else
print *,gettime(),"[100] no value specified for ",trim(options(i))
call exit(100)
endif
endif
if ((trim(options(i))=="-of" .or. trim(options(i))=="--output-format")) then
if ((i+1) .le. Narg) then
read (options(i+1),*) outputformat
options(i:i+1)=""
if (outputformat.ne."text".and.outputformat.ne."fits".and.outputformat.ne."latex".and.outputformat.ne."csv") then
print *,gettime(),"[100] invalid output format. valid formats are fits, text, csv and latex"
call exit(100)
endif
else
print *,gettime(),"[100] no value specified for ",trim(options(i))
call exit(100)
endif
endif
if ((trim(options(i))=="-cl" .or. trim(options(i))=="--clobber")) then
clobber=.true.
options(i)=""
endif
enddo
nargused=narg-count(options.ne."")
do i=1,narg
if (len(trim(options(i))).gt.0) then
spectrumfile=options(i)
exit
endif
enddo
!check that an input file was specified and that no unrecognised options are present
if (len(trim(spectrumfile)).eq.0) then
print *,gettime(),"[100] no input file specified"
call exit(100)
elseif (narg - nargused .gt. 1) then
print *,gettime(),"[100] some input options were not recognised:"
do i=1,narg
if (len(trim(options(i))).gt.0) then
print *,trim(options(i))
endif
enddo
call exit(100)
endif
!deal with image sections
if (index(spectrumfile,"[") .gt. 0) then !image section specified
imagesection=spectrumfile(index(spectrumfile,"["):)
spectrumfile=spectrumfile(1:index(spectrumfile,"[")-1)
endif
!check if input file exists
inquire(file=spectrumfile, exist=file_exists) ! see if the input file is present
if (.not. file_exists) then
print *,gettime(),"[101] input spectrum ",trim(spectrumfile)," does not exist"
call exit(101)
endif
deallocate(options)
!display the settings
print *,gettime(),"ALFA is running with the following settings:"
print *," file: ",trim(spectrumfile)
if (len(trim(imagesection)).gt.0) print *," fitting section: ",imagesection
if (.not.normalise) then
print *," normalisation: using measured value of Hb"
else
if (normalisation.eq.0.d0) then
print *," normalisation: no normalisation"
else
print *," normalisation: to Hb=",normalisation
endif
endif
if (subtractcontinuum) then
print *," continuum fitting: enabled"
print *," continuum window: ",continuumwindow
else
print *," continuum fitting: disabled"
endif
print *," spectrum fitted if max value > ",baddata
print *," Angstroms per wavelength unit: ",wavelengthscaling
if (tablewavelengthcolumn.ne.1) then
print *," table wavelength column: ",tablewavelengthcolumn
endif
if (tablefluxcolumn.ne.2) then
print *," table flux column: ",tablefluxcolumn
endif
if (collapse) then
print *," multiple spectra: collapsed to 1D"
else
print *," multiple spectra: fitted individually"
endif
print *," velocity guess: ",redshiftguess_initial
print *," resolution guess: ",resolutionguess_initial
print *," first pass velocity tolerance: ",vtol1*c
print *," second pass velocity tolerance: ",vtol2*c
print *," first pass resolution tolerance: ",rtol1
print *," second pass resolution tolerance: ",rtol2
if (subtractsky) then
print *," sky line fitting: enabled"
print *," sky line catalogue: ",trim(skylinelistfile)
else
print *," sky line fitting: disabled"
endif
print *," strong line catalogue: ",trim(stronglinelistfile)
print *," deep line catalogue: ",trim(deeplinelistfile)
if (exclusioncount .gt. 0) then
print *," lines excluded from fitting: ",exclusions
endif
if (rebinfactor .gt. 1) then
print *," spectra rebinned by factor of: ",rebinfactor
endif
print *," number of generations: ",generations
print *," population size: ",popsize
print *," pressure factor: ",pressure
print *," output directory: ",trim(outputdirectory)
print *," output format: ",outputformat
end subroutine readcommandline
end module mod_commandline
|
