forked from cdominik/optool
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathoptool.f90
3156 lines (2926 loc) · 124 KB
/
optool.f90
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
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
!! **** IO units handling
module IOunits
#ifdef NO_F2003
! ----------------------------------------------------------------------
! IO unit numbers. 5,6,0 should work, but not defined in standard
! ----------------------------------------------------------------------
integer :: stdi = 5
integer :: stdo = 6
integer :: stde = 0
#else
! ----------------------------------------------------------------------
! Use the 2003 standard intrinsic module, modern compilers have it
! ----------------------------------------------------------------------
use, intrinsic :: iso_fortran_env, only : input_unit, output_unit, error_unit
integer :: stdi = input_unit
integer :: stdo = output_unit
integer :: stde = error_unit
#endif
end module IOunits
!!! **** The Usage routine and the module for shared variables
subroutine usage()
use IOunits
write(stdo,'("")')
write(stdo,'("===============================================================================")')
write(stdo,'("optool - dust opacities from the command line")')
write(stdo,'(" Dominik, Min, Tazaki 2021, https://ascl.net/2104.010, version 1.9.13")')
write(stdo,'("")')
write(stdo,'("-c List available materials")')
write(stdo,'("-c KEY-or-FILE [Mfrac] Add material with mass fraction. -c may be omitted")')
write(stdo,'("-m KEY-or-FILE [Mfrac] Add material with mass fraction in mantle")')
write(stdo,'("-p POROSITY [PMANTLE] Set porosity, possibly different for core and mantle")')
write(stdo,'("-dhs [FMAX] Maximum volume fraction of vacuum in DHS computation")')
write(stdo,'("-mmf [A0 [DF-or-FILL]] Use MMF with monomer size A0 and fractal dim or fill")')
write(stdo,'("-a AMIN [AMAX [SD [NA]]] Grain size [micron] dist. SD=P (>PL) or M:S (>log-n)")')
write(stdo,'("-l LMIN [LMAX [NLAM]] Set up wavelength grid (units for a and l: microns)")')
write(stdo,'("-a FILE; -l FILE Read size distribution or lambda grid from a file")')
write(stdo,'("-d [NSUB] Write NA files for specific grain sizes")')
write(stdo,'("-s [NANG] Add scattering matrix for NANG angles to output")')
write(stdo,'("-chop [NDEG] Remove NDEG degrees from the forward-scattering peak")')
write(stdo,'("-o [DIRECTORY] Output to DIRECTORY instead of current working dir")')
write(stdo,'("-radmc [LABEL] Output as RADMC-3D input file")')
write(stdo,'("-fits; -print [Var] Output to FITS file, or to STDOUT")')
write(stdo,'("-h [OPT]; -man Show this msg or help on -OPT; Show the full manual")')
write(stdo,'("-q; -v Quiet or more verbose on STDOUT")')
write(stdo,'("-xlim XLIM Switch to Mie for speed at x=XLIM")')
write(stdo,'("===============================================================================")')
end subroutine usage
module Defs
! Module with constants and the basic data structures used in the program
implicit none
integer, public, parameter :: dp = selected_real_kind(P=15)
! ----------------------------------------------------------------------
! Physics and math constants
! ----------------------------------------------------------------------
real (kind=dp),public,parameter :: pi = 3.1415926535897932384_dp
real (kind=dp),public,parameter :: clight = 2.99792458d10
! ----------------------------------------------------------------------
! Some global switches
! ----------------------------------------------------------------------
logical, public :: blendonly = .false. ! only blend materials and write result out
logical, public :: split = .false. ! split to many files
logical, public :: quiet = .false. ! reduce output to STDOUT
logical, public :: verbose = .false. ! additional output to STDOUT
logical, public :: debug = .false. ! Additional info to STDOUT
logical, public :: write_grd = .false. ! Write out the size distribution and wavelength grid
logical :: mmfss = .false. ! Force single scattering result if phase shift is too large
! ----------------------------------------------------------------------
! Lambda is shared, because multiple routines need it
! ----------------------------------------------------------------------
real (kind=dp), allocatable :: lam(:) ! wavelength array
integer :: nlam ! nr of wavelength points
! ----------------------------------------------------------------------
! Number of angles to be considered
! ----------------------------------------------------------------------
integer :: nang
real (kind=dp) :: chopangle
! ----------------------------------------------------------------------
! Control for sparse scattering files
! ----------------------------------------------------------------------
real (kind=dp) :: scatlammin(30),scatlammax(30) ! 30 should be plenty?
integer :: nsparse=0 ! nr of lam intervals
integer (kind=dp), allocatable :: iscatlam(:) ! flag for scatmat
! ----------------------------------------------------------------------
! Control tricks at large size parameters
! ----------------------------------------------------------------------
real (kind=dp) :: xlim = 1d8 ! size parameter to switch to Mie for speed
real (kind=dp) :: xlim_dhs = 1d4 ! size parameter limit in DHS for stability
! ----------------------------------------------------------------------
! Material properties
! ----------------------------------------------------------------------
integer :: mat_nm ! number of materials specified
integer :: mat_nmc ! number of core materials specified
integer :: mat_nmm ! number of mantle materials specified
character*500 :: mat_loc(21) ! either 'core' or 'mantle'
character*500 :: mat_lnk(21) ! the lnk key of file path
logical :: mat_cmd(21) ! are n, k, and rho from the command line?
real (kind=dp) :: mat_rfn(21) ! real refractive index if specified on command line
real (kind=dp) :: mat_rfk(21) ! imag refractive index if specified on command line
real (kind=dp) :: mat_rho(21) ! specific mass density of material
real (kind=dp) :: mat_mfr(21) ! mass fraction of each component
real (kind=dp), allocatable :: mat_e1(:,:) ! Real part of refractive index
real (kind=dp), allocatable :: mat_e2(:,:) ! Imaginary part of refractive index
! ----------------------------------------------------------------------
! Mueller matrix structure, records only non-zero elements of the matrix
! ----------------------------------------------------------------------
type mueller
real (kind=dp),allocatable :: F11(:),F12(:),F22(:),F33(:),F44(:),F34(:)
end type mueller
! ----------------------------------------------------------------------
! The particle structure, contains particle and scattering properties
! ----------------------------------------------------------------------
type particle
real (kind=dp) :: rv,rvmin,rvmax ! grain radius with min and max
real (kind=dp) :: rho ! mass density in g.cm^-3
real (kind=dp), allocatable :: Kabs(:),Ksca(:),Kext(:) ! Opacities
real (kind=dp), allocatable :: g(:) ! asymmetry parameter
TYPE(MUELLER), allocatable :: F(:) ! Mueller matrix elements
logical , allocatable :: testscat(:) ! Can we trust the scattering matrix?
logical :: scat_ok ! Are F11... and g_asym usable?
real (kind=dp) :: scat_ok_lmin ! last lambda with bad scattering
end type particle
! ----------------------------------------------------------------------
! The output directory and other strings
! ----------------------------------------------------------------------
character*500 :: outdir = '' ! Output directory
character*500 :: sdfile = '' ! Size distribution file
character*500 :: sdoutfile = '' ! Size distribution file
character*500 :: lamoutfile = '' ! Size distribution file
character*500 :: blendoutfile = '' ! Blended optical constants file
real (kind=dp) :: ameans_file(3) ! for size means from file
character*3 :: method ! DHS or MMF
character*4 :: sdkind ! apow, lgnm, norm, or file
character*1 :: justnum = ' ' ! What to print to STDOUT
end module Defs
!!! **** Main program and ComputePart
program optool
use Defs
use IOunits
use omp_lib
implicit none
integer :: na ! nr of sizes for size distribution
real (kind=dp) :: amin,amax ! min and max size of grains
real (kind=dp) :: apow ! power law index f(a) ~ a^(-apow)
real (kind=dp) :: amean,asig ! mean and standard deviation for lognormal f(a)
real (kind=dp) :: fmax ! maximum fraction of vaccum for DHS
real (kind=dp) :: pcore, pmantle ! porosity for core and mantle
real (kind=dp) :: lmin,lmax,l1,l2 ! min and max wavelength
logical :: write_scatter ! Should the scattering matrix be written?
logical :: write_fits ! Should a fits file be written?
logical :: for_radmc ! Should the scattering matrix use RADME-3D convention?
integer :: nm ! nr of grain materials
integer :: nmant ! nr of mantle materials
integer :: it,il
type(particle) :: p
integer :: i,ndone ! counter
integer :: im,ia ! for material, radius
character*1000 :: tmp,value,sub ! for processing args
character*100 :: feature ! for processing args
! String scanning functions, mostly used in argument checking
logical :: string_starts_like_number
logical :: string_is_number
logical :: string_is_n_numbers
integer :: count_colons
! Functions to test command line arguments
logical :: arg_is_present ! functions to test arguments
logical :: arg_is_switch ! functions to test arguments
logical :: arg_is_value ! functions to test arguments
logical :: arg_is_number ! functions to test arguments
logical :: arg_is_1_number ! functions to test arguments
logical :: arg_is_n_numbers! functions to test arguments
logical :: is_key_or_file ! functions to test arguments
logical :: is_file ! functions to test arguments
logical :: file_exists ! return value
! Store and check for file names
character*500 :: fitsfile ! file name for FITS output
character*500 :: meanfile ! file name for mean opacity output
character*500 :: make_file_path ! function
character*50 :: radmclbl = "" ! file label for RADMC-3D compatible
character*50 :: label ! for use in file names
! Temporary storage vaiables
character*500 :: dumc ! temporary storage
real (kind=dp) :: dum ! temporary storage
logical :: duml ! temporary storage
real (kind=dp), allocatable :: e1d(:),e2d(:)
! Variables for splitting size range
real (kind=dp) :: asplit,afact,afsub,amaxsplit,aminsplit
integer :: nsubgrains = 5,nsub
! MMF implementation
real(kind=dp) :: mmf_a0,mmf_struct,mmf_kf
! Variables for fine-tuning the -pring mechaanism
integer iop,icp
real (kind=dp) justnum_angle
! ----------------------------------------------------------------------
! Defaults values for parameters and switches
! ----------------------------------------------------------------------
amin = 0.05_dp ! micrometer
amax = 3000._dp ! micrometer
apow = 3.50_dp ! a minus sign will be added internally
amean = 0. ! a0 in micrometer for log-normal distribution
asig = 0. ! Sigma, not units, for log-normal distribution
na = 0 ! will be computed to 10 per decade
sdkind = 'apow' ! size distribution method
lmin = 0.05_dp ! micrometer
lmax = 10000.0_dp ! micrometer
nlam = 300 ! size of wavelengths grid
nang = 180 ! Number of angular point, has to be even
chopangle = 0.d0 ! Angle in degree to chop forward peak
pcore = 0.0_dp ! porosity core
pmantle = 0.0_dp ! porosity mantle
nm = 0 ! number of materials - zero to start with
nmant = 0 ! number of mantle materials - zero to start with
method = 'DHS' ! method to compute the opacities
fmax = 0.8_dp ! maximum volume fraction DHS
write_fits = .false. ! Default is to write ASCII output
write_scatter = .false. ! Default is to not write scattering matrix
for_radmc = .false. ! Default is to use optool conventions.
! ----------------------------------------------------------------------
! Initialize rho, because we need the fact that it has not been set
! to decide what to do with lnk files where it is missing
! ----------------------------------------------------------------------
mat_rho(:) = 0.d0
! ----------------------------------------------------------------------
! Process the command line arguments
! ----------------------------------------------------------------------
! Prescan to get -q,-v,-debug, so they can be active during arg processing
i = 1; call getarg(i,tmp)
do while(tmp.ne.' ')
if (tmp.eq.'-q') quiet = .true.
if (tmp.eq.'-v') verbose = .true.
if (tmp.eq.'-debug') debug = .true.
i = i+1; call getarg(i,tmp)
enddo
! Loop over all command line arguments
i = 1; call getarg(i,tmp)
do while(tmp.ne.' ')
! If the are two dashes, keep only one.
if (tmp(1:2).eq.'--') tmp = tmp(2:)
select case(tmp)
case('?','-h','-help','help')
if (.not. arg_is_present(i+1)) then
call usage(); stop
endif
i=i+1
call getarg(i,value)
if (value(1:2).eq.'--') value = value(2:)
if (.not. (value(1:1).eq.'-')) value = '-' // value
call manual(trim(value))
stop
case('??','-man','-manual')
call manual('all'); stop
case('-version')
write(stdo,*) "OpTool version 1.9.13, June 2023, (c) C. Dominik, M. Min & R. Tazaki"
stop
! ----------------------------------------------------------------------
! Definition of a material for the mix
! ----------------------------------------------------------------------
case('-c','-m')
i = i+1
nm = nm+1;
if (nm .gt. 20) then
write(stde,*) 'ERROR: too many materials'; stop
endif
! First value is the material key or refindex file path
call getarg(i,value); read(value,'(A)') mat_lnk(nm)
if ((value .eq. '?') .or. (value .eq. '')) then
call ListBuiltinMaterials(); stop
endif
! Set the type
if (tmp.eq.'-m') then
! This is the mantle material
mat_loc(nm) = 'mantle'; nmant = nmant+1
else
! This is a core material.
mat_loc(nm) = 'core'
endif
if (arg_is_n_numbers(i,3)) then
! This is a material with specified n and k and rho
it = index(value,':')
if (.not. quiet) write(stde,*) &
"WARNING: 3-value list read as material with n:k:rho: ",trim(value)
read(value(1:it-1),*) mat_rfn(nm)
value=value(it+1:)
it = index(value,':')
if (it.gt.0) then
read(value(1:it-1),*) mat_rfk(nm)
read(value(it+1:len(value)),*) mat_rho(nm)
if (mat_rho(nm) .le. 0.d0) then
write(stde,*) 'ERROR: Density must be larger than zero: ',mat_rho(nm);
stop
endif
else
write(stde,*) 'ERROR: Please specify a density for artificial material: -c n:k:rho';
stop
endif
mat_cmd(nm) = .true.
else
mat_cmd(nm) = .false.
! Check if this is a valid material key or a file
if (.not. is_key_or_file(trim(value),.true.)) then
write(stde,*) "ERROR: not a material key or lnk file: ",trim(value)
stop
endif
endif
! Second value is the mass fraction
if (.not. arg_is_1_number(i+1)) then
if (.not. quiet) write(stde,*) "WARNING: 1.0 used for missing mass fraction of material: ",trim(mat_lnk(nm))
mat_mfr(nm) = 1.0d0
goto 997
else
i = i+1; call getarg(i,value); read(value,*) mat_mfr(nm)
endif
if (mat_mfr(nm).eq.0d0) then
if (.not. quiet) write(stde,*) &
"WARNING: Ignoring material with zero mass fraction: ",trim(mat_lnk(nm))
nm = nm-1
endif
! There might be a density, as a third argument
if (arg_is_1_number(i+1)) then
i = i+1; call getarg(i,value); read(value,*) mat_rho(nm)
if (mat_cmd(nm)) then
write(stde,*) "ERROR: rho specified twice (in name and after MFRAC) for material: ", &
trim(mat_lnk(nm))
stop
endif
if (mat_rho(nm) .le. 0.d0) then
write(stde,*) 'ERROR: Density must be larger than zero: ',mat_rho(nm);
stop
endif
endif
997 continue
! ----------------------------------------------------------------------
! Special compositions known in the literature
! ----------------------------------------------------------------------
case('-diana','-dsharp','-dsharp-no-ice')
if (nm.gt.0) then
write(stde,*) "ERROR: Standard mixtures must be specified before any additional materials"
stop
endif
if (tmp.eq.'-diana') then
nm = 2
mat_lnk(1) = 'pyr-mg70'; mat_loc(1) = 'core'; mat_mfr(1) = 0.87d0
mat_lnk(2) = 'c-z' ; mat_loc(2) = 'core'; mat_mfr(2) = 0.1301d0
pcore =0.25d0
elseif (tmp.eq.'-dsharp') then
nm = 4
mat_lnk(1) = 'astrosil'; mat_loc(1) = 'core'; mat_mfr(1) = 0.3291d0
mat_lnk(2) = 'c-org' ; mat_loc(2) = 'core'; mat_mfr(2) = 0.3966d0
mat_lnk(3) = 'fes' ; mat_loc(3) = 'core'; mat_mfr(3) = 0.0743d0
mat_lnk(4) = 'h2o-w' ; mat_loc(4) = 'core'; mat_mfr(4) = 0.2000d0
pcore = 0.d0
elseif (tmp.eq.'-dsharp-no-ice') then
nm = 3
mat_lnk(1) = 'astrosil'; mat_loc(1) = 'core'; mat_mfr(1) = 0.3291d0
mat_lnk(2) = 'c-org' ; mat_loc(2) = 'core'; mat_mfr(2) = 0.3966d0
mat_lnk(3) = 'fes' ; mat_loc(3) = 'core'; mat_mfr(3) = 0.0743d0
pcore = 0.d0
endif
! ----------------------------------------------------------------------
! Grain size setup
! ----------------------------------------------------------------------
case('-a')
! ----------------------------------------------------------------------
! -a expects 1, 2, 3, or 4 values: amin [amax [apow [na]]]
! amin amax amean:asig [na]
! sizedist_file
! ----------------------------------------------------------------------
if (.not. arg_is_1_number(i+1)) then
if (arg_is_present(i+1)) then
if (arg_is_switch(i+1)) then
write(stde,*) "ERROR: -a needs 1-4 values: amin [amax [na [apow]]]"; stop
endif
i=i+1
call getarg(i,sdfile)
inquire (file=trim(sdfile),exist=file_exists)
if (file_exists) then
call checkout_sdfile(sdfile,amin,amax,na)
sdkind = 'file'
else
write(stde,*) "Size distribution file does not exist: ",trim(sdfile)
stop
endif
else
write(stde,*) "ERROR: -a needs 1-4 values: amin [amax [na [apow]]]"; stop
endif
endif
if (trim(sdfile).eq.'') then
! There are numbers present, get the values
i=i+1; call getarg(i,value); call uread(value,amin)
! Let's see if there is more, we expect amax
if (arg_is_1_number(i+1)) then
i=i+1; call getarg(i,value); call uread(value,amax)
if (amax .lt. 0d0) then
! FIXME: Take this out?
if (amin+amax .le. 0d0) then
write(stde,'(" ERROR: delta a cannot be larger than a: ",F10.2,F10.2)') amin,amax
stop
endif
amin = amin+amax; amax = amin-2d0*amax
apow = 0.d0
endif
! Let's see if there is more, we expect apow, or other sizedistribution parameters.
! First we exclude switches and material specifications
if (arg_is_present(i+1) .and. (.not. arg_is_switch(i+1)) &
.and. arg_is_number(i+1) .and. (.not. arg_is_n_numbers(i+1,3))) then
! OK, we have something describing a size distribution
i=i+1; call getarg(i,value);
it = index(value,':')
if (it .gt. 0) then
! We seem to have 2 numbers; so we have a (log-)normal size distribution
read(value(1:it-1),*) amean
! The following check should no longer be necessary.
if (index(value(it+1:len(value)),':').gt.0) then
write(stde,'(" ERROR: Problems interpreting argument: ",A)') trim(value)
stop
endif
read(value(it+1:len(value)),*) asig
sdkind = 'norm' ! could still also be lgnm, decide later
else if (arg_is_1_number(i)) then
! Just one number - we have a powerlaw
read(value,*) apow
sdkind = 'apow'
endif
! Let's see if there is more, we expect na
if (arg_is_1_number(i+1)) then
i=i+1; call getarg(i,value); read(value,*) na
endif
endif
else
! There was only 1 number after -a. Set up single grain size computation
! If na has already been set, do not change it - we will check this later
amax = amin;
if (na.eq.0) na = 1;
endif
endif
case('-amin')
i = i+1; call getarg(i,value); call uread(value,amin)
case('-amax')
i = i+1; call getarg(i,value); call uread(value,amax)
case('-apow')
i = i+1; call getarg(i,value); read(value,*) apow
case('-amean')
i = i+1; call getarg(i,value); read(value,*) amean
sdkind = 'norm'
case('-asig')
i = i+1; call getarg(i,value); read(value,*) asig
sdkind = 'norm'
case('-na')
i = i+1; call getarg(i,value); read(value,*) na
! ----------------------------------------------------------------------
! Wavelength setup
! ----------------------------------------------------------------------
case('-l')
! ----------------------------------------------------------------------
! -l expects a file name, or 1-3 numbers: lmin [lmax [nlam]]
! ----------------------------------------------------------------------
if (.not. arg_is_value(i+1)) then
write(stde,*) "ERROR: -l needs a file or numbers as values"; stop
else if (.not. arg_is_number(i+1)) then
! Could be a file name. If yes, read the lambda grid from it
call getarg(i+1,value)
call require_file(trim(value))
i=i+1
call read_lambda_grid(trim(value))
lmin = minval(lam); lmax = maxval(lam)
else
! We have a number, should be lmin
i = i+1; call getarg(i,value); call uread(value,lmin)
! Let's see if there is more, we expect lmax
if (arg_is_1_number(i+1)) then
i = i+1; call getarg(i,value); call uread(value,lmax)
! Let's see if there is more, we expect nlam
if (arg_is_1_number(i+1)) then
i=i+1; call getarg(i,value); if (value.ne.'') read(value,*) nlam
endif
else
! only lmin was given, set up single lambda computation
lmax = lmin; nlam = 1;
endif
endif
case('-lmin')
i = i+1; call getarg(i,value); call uread(value,lmin)
case('-lmax')
i = i+1; call getarg(i,value); call uread(value,lmax)
case('-nlam','-nl')
i = i+1; call getarg(i,value); read(value,*) nlam
! ----------------------------------------------------------------------
! Grain geometry, including method DHS versus MMF
! ----------------------------------------------------------------------
case('-p','-porosity')
i = i+1; call getarg(i,value); read(value,*) pcore
if (arg_is_1_number(i+1)) then
i=i+1; call getarg(i,value); read(value,*) pmantle
else
pmantle = pcore
endif
case('-dhs','-fmax','-mie')
method = 'DHS'
if (tmp.eq.'-mie') then
fmax = 0.
else if (arg_is_1_number(i+1)) then
i = i+1; call getarg(i,value); read(value,*) fmax
else
fmax = 0.8
endif
case('-xlim')
if (arg_is_1_number(i+1)) then
i = i+1; call getarg(i,value); read(value,*) xlim
else
write(stde,*) "ERROR: -xlim needs a numeric value"; stop
endif
case('-xlim_dhs')
if (arg_is_1_number(i+1)) then
i = i+1; call getarg(i,value); read(value,*) xlim_dhs
else
write(stde,*) "ERROR: -xlim_dhs needs a numeric value"; stop
endif
case('-mmf','-mmfss')
method = 'MMF'
if (tmp.eq.'-mmfss') then
! FIXME: Should this be wrapped into a test for -q ?
write(stde,*) "WARNING: We will use the assumption of single scattering to compute the"
write(stde,*) " MMF matrix elements when the phase shift is too large. See UserGuide."
mmfss = .true.
endif
! Defaults
mmf_a0 = 0.1d0
mmf_struct = 0.2d0
mmf_kf = 0.d0
if (arg_is_1_number(i+1)) then
i=i+1; call getarg(i,value); call uread(value,mmf_a0)
if (arg_is_1_number(i+1)) then
i=i+1; call getarg(i,value); read(value,*) mmf_struct
if (arg_is_1_number(i+1)) then
i=i+1; call getarg(i,value); read(value,*) mmf_kf
endif
endif
endif
case('-cde')
method = 'CDE'
! ----------------------------------------------------------------------
! Various other switches
! ----------------------------------------------------------------------
case('-o')
if (arg_is_value(i+1)) then
i=i+1; call getarg(i,outdir)
else
outdir = 'output'
endif
case('-s','-scatter','-scat')
write_scatter=.true.
if (arg_is_1_number(i+1)) then
! Change size of the angle grid
i=i+1; call getarg(i,value); read(value,*) nang
endif
case('-sp','-sparse')
write_scatter=.true.
if (arg_is_1_number(i+1) .and. (arg_is_1_number(i+2))) then
! Two numbers. This is a wavelength range for a sparse file
nsparse = nsparse+1
if (nsparse.gt.10) then
write(stde,*) "ERROR: To many sparse file ranges (10 is max)"; stop
endif
i=i+1; call getarg(i,value); call uread(value,l1)
i=i+1; call getarg(i,value); call uread(value,l2)
if (l1.lt.l2) then
scatlammin(nsparse) = l1; scatlammax(nsparse) = l2
else
scatlammin(nsparse) = l2; scatlammax(nsparse) = l1
endif
else if (arg_is_1_number(i+1)) then
! One number. One specific wavelength
i=i+1; call getarg(i,value); call uread(value,l1)
nsparse = nsparse+1
scatlammin(nsparse) = l1; scatlammax(nsparse) = l1
else
write(stde,*) "ERROR: -sparse needs one or two wavelengths"
stop
endif
case('-chop')
if (.not. arg_is_value(i+1)) then
chopangle = 2.d0
else
i=i+1; call getarg(i,value); read(value,*) chopangle
endif
case('-radmc','-radmc3d')
for_radmc = .true.
if (arg_is_value(i+1)) then
call getarg(i+1,value)
if ((.not. quiet) .and. &
(is_key_or_file(trim(value),.false.) .or. &
string_is_n_numbers(trim(value),3))) then
! The optional -radmc label is also a material key, this is ambiguous
write(stde,*) "WARNING: Ambiguous argument could be meant as (another) material key"
write(stde,*) " ... but is read as optional RADMC-3D label: ",trim(tmp)," ",trim(value)
write(stde,*) " ... Use -c or reorder args to disambiguate"
endif
if (is_file(trim(value)) .or. (scan(value,"/").gt.0)) then
! It is a file, treat as core material
tmp = 'FoRcE_-c_FoRcE'
i=i+1
else
i=i+1
call getarg(i,radmclbl)
endif
endif
case('-fits')
write_fits = .true.
case('-d')
split = .true.
if (arg_is_value(i+1)) then
i=i+1; call getarg(i,value); read(value,*) nsubgrains
endif
case('-q')
! Be less noisy
quiet = .true.
case('-v')
! Be more noisy
verbose = .true.
case('-print')
! FIXME: quiet = .true.
justnum = 'x'
justnum_angle = -1.
if (arg_is_value(i+1)) then
i=i+1; call getarg(i,value)
iop = index(value,'[')
icp = index(value,']')
if ((iop.gt.0).and.(icp.gt.0).and.(iop.lt.icp-1)) then
read(value(iop+1:icp-1),*) justnum_angle
value = trim(value(1:iop-1))
endif
select case(trim(value))
case('?')
! Show what option are available for -print
write(stdo,*) 'Keywords for -print KEYWORD. Output to STDOUT, for each lambda.'
write(stdo,*) 'Headers will be printed to STDERR for viewing convenience - '
write(stdo,*) ' they will not show up in a pipe or redirect.'
write(stdo,*) ' KEYS | ACTION'
write(stdo,*) '-----------------------------------------------------------------------------'
write(stdo,*) ' : print LAM KABS KSCA KEXT G'
write(stdo,*) 'a | kabs : print the absorption opacity'
write(stdo,*) 's | ksca : print the scattering opacity'
write(stdo,*) 'e | kext : print the extinction opacity'
write(stdo,*) 'g | gsca : print the asymmetry factor G'
write(stdo,*) 'm | lnk : print mixed refractive index LAM N K'
write(stdo,*) 'f | fmat : print ANG F11 F12 F22 F33 F34 F44. f[67] selects just ~67 deg.'
write(stdo,*) 'p | pol : print ANG F11 F12 |F12/F11| . p[21] selects just ~21 deg.'
stop
! Select the right single-character key
case('all') ; justnum = 'x'
case('kabs','a') ; justnum = 'a'
case('ksca','kscat','s') ; justnum = 's'
case('kext','e') ; justnum = 'e'
case('g','gsca','gscat') ; justnum = 'g'
case('lnk','m') ; justnum = 'l'
case('f','F') ; justnum = 'f'
case('ip','p') ; justnum = 'p'
case default
i=i-1
write(stde,*) 'WARNING: "',trim(value),'" is not a -print key. Trying core material...';
endselect
endif
case ('-tex')
! run optool2tex with the same command line arguments
call run_optool2tex()
case('-debug')
! More info to STDOUT
debug = .true.
case('-w')
! Write some files and exit without doing a computation
write_grd = .true.
blendonly = .true.
case('-wgrid')
! Write the sitze distribution sizedist.dat
write_grd = .true.
case('-b','-blendonly')
! Write blended refractive index to file and exit
blendonly = .true.
case('-feature')
if (.not. arg_is_value(i+1)) then
write(stde,*) "ERROR: -feature switch needs value"
stop
endif
i=i+1; call getarg(i,value); read(value,*) feature
select case(trim(feature))
case('multi')
#ifdef USE_MULTI
write(stdo,'(A)') 'True'
#else
write(stdo,'(A)') 'False'
#endif
stop
case('fits')
#ifdef USE_FITSIO
write(stdo,'(A)') 'True'
#else
write(stdo,'(A)') 'False'
#endif
stop
case('gfortran')
#ifdef USE_GFORTRAN
write(stdo,'(A)') 'True'
#else
write(stdo,'(A)') 'False'
#endif
stop
case('ifort')
#ifdef USE_IFORT
write(stdo,'(A)') 'True'
#else
write(stdo,'(A)') 'False'
#endif
stop
case default
write(stde,*) "ERROR: No information for feature: ",trim(feature)
stop
end select
case default
if (arg_is_switch(i)) then
write(stde,*) "ERROR: Option or Arg: >",trim(tmp),'> not recognized'
write(stde,*) "For help, try: optool -h ... or find the user guide OpTool.pdf"
stop
else
if (debug) write(stde,*) trim(tmp),' will be interpreted as a material'
tmp = 'FoRcE_-c_FoRcE'
endif
end select
if (tmp .eq. 'FoRcE_-c_FoRcE') then
tmp = '-c'
i = i-1
else
i = i+1
call getarg(i,tmp) ! Get the next argument for the loop.
endif
enddo
! ----------------------------------------------------------------------
! Sanity checks and preparations
! ----------------------------------------------------------------------
! *** Materials
if (nm .ge. 20) then
write(stde,*) 'ERROR: Too many materials'; stop
endif
if ( (nm.eq.nmant) .and. (nm.gt.0) ) then
write(stde,*) "ERROR: at least one core material must be specified"; stop
endif
! *** Porosity
if ( (pcore.lt.0d0).or.(pcore.ge.1d0).or.(pmantle.lt.0d0).or.(pmantle.ge.1d0) ) then
write(stde,*) "ERROR: prosities must be 0 <= p < 1"; stop
endif
if ((justnum .ne. ' ') .and. split) then
write(stde,*) "ERROR: -print does not work in connection with -d splitting"; stop
endif
! *** Grain size distribution
if ( (amin.le.0d0) .or. (amax.le.0d0) ) then
write(stde,*) 'ERROR: Both amin and amax need to be positive numbers',amin,amax; stop
endif
if (amin .gt. amax) then
! Swap min and max values
dum = amin; amin = amax; amax = dum
endif
if (na .eq. 0) then
! set sampling of the grain radius: 15 per decade, min 5
na = max(5,int((log10(amax)-log10(amin))*15d0+1d0))
endif
if ( (amin.eq.amax) .and. (na.ne.1) ) then
write(stde,*) 'WARNING: Setting na=1 because amin=amax'
na = 1
endif
if (sdkind .eq. 'apow') then
if (apow .lt. 0d0) then
write(stde,*) 'WARNING: Unusual negative value for apow. apow=-3 means f(a)~a^(+3)'
endif
amean = 0.d0; asig=0.d0
else if (sdkind .eq. 'file') then
!amean = 0.d0; asig=0.d0; apow = 0.d0
else if (sdkind .eq. 'norm') then
!apow = 0.d0
if (amean .le. 0.d0) then
write(stde,*) "ERROR: amean must be positive for (log-)normal distribution"
stop
endif
if (asig .eq. 0.d0) then
write(stde,*) "ERROR: asig cannot be zero for (log-)normal distribution"
stop
else if (asig .gt. 0.d0) then
sdkind = 'lgnm'
else if (asig .lt. 0.d0) then
sdkind = 'norm' ! redundant, but for clarity
endif
endif
! *** Wavelength grid
if ( (lmin.le.0d0) .or. (lmax.le.0d0) ) then
write(stde,*) 'ERROR: Both lmin and lmax need to be positive numbers',lmin,lmax; stop
endif
if (lmin .gt. lmax) then
! Swap min and max values
dum = lmin; lmin = lmax; lmax = dum
endif
if ( (nlam.le.1) .and. (lmin.ne.lmax)) then
write(stde,*) 'ERROR: More than one wavelength point needed to sample a range',lmin,lmax,nlam; stop
endif
if ( (lmin.eq.lmax) .and. (nlam.ne.1) ) then
write(stde,*) 'WARNING: Setting nlam=1 because lmin=lmax'
nlam = 1
endif
if ((nsparse.gt.0) .and. (.not. quiet)) then
write(stde,*) 'WARNING: Creating a sparse scattering matrix file'
endif
! *** DHS
if (method .eq. 'DHS') then
if ((fmax.lt.0.d0) .or. (fmax.ge.1.d0)) then
write(stde,*) 'ERROR: fmax for DHS must be >0 and <1'; stop
endif
endif
! *** MMF
if (method .eq. 'MMF') then
if (mmf_struct .gt. 3.0d0) then
write(stde,*) 'ERROR: Fractal dimension needs to be between 1 and 3'; stop
endif
if (mmf_struct .le. 0d0) then
write(stde,*) 'ERROR: MMF structure parameter needs to be positive'; stop
endif
if (mmf_a0 .ge. amin) then
write(stde,*) 'ERROR: Minimum grain size cannot be smaller than monomer size'; stop
endif
endif
! *** CDE
if (method .eq. 'CDE') then
if (lmin .le. 2.d0*pi*amax) then
write(stde,'("WARNING: CDE requires Rayleigh limit, but 2 pi a_max/lambda_min =",1p,e8.1)') &
2.d0*pi*amax/lmin
endif
endif
! *** Angular grid
if (mod(nang,2) .eq. 1) then
write(stde,*) 'ERROR: The number of angles in -s NANG must be even'
stop
endif
! *** Other checks
if (split .and. blendonly) then
if (.not. quiet) write(stde,*) 'WARNING: Turning off -s for -blendonly'
split = .false.
endif
if (split .and. (sdkind .ne. 'apow')) then
write(stde,*) "ERROR: Please only use -d with a powerlaw size distribution"
stop
endif
! *** Output files
#ifndef USE_FITSIO
if (write_fits) then
write(stde,*) 'ERROR: Support for writing FITS files needs to be compiled in.'
write(stde,*) ' If you want FITS output, make sure cfitsio library is installed."'
write(stde,*) ' Then recompile with: "make clean", and then "make fitsio=true"'
stop
endif
#endif
if (trim(outdir) .ne. '') then
call make_directory(outdir)
endif
meanfile = make_file_path(outdir,"dustkapmean.dat")
fitsfile = make_file_path(outdir,"dustkappa.fits")
sdoutfile = make_file_path(outdir,"optool_sd.dat")
lamoutfile = make_file_path(outdir,"optool_lam.dat")
blendoutfile = make_file_path(outdir,"optool_mix.lnk")
! ----------------------------------------------------------------------
! Default grain composition if nothing is specified
! ----------------------------------------------------------------------
if (nm.eq.0) then
! Set default composition will set a default composition here, the DIANA opacities
if (.not. quiet) then
write(stde,'("No materials specified, using DIANA standard")')
endif
nm = 2
mat_lnk(1) = 'pyr-mg70' ; mat_loc(1) = 'core' ; mat_mfr(1) = 0.87d0
mat_lnk(2) = 'c-z' ; mat_loc(2) = 'core' ; mat_mfr(2) = 0.1301d0
pcore = 0.25d0
endif
mat_nm = nm
! ----------------------------------------------------------------------
! Sort the materials to make sure that core materials come first
! ----------------------------------------------------------------------
it = nm ! target: where to move the next mantle material
mat_nmc = nm; mat_nmm = 0
do il = nm,1,-1 ! il loops down, checking all materials
if (mat_loc(il).eq.'mantle') then
mat_nmc = mat_nmc-1
mat_nmm = mat_nmm+1
if (il.lt.it) then
dumc = mat_lnk(il); mat_lnk(il)=mat_lnk(it); mat_lnk(it)=dumc;
dumc = mat_loc(il); mat_loc(il)=mat_loc(it); mat_loc(it)=dumc;
dum = mat_mfr(il); mat_mfr(il)=mat_mfr(it); mat_mfr(it)=dum;
dum = mat_rho(il); mat_rho(il)=mat_rho(it); mat_rho(it)=dum;
duml = mat_cmd(il); mat_cmd(il)=mat_cmd(it); mat_cmd(it)=duml;
dum = mat_rfn(il); mat_rfn(il)=mat_rfn(it); mat_rfn(it)=dum;
dum = mat_rfk(il); mat_rfk(il)=mat_rfk(it); mat_rfk(it)=dum;
it = it-1
endif
endif
enddo
! ----------------------------------------------------------------------
! Make a logarithmic lambda grid, unless read in from file
! ----------------------------------------------------------------------
if (allocated(lam)) then
! Lam was allocated by reading from a file
continue
else
allocate(lam(nlam))
if (nlam.eq.1) then
lam(1) = lmin
else
do i = 1,nlam
lam(i)=10.0_dp**(log10(lmin)+log10(lmax/lmin)*(i-1)/(nlam-1))
enddo
endif
endif
allocate(iscatlam(nlam)); iscatlam=0
if (nsparse.gt.0) then
call prepare_sparse()
endif
if (write_grd) then
if (.not. quiet) write(stde,'("Writing wavelength grid to file optool_lam.dat")')
open(unit=20,file=lamoutfile)
write(20,'("# Wavelength grid written by optool, can be read back in with -l optool_lam.dat")')
write(20,'("# First line: number of wavelengths")')
write(20,'("# Then one lambda per line, in micrometer")')
write(20,*) nlam
do i=1,nlam
write(20,'(e18.5)') lam(i)
enddo
close(unit=20)
endif
! Allocate space for the refractive indices
allocate(mat_e1(nm+1,nlam),mat_e2(nm+1,nlam))
! ----------------------------------------------------------------------
! Get the refractory index data for all materials
! ----------------------------------------------------------------------
allocate(e1d(nlam),e2d(nlam))
do im=1,nm
if (mat_cmd(im)) then
! This is a case where n and k were given on the command line