forked from cdominik/optool
-
Notifications
You must be signed in to change notification settings - Fork 0
/
optool_manual.f90
1046 lines (1046 loc) · 66.4 KB
/
optool_manual.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
subroutine manual(what)
character*(*) what
write(*,*)
if (what.eq.'all') then
write(*,'("1 Introduction")')
write(*,'("==============")')
write(*,'("")')
write(*,'(" This tool produces complex dust particle opacities right from the")')
write(*,'(" command line. It is derived from Michiel Min s DHS [OpacityTool] and")')
write(*,'(" also implements Ryo Tazaki s MMF theory for highly porous aggregates.")')
write(*,'("")')
write(*,'("1.1 Capabilities")')
write(*,'("~~~~~~~~~~~~~~~~")')
write(*,'("")')
write(*,'(" - stand-alone tool, fully command line driven, no input files need to")')
write(*,'(" be edited")')
write(*,'(" - full scattering matrix output in several formats, including for")')
write(*,'(" RADMC-3D")')
write(*,'(" - combining materials through mixing into a complex grain with")')
write(*,'(" porosity")')
write(*,'(" - /built-in/: a curated collection of materials for applications in")')
write(*,'(" astronomy")')
write(*,'(" - external refractive index data can be used just as easily")')
write(*,'(" - computational methods: (i) *DHS (Distribution of Hollow Spheres)*")')
write(*,'(" for /irregular grains/ and /low-porosity/ aggregates. Standard *Mie")')
write(*,'(" theory* for /perfect spheres/ is available as a limiting case. (ii)")')
write(*,'(" *MMF (Modified Mean Field)* theory for /high-porosity/fractal")')
write(*,'(" aggregates/. (iii) *CDE* approximation in the Rayleigh limit.")')
write(*,'(" - Python interface module for plotting and post-processing")')
write(*,'("")')
write(*,'("1.2 Terms of use")')
write(*,'("~~~~~~~~~~~~~~~~")')
write(*,'("")')
write(*,'(" optool is distributed under the [MIT license] and can be used,")')
write(*,'(" changed and redistributed freely. But we do ask you to provide a")')
write(*,'(" reference to optool when using it. Relevant references are listed")')
write(*,'(" below and the corresponding BibTeX entries are available in the file")')
write(*,'(" optool.bib . optool is [hosted on github].")')
write(*,'("")')
write(*,'(" - *optool:* [Dominik, C., Min, M. & Tazaki, R. 2021, Optool, 1.9,")')
write(*,'(" Astrophysics Source Code Library, ascl:2104.010]")')
write(*,'(" - *DHS model for irregular grains:* [Min, M. et al. 2005, A&A, 432,")')
write(*,'(" 909]")')
write(*,'(" - *MMF model for aggregates:* [Tazaki, R. & Tanaka, H. 2018, ApJ 860,")')
write(*,'(" 79]")')
write(*,'(" - *DIANA standard opacities:* [Woitke, P. et al. 2016, A&A 586, 103]")')
write(*,'(" - References to [refractive index data] used in your particular")')
write(*,'(" application.")')
write(*,'("")')
write(*,'("1.3 Physical units used by optool")')
write(*,'("~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~")')
write(*,'("")')
write(*,'(" Due to conventions in our field, the input and output of optool uses")')
write(*,'(" the following units:")')
write(*,'(" grain sizes and wavelengths.[1] mum ")')
write(*,'(" mass densities of materials g cm^-3 ")')
write(*,'(" opacities kappa_abs, kappa_sca, kappa_ext cm^2 g^-1 ")')
write(*,'(" scattering matrix, see [App. B.1] sr^-1 /or/ cm^2 g^-1 sr^-1 ")')
write(*,'("")')
write(*,'("2 Examples")')
write(*,'("==========")')
write(*,'("")')
write(*,'(" A simple grain made only of the default pyroxene, for the default")')
write(*,'(" grain size distribution ($a^{-3.5}$ powerlaw from 0.05 to 3000mum), on")')
write(*,'(" the default wavelength grid (0.05mum to 1cm).")')
write(*,'("")')
write(*,'(" ,----")')
write(*,'(" | optool pyr")')
write(*,'(" ----")')
write(*,'("")')
write(*,'(" Include the scattering matrix in the produced output")')
write(*,'("")')
write(*,'(" ,----")')
write(*,'(" | optool pyr -s")')
write(*,'(" ----")')
write(*,'("")')
write(*,'(" Reproduce the DIANA standard dust model, using a specific pyroxene")')
write(*,'(" (70% Mg) and carbon, in a mass ratio 0.87/0.13, and with a porosity of")')
write(*,'(" 25%.")')
write(*,'("")')
write(*,'(" ,----")')
write(*,'(" | optool pyr-mg70 0.87 c 0.13 -p 0.25")')
write(*,'(" ----")')
write(*,'("")')
write(*,'(" List the built-in materials")')
write(*,'("")')
write(*,'(" ,----")')
write(*,'(" | optool -c")')
write(*,'(" ----")')
write(*,'("")')
write(*,'(" Add a water ice mantle (built-in data from Warren+08) that is 20% of")')
write(*,'(" the core mass")')
write(*,'("")')
write(*,'(" ,----")')
write(*,'(" | optool pyr-mg70 0.87 c 0.13 -m h2o-w 0.2 -p 0.25")')
write(*,'(" ----")')
write(*,'("")')
write(*,'(" Like the previous example, but use ice refractive index data from a")')
write(*,'(" separate file.")')
write(*,'("")')
write(*,'(" ,----")')
write(*,'(" | optool pyr-mg70 0.87 c 0.13 -p 0.25 -m data/ice_hudgins.dat 0.2")')
write(*,'(" ----")')
write(*,'("")')
write(*,'(" Pure water ice grains in a narrow size distribution from 1 to 3")')
write(*,'(" microns, with 15 sample sizes following an $f(a)\propto a^{-2.5}$")')
write(*,'(" powerlaw size distribution. Also, restrict the wavelength range to")')
write(*,'(" 10-100mum, and turn off DHS to get perfect spheres (Mie).")')
write(*,'("")')
write(*,'(" ,----")')
write(*,'(" | optool h2o -a 1 3 2.5 15 -l 10 100 -mie")')
write(*,'(" ----")')
write(*,'("")')
write(*,'(" Use a log-normal size distribution around 2 mum with sigma=0.7")')
write(*,'(" instead.")')
write(*,'(" ,----")')
write(*,'(" | optool h2o -a 0.1 30 2.0:0.7 -l 10 100 -mie")')
write(*,'(" ----")')
write(*,'("")')
write(*,'(" For silicon carbide, compute the opacity of a single grain size")')
write(*,'(" (2.5mum) at lambda=8.9mum.")')
write(*,'("")')
write(*,'(" ,----")')
write(*,'(" | optool -a 2.5 -l 8.9 sic")')
write(*,'(" ----")')
write(*,'("")')
write(*,'(" Represent the default dust model (DIANA, you also get this when you do")')
write(*,'(" not give any materials at all) in 42 grain sizes, and produce input")')
write(*,'(" files for RADMC-3D, one for each grain size, with full scattering")')
write(*,'(" matrix, chopping 3 degrees from the scattering peak.")')
write(*,'("")')
write(*,'(" ,----")')
write(*,'(" | optool -na 42 -d -s -radmc -chop 3")')
write(*,'(" ----")')
write(*,'("")')
write(*,'(" Use MMF to compute the opacities of dust aggregates made of pyroxene")')
write(*,'(" monomers. Use a monomer radius of 0.3 mum to construct aggregates")')
write(*,'(" with compact-volume radii between 10 and 30 mum, and a fractal")')
write(*,'(" dimension of 1.9.")')
write(*,'("")')
write(*,'(" ,----")')
write(*,'(" | optool pyr -a 10 30 -mmf 0.3 1.9")')
write(*,'(" ----")')
write(*,'("")')
write(*,'(" Compute CDE for small graphite grains")')
write(*,'(" ,----")')
write(*,'(" | optool gra -a 0.01 0.1 -l 1 30 -cde")')
write(*,'(" ----")')
write(*,'("")')
write(*,'("3 Installation")')
write(*,'("==============")')
write(*,'("")')
write(*,'(" You can download, compile, and install optool with these simple")')
write(*,'(" steps, using the freely available GNU FORTRAN compiler [ gfortran ].")')
write(*,'("")')
write(*,'(" ,----")')
write(*,'(" | git clone https://github.com/cdominik/optool.git # clone repository")')
write(*,'(" | cd optool # enter code directory")')
write(*,'(" | make multi=true # compile with multicore support")')
write(*,'(" ----")')
write(*,'("")')
write(*,'(" ,----")')
write(*,'(" | make install bindir=~/bin/ # optional: copy binaries to binary path")')
write(*,'(" | pip install -e . # optional: install the python module")')
write(*,'(" ----")')
write(*,'("")')
write(*,'(" In the compilation step, use multi=true to add multicore support")')
write(*,'(" (recommended!), ifort=true to use the [Intel fortran compiler],")')
write(*,'(" fits=true to support FITS files[2], and oldio=true if your")')
write(*,'(" compiler does not have the ISO_FORTRAN_ENV module.")')
write(*,'("")')
write(*,'(" The executable is called optool , run it with ./optool or move it")')
write(*,'(" onto your binary path.")')
write(*,'("")')
write(*,'("4 Command line arguments")')
write(*,'("========================")')
write(*,'("")')
endif
if((what.eq.'-h').or.(what.eq.'all')) then
write(*,'(" -h [OPT] ")')
write(*,'(" Show command line option summary, or specific information about")')
write(*,'(" option * -OPT *.")')
endif
if((what.eq.'-q').or.(what.eq.'all')) then
write(*,'(" -q ")')
write(*,'(" Reduce output to STDOUT to essential warnings and errors.")')
endif
if((what.eq.'-v').or.(what.eq.'all')) then
write(*,'(" -v ")')
write(*,'(" More verbose output to STDOUT.")')
write(*,'("")')
endif
if(what.eq.'all') then
write(*,'("4.1 Grain composition")')
write(*,'("~~~~~~~~~~~~~~~~~~~~~")')
write(*,'("")')
write(*,'(" If no composition is specified, the (DIANA) default is *-c pyr 0.87 -c")')
write(*,'(" c 0.13 -p 0.25*.")')
write(*,'("")')
endif
if((what.eq.'-c').or.(what.eq.'all')) then
write(*,'(" -c ")')
write(*,'(" List available built-in materials (the keys for the *-c* and")')
write(*,'(" *-m* options).")')
write(*,'("")')
endif
if((what.eq.'-c').or.(what.eq.'all')) then
write(*,'(" [-c] MATERIAL [MFRAC] ")')
write(*,'(" Specify a material to include in the grain. MATERIAL can be")')
write(*,'(" the [key for a builtin material], the [path to an lnk file],")')
write(*,'(" or colon-separated numbers n:k:rho [3]. MFRAC is the /mass/")')
write(*,'(" fraction (default 1.0) of the material. You can give up to 20")')
write(*,'(" materials to build up the grain. Mass fractions do not have to")')
write(*,'(" add up to one, they will be renormalized. All materials will be")')
write(*,'(" mixed together using the /Bruggeman/ rule, and vacuum can be")')
write(*,'(" added through the porosity. *-c* in front of each material is")')
write(*,'(" optional.")')
write(*,'("")')
endif
if((what.eq.'-m').or.(what.eq.'all')) then
write(*,'(" -m MATERIAL [MFRAC] ")')
write(*,'(" Like *-c*, but place this material into the grain")')
write(*,'(" mantle. Multiple mantle materials will be mixed using the")')
write(*,'(" Bruggeman rule, and then that mix will be added to the core")')
write(*,'(" using the /Maxwell-Garnett/ rule. The *-m* is /not/ optional,")')
write(*,'(" it must be present.")')
write(*,'("")')
endif
if((what.eq.'-p').or.(what.eq.'all')) then
write(*,'(" -p POROSITY [P_MANTLE] ")')
write(*,'(" Porosity, the /volume/ fraction of vacuum, a number smaller than")')
write(*,'(" 1. The default is 0. A single value will apply to both core")')
write(*,'(" and mantle, but a second value will be specific for the mantle")')
write(*,'(" (and may be 0).")')
write(*,'("")')
endif
if((what.eq.'-diana').or.(what.eq.'all')) then
write(*,'(" -diana, -dsharp , -dsharp-no-ice ")')
write(*,'(" Use DIANA (Woitke+2016) or DSHARP (Birnstiel+2018) compositions.")')
write(*,'("")')
endif
if(what.eq.'all') then
write(*,'("4.2 Grain geometry and computational method")')
write(*,'("~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~")')
write(*,'("")')
write(*,'(" If no method is explicitly specified, the default is *-dhs 0.8*, i.e.")')
write(*,'(" DHS with f_max=0.8.")')
write(*,'("")')
endif
if((what.eq.'-dhs').or.(what.eq.'all')) then
write(*,'(" -dhs [FMAX] ")')
write(*,'(" Use the /Distribution of Hollow Spheres/ (DHS, Min+ 2005)")')
write(*,'(" approach to model deviations from perfect spherical symmetry and")')
write(*,'(" low-porosity aggregates. Spheres with inner holes with volume")')
write(*,'(" fractions between 0 and f_max (default 0.8) are averaged to")')
write(*,'(" mimic irregularities. f_max=0 means to use solid spheres (Mie")')
write(*,'(" theory), i.e. perfectly regular grains. For backward")')
write(*,'(" compatibility, *-fmax* can be used instead of *-dhs*.")')
write(*,'("")')
endif
if((what.eq.'-mmf').or.(what.eq.'all')) then
write(*,'(" -mmf [A0 [DFRAC-OR-FILL [KF]]] ")')
write(*,'(" Use /Modified Mean Field/ theory (MMF, Tazaki & Tanaka 2018) to")')
write(*,'(" compute opacities of highly porous or fractal aggregates. *-c*,")')
write(*,'(" *-m*, and *-p* determine the composition of monomers with radius")')
write(*,'(" A0 (default 0.1mum). Particles will be aggregates with a")')
write(*,'(" /compact size/ given by the *-a* switch, giving rise to")')
write(*,'(" $N=a^3/a_0^3$ monomers. DFRAC-OR-FILL specifies either the")')
write(*,'(" fractal dimension (if >1) or the /volume filling factor/ (if")')
write(*,'(" <1). The default is 0.2. KF may be used to change the default")')
write(*,'(" prefactor.")')
write(*,'("")')
endif
if((what.eq.'-mie').or.(what.eq.'all')) then
write(*,'(" -mie ")')
write(*,'(" Do a standard Mie calculation for perfect spheres. This is")')
write(*,'(" short for *-dhs 0*.")')
write(*,'("")')
endif
if((what.eq.'-cde').or.(what.eq.'all')) then
write(*,'(" -cde ")')
write(*,'(" Compute CDE (continuous distribution of ellipsoids) Rayleigh")')
write(*,'(" limit opacities.")')
write(*,'("")')
endif
if(what.eq.'all') then
write(*,'("4.3 Grain size distribution")')
write(*,'("~~~~~~~~~~~~~~~~~~~~~~~~~~~")')
write(*,'("")')
endif
if((what.eq.'-a').or.(what.eq.'-amin').or.(what.eq.'-amax').or.(what.eq.'-apow').or.(what.eq.'-na').or.(what.eq.'all')&
.or.(what.eq.'-amean').or.(what.eq.'-asig')) then
write(*,'(" -a AMIN [AMAX [APOW [NA]]] (powerlaw size distribution)")')
write(*,'(" Specify (minimum) grain radius, and optionally maximum grain")')
write(*,'(" radius, the [size distribution powerlaw] and the number of size")')
write(*,'(" bins. You may also use options to set individual values with")')
write(*,'(" *-amin*, *-amax*, *-apow*, *-na*. The defaults are 0.05 mum,")')
write(*,'(" 3000 mum, 3.5, and /15 per size decade with a fixed minimum of")')
write(*,'(" 5/, respectively.")')
write(*,'(" > If only a single size is specified with *-a*, then")')
write(*,'(" a_max=a_min and n_a=1 are implied.")')
write(*,'("")')
endif
if((what.eq.'-a').or.(what.eq.'-amin').or.(what.eq.'-amax').or.(what.eq.'-apow').or.(what.eq.'-na').or.(what.eq.'all')&
.or.(what.eq.'-amean').or.(what.eq.'-asig')) then
write(*,'(" -a AMIN AMAX AMEAN:ASIG [NA] ([log-]normal size distribution)")')
write(*,'(" Specify the centroid size and the logarithmic width for a")')
write(*,'(" [log-normal size distribution]. You may also use *-amean* and")')
write(*,'(" *-asig* options to set these values. If ASIG is negative,")')
write(*,'(" create a [normal distribution] with that width (in mum) around")')
write(*,'(" AMEAN .")')
write(*,'("")')
endif
if((what.eq.'-a').or.(what.eq.'-amin').or.(what.eq.'-amax').or.(what.eq.'-apow').or.(what.eq.'-na').or.(what.eq.'all')&
.or.(what.eq.'-amean').or.(what.eq.'-asig')) then
write(*,'(" -a FILE ")')
write(*,'(" Read the size distribution from a file. The file format is")')
write(*,'(" described in [appendix A]. To get an example file")')
write(*,'(" optool_sd.dat , run optool with the option *-w*.")')
write(*,'("")')
endif
if(what.eq.'all') then
write(*,'("4.4 Wavelength grid")')
write(*,'("~~~~~~~~~~~~~~~~~~~")')
write(*,'("")')
endif
if((what.eq.'-l').or.(what.eq.'-lmin').or.(what.eq.'-lmax').or.(what.eq.'-nl').or.(what.eq.'-nlam').or.(what.eq.'all')) then
write(*,'(" -l LMIN [LMAX [NLAM]] ")')
write(*,'(" Specify the (minimum) wavelength, and optionally the maximum")')
write(*,'(" wavelength and the number of wavelengths points for the")')
write(*,'(" construction of the wavelength grid. The default values are")')
write(*,'(" 0.05 mum, 10000 mum, and 300, respectively. You may also use")')
write(*,'(" the options *-lmin*, *-lmax*, and *-nlam* (or *-nl*) to set")')
write(*,'(" individual values.")')
write(*,'(" > If only one wavelength is specified with *-l*, then")')
write(*,'(" lambda_max=lambda_min and n_lambda=1 are implied.")')
write(*,'("")')
endif
if((what.eq.'-l').or.(what.eq.'-lmin').or.(what.eq.'-lmax').or.(what.eq.'-nl').or.(what.eq.'-nlam').or.(what.eq.'all')) then
write(*,'(" -l FILE ")')
write(*,'(" Read the wavelength grid from FILE . To get an example file")')
write(*,'(" optool_lam.dat , run optool with the option *-w*. An [ lnk ]")')
write(*,'(" file could be used here as well!")')
write(*,'("")')
endif
if(what.eq.'all') then
write(*,'("4.5 Controlling the output")')
write(*,'("~~~~~~~~~~~~~~~~~~~~~~~~~~")')
write(*,'("")')
write(*,'(" The standard output is the file [ dustkappa.dat ], with the opacities")')
write(*,'(" and the asymmetry parameter /g/. The following options control and")')
write(*,'(" extend the [output].")')
write(*,'("")')
endif
if((what.eq.'-o').or.(what.eq.'all')) then
write(*,'(" -o [DIR] ")')
write(*,'(" Put the output files in directory DIR instead of the current")')
write(*,'(" working directory. ./output will be used if *-o* is present")')
write(*,'(" but DIR is not specified.")')
write(*,'("")')
endif
if((what.eq.'-s').or.(what.eq.'all')) then
write(*,'(" -s [NANG] ")')
write(*,'(" Include the scattering matrix in the output. NANG may optionally")')
write(*,'(" change the number of equally-spaced [angular grid points] to")')
write(*,'(" cover the range of angles between 0 and 180 degrees. The")')
write(*,'(" default for NANG is 180 and should normally be just fine.")')
write(*,'("")')
endif
if((what.eq.'-d').or.(what.eq.'all')) then
write(*,'(" -d [NSUB] ")')
write(*,'(" Divide the computation up into n_a parts to produce a file for")')
write(*,'(" each grain size. Each size will be an average over a range of")')
write(*,'(" NSUB (default 5) grains around the real size.")')
write(*,'("")')
endif
if((what.eq.'-chop').or.(what.eq.'all')) then
write(*,'(" -chop [NDEG] ")')
write(*,'(" Cap the first NDEG (2 if unspecified) degrees of the [forward")')
write(*,'(" scattering peak].")')
write(*,'("")')
endif
if((what.eq.'-fits').or.(what.eq.'all')) then
write(*,'(" -fits ")')
write(*,'(" Write [ dustkappa.fits ] instead of ASCII output. With -d ,")')
write(*,'(" write n_a files.")')
write(*,'("")')
endif
if((what.eq.'-radmc').or.(what.eq.'all')) then
write(*,'(" -radmc [LABEL] ")')
write(*,'(" RADMC-3D uses a different angular grid and [scattering matrix]")')
write(*,'(" normalization. File names will contain LABEL if specified and")')
write(*,'(" have the extension .inp .")')
write(*,'("")')
endif
if((what.eq.'-print').or.(what.eq.'all')) then
write(*,'(" -print [KEY] ")')
write(*,'(" Write to STDOUT instead of files. The default is to write")')
write(*,'(" lambda, kappa_abs, kappa_sca, kappa_ext, and g. Many other")')
write(*,'(" outputs are possible, run optool -print ? for a full list. For")')
write(*,'(" readability, a header line may be printed to STDERR , but")')
write(*,'(" STDOUT gets only numbers which can be used in pipes and for")')
write(*,'(" redirection. You can use this to extract a single value, for")')
write(*,'(" example the 850mum extinction opacity of grains between 1 and")')
write(*,'(" 3mm: optool -a 1000 3000 -l 850 -print kext .")')
write(*,'("")')
endif
if((what.eq.'-w').or.(what.eq.'all')) then
write(*,'(" -w ")')
write(*,'(" Write the files optool_sd.dat and optool_lam.dat with the")')
write(*,'(" grain size distribution and the wavelength grid,")')
write(*,'(" respectively. Also, write optool_mix.lnk with the result of")')
write(*,'(" mixing refractive index data. Exit without doing a computation.")')
write(*,'("")')
endif
if(what.eq.'all') then
write(*,'("5 Material properties")')
write(*,'("=====================")')
write(*,'("")')
write(*,'(" optool needs refractive index data to work. For your convenience, a")')
write(*,'(" useful list of materials is compiled into optool . You can also find")')
write(*,'(" and use other data.")')
write(*,'("")')
write(*,'("5.1 Built-in materials")')
write(*,'("~~~~~~~~~~~~~~~~~~~~~~")')
write(*,'("")')
write(*,'(" To access one of the built-in materials, specify the corresponding key")')
write(*,'(" string like pyr-mg70 . In each material class we have selected a")')
write(*,'(" useful default, accessible with an even simpler generic key (for")')
write(*,'(" example, pyr is an alias for pyr-mg70 ). Most of the built-in")')
write(*,'(" refractive index datasets have a reasonably wide wavelength coverage -")')
write(*,'(" the few exceptions are highlighted by bold-face numbers. If a")')
write(*,'(" material is being used outside of the measured region, optool will")')
write(*,'(" still function, using extrapolated optical properties.")')
write(*,'("")')
write(*,'(" Even the limited number of materials we have selected to include with")')
write(*,'(" optool can be daunting. To get started with some kind of standard")')
write(*,'(" opacity, we recommend to work with pyroxene \fbox{pyr}, carbon")')
write(*,'(" \fbox{c}, and, at low temperatures, water ice \fbox{h2o} (Woitke+")')
write(*,'(" 2016). If you need to account for sulfur, you may want to include")')
write(*,'(" troilite \fbox{tro} (Birnstiel+ 2016).")')
write(*,'("")')
write(*,'(" -c Key -c Key Material State Reference ")')
write(*,'(" generic full key ")')
write(*,'("---------------------------------------------------------- --------------")')
write(*,'(" pyr-mg100 MgSiO_3 amorph [Dorschner+95")')
write(*,'(" pyr-mg95 Mg_0.95 Fe_0.05 SiO_3 amorph [Dorschner+95")')
write(*,'(" pyr-mg80 Mg_0.8 Fe_0.2 SiO_3 amorph [Dorschner+95")')
write(*,'(" pyr pyr-mg70 Mg_0.7 Fe_0.3 SiO_3 amorph [Dorschner+95")')
write(*,'(" pyr-mg60 Mg_0.6 Fe_0.4 SiO_3 amorph [Dorschner+95")')
write(*,'(" pyr-mg50 Mg_0.5 Fe_0.5 SiO_3 amorph [Dorschner+95")')
write(*,'(" pyr-mg40 Mg_0.4 Fe_0.6 SiO_3 amorph [Dorschner+95")')
write(*,'(" ens pyr-c-mg96 Mg_0.96 Fe_0.04 SiO3 cryst[4 [Jäger+98")')
write(*,'("---------------------------------------------------------- --------------")')
write(*,'(" ol ol-mg50 MgFeSiO_4 amorph [Dorschner+95")')
write(*,'(" ol-mg40 Mg_0.8 Fe_1.2 SiO_4 amorph [Dorschner+95")')
write(*,'(" for ol-c-mg100 Mg_2 SiO_4 cryst[4 [Suto+06")')
write(*,'(" ol-c-mg95 Mg_1.9 Fe_0.1 SiO_4 cryst[4 [Fabian+01")')
write(*,'(" fay ol-c-mg00 Fe_2 SiO_4 cryst[4 [Fabian+01")')
write(*,'("---------------------------------------------------------- --------------")')
write(*,'(" astrosil MgFeSiO_4 mixed [Draine+03")')
write(*,'("---------------------------------------------------------- --------------")')
write(*,'(" c c-z C amorph? [Zubko+96")')
write(*,'(" c-p C amorph [Preibisch+93")')
write(*,'(" gra c-gra C graphite cryst[4 [Draine+03")')
write(*,'(" org c-org CHON organics amorph [Henning+96")')
write(*,'(" c-nano C nano-diamond cryst [Mutschke+04")')
write(*,'("---------------------------------------------------------- --------------")')
write(*,'(" iron fe-c Fe metal [Henning+96")')
write(*,'(" tro fes FeS metal [Henning+96")')
write(*,'(" sic SiC cryst [Laor93")')
write(*,'("---------------------------------------------------------- --------------")')
write(*,'(" qua sio2 SiO_2 amorph [Kitamura+07")')
write(*,'(" cor cor-c Al_2 O_3 cryst [Koike+95")')
write(*,'("---------------------------------------------------------- --------------")')
write(*,'(" h2o h2o-w Water ice cryst [Warren+08")')
write(*,'(" h2o-a Water ice amorph [Hudgins+93")')
write(*,'(" co2 co2-w CO_2 ice cryst [Warren+86")')
write(*,'(" nh3 nh3-m NH_3 ice cryst [Martonchik+83")')
write(*,'(" co co-a CO ice amorph [Palumbo+06")')
write(*,'(" co2-a / c CO_2 ice am / cr [Gerakines+20")')
write(*,'(" ch4-a / c CH_4 ice am / cr [Gerakines+20")')
write(*,'(" ch3oh-a / c CH_3 OH ice am / cr [Gerakines+20")')
write(*,'("")')
write(*,'("5.2 External refractory index files ( lnk files)")')
write(*,'("~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~")')
write(*,'("")')
write(*,'(" optool can use external refractive index data in files with the")')
write(*,'(" following format[5]:")')
write(*,'(" - The file may start with several comment lines (lines starting with")')
write(*,'(" ! , # , or * ).")')
write(*,'(" - The next line contains two numbers, the number of wavelengths")')
write(*,'(" $n_\lambda$ and the specific density rho of the material in")')
write(*,'(" g/cm^{3}.")')
write(*,'(" - The remaining lines should form three columns of data: lambda[mum]")')
write(*,'(" (sorted either up or down), and the real and imaginary parts of the")')
write(*,'(" refractive index, $n$ and $k$.")')
write(*,'("")')
write(*,'(" We provide additional data ready for use with optool in [a separate")')
write(*,'(" repository]. Other resources are the [Jena database], [ARIA] and")')
write(*,'(" original papers in the literature. Don t forget to add the line with")')
write(*,'(" $n_\lambda$ and rho! If that is not possible, optool will count the")')
write(*,'(" lines and you can specify the density after the mass fraction, like")')
write(*,'(" this: optool -c path/to/file.lnk 0.7 3.42 . Please include")')
write(*,'(" references for any optical properties used in your study.")')
write(*,'("")')
write(*,'("5.3 One-off materials")')
write(*,'("~~~~~~~~~~~~~~~~~~~~~")')
write(*,'("")')
write(*,'(" For a calculation at a single wavelength you can give the refractive")')
write(*,'(" index on the command line, like this: optool 1.57:0.56:2.08 -l 0.74")')
write(*,'(" -s . This example specifies the refractive index of $m=1.57+0.56i$ for")')
write(*,'(" a material with a density of 2.08g/cm^3, and the computation of the")')
write(*,'(" scattering matrix will be done at a wavelength of 0.74mum.")')
write(*,'("")')
write(*,'("6 Output files")')
write(*,'("==============")')
write(*,'("")')
write(*,'(" dustkappa.dat")')
write(*,'(" This is an ASCII file containing the basic opacity results. It")')
write(*,'(" starts with a comment section describing the dust model and also")')
write(*,'(" showing the exact command line that was used to produce the")')
write(*,'(" file. The header is followed by the format number (3,")')
write(*,'(" currently), followed by the number of wavelengths in the grid,")')
write(*,'(" both on lines by themselves. This is followed by a block with")')
write(*,'(" these columns:")')
write(*,'("")')
write(*,'(" 1. wavelength lambda [micron]")')
write(*,'(" 2. mass absorption cross section kappa_abs [cm^2/g]")')
write(*,'(" 3. mass scattering cross section kappa_sca [cm^2/g]")')
write(*,'(" 4. asymmetry parameter /g/")')
write(*,'("")')
write(*,'(" dustkapscatmat.dat")')
write(*,'(" ASCII file with cross sections and full scattering matrix. It is")')
write(*,'(" an extended version of the dustkappa.dat file. This file has")')
write(*,'(" a format number (0), the number of wavelengths and then the")')
write(*,'(" number of angular points after the comment section. After an")')
write(*,'(" empty line, the same opacity block as in dustkappa.dat is")')
write(*,'(" present. Another empty line is followed by a list of the grid")')
write(*,'(" angles, another empty line, and then the scattering matrix")')
write(*,'(" elements for all wavelengths and all angles. The comment section")')
write(*,'(" at the start of the file shows the structure in a formal way.")')
write(*,'(" See [appendix B.1] for information about the normalization of")')
write(*,'(" the scattering matrix and about the angular grid that is used")')
write(*,'(" for it. Also, see the -radmc switch which will modify[6] the")')
write(*,'(" output to make sure it can be used as an input file for")')
write(*,'(" [RADMC-3D].")')
write(*,'("")')
write(*,'(" To save space, optool can write a /sparse file/ (iformat=100)")')
write(*,'(" that stores the full scattering matrix only for selected")')
write(*,'(" wavelengths (for example, the ones that will be used for image")')
write(*,'(" generation). Use -sp LAM or -sp LAM1 LAM2 to define a")')
write(*,'(" wavelength (interval)[7] for the matrix to be stored. Multiple")')
endif
if((what.eq.'-sp').or.(what.eq.'all')) then
write(*,'(" -sp switches are allowed.")')
write(*,'("")')
write(*,'(" dustkappa.fits")')
write(*,'(" The FITS-file is written when using the -fits switch. It has")')
write(*,'(" two HDU blocks. The first contains the cross sections per unit")')
write(*,'(" mass (units cm^2/g). It is a n_lambda * 4 matrix with these")')
write(*,'(" columns: wavelength in micron, kappa_ext, kappa_abs, kappa_sca.")')
write(*,'(" The second block is a n_lambda * 6 * n_ang matrix, containing")')
write(*,'(" the 6 elements of the scattering matrix (F_11, F_12, F_22, F_33,")')
write(*,'(" F_34, and F_44) for n_ang equidistant scattering angles from")')
write(*,'(" forward scattering (element 0) to backward scattering (element")')
write(*,'(" n_ang-1), for each lambda.")')
write(*,'("")')
write(*,'(" optool.tex")')
write(*,'(" As a little gimmick, you can run optool2tex with the exact")')
write(*,'(" same command line arguments as used in an optool ")')
write(*,'(" run. optool.tex then contains text and a table, describing the")')
write(*,'(" methods used for the opacity computation and listing the")')
write(*,'(" composition of the grains. All relevant references are given -")')
write(*,'(" the BibTeX file optool.bib is required for the file to be")')
write(*,'(" processed properly. You can rework this text to include it into")')
write(*,'(" your paper. For more details, read the comment section in")')
write(*,'(" optool2tex .")')
write(*,'("")')
write(*,'(" optool_mix.lnk")')
write(*,'(" When using the -w switch, optool will write the result of")')
write(*,'(" mixing to this file.")')
write(*,'("")')
write(*,'(" optool_sd.dat & optool_lam.dat")')
write(*,'(" When using the -w switch, optool will write the grain size")')
write(*,'(" grid and the wavelength grid to two files. The files serve as")')
write(*,'(" example files for what the structure of files need to be to be")')
write(*,'(" read in with -a file.dat and -l file.dat for user-provided")')
write(*,'(" size and wavelength grids.")')
write(*,'("")')
endif
if(what.eq.'all') then
write(*,'("7 Python interface")')
write(*,'("==================")')
write(*,'("")')
write(*,'(" optool comes with a [ python ] module[8] optool.py that runs")')
write(*,'(" optool in the background[9] and puts all computed quantities as")')
write(*,'(" numpy arrays into a python object. This makes it straight forward")')
write(*,'(" to inspect and further process the output. Here is how to use it:")')
write(*,'("")')
write(*,'(" ,----")')
write(*,'(" | import optool")')
write(*,'(" | p = optool.particle(QQ~/bin/optool pyr 0.8 -m h2o 0.2 -na 24 -dQQ)")')
write(*,'(" ----")')
write(*,'("")')
write(*,'(" The argument to optool.particle() must be a valid shell command[10]")')
write(*,'(" to run optool , if necessary with the full path to the optool ")')
write(*,'(" binary. Depending on the presence of the optool s *-d* switch, the")')
write(*,'(" command will produce opacities either for $n_p=1$ particle, or for")')
write(*,'(" $n_p=n_a$ particles. Most of the attributes (with the exception of the")')
write(*,'(" global wavelength and angular grids) will therefore be arrays with the")')
write(*,'(" first dimension equal to $n_p$, even if $n_p=1$. The resulting object")')
write(*,'(" will have the following attributes:")')
write(*,'("")')
write(*,'(" *Attribute* *Type/Shape* *Quantity* ")')
write(*,'(" --------------------------------------------------------------------------------------------")')
write(*,'(" cmd string The full command given in the particle() call ")')
write(*,'(" --------------------------------------------------------------------------------------------")')
write(*,'(" radmc boolean Output follows RADMC conventions ")')
write(*,'(" scat boolean Scattering matrix is available ")')
write(*,'(" --------------------------------------------------------------------------------------------")')
write(*,'(" nlam int Number of wavelength points ")')
write(*,'(" lam float[nlam] The wavelength grid ")')
write(*,'(" nang int Number of scattering angles ")')
write(*,'(" scatang float[nang] The angular grid ")')
write(*,'(" --------------------------------------------------------------------------------------------")')
write(*,'(" materials [[[...]...]... ] Lists with [location,m_{frac},rho,material] ")')
write(*,'(" --------------------------------------------------------------------------------------------")')
write(*,'(" np int Number of particles, either 1 or (with -d) n_a ")')
write(*,'(" --------------------------------------------------------------------------------------------")')
write(*,'(" fmax float[np] Maximum volume fraction of vacuum for DHS ")')
write(*,'(" pcore , pmantle float[np] Porosity of the core/mantle material ")')
write(*,'(" --------------------------------------------------------------------------------------------")')
write(*,'(" amin , amax float[np] min/max grain size used for each particle ")')
write(*,'(" nsub int[np] Number of sizes averaged for each particle ")')
write(*,'(" apow float[np] Negative size distribution power law (e.g. 3.5) ")')
write(*,'(" amean , asig float[np] Centroid & width of (log-)normal distrbution ")')
write(*,'(" a1 , a2 , a3 float[np] Mean <a>, $\sqrt{<a^2>}$, and $\sqrt[3]{<a^3>}$ ")')
write(*,'(" rho float[np] Specific density of grains ")')
write(*,'(" --------------------------------------------------------------------------------------------")')
write(*,'(" kabs,ksca,kext float[np,nlam] Absorption,scattering,extinction cross section ")')
write(*,'(" gsca float[np,nlam] Asymmetry parameter ")')
write(*,'(" --------------------------------------------------------------------------------------------")')
write(*,'(" f11 , ..., f44 float[np,nlam,nang] Scattering matrix element F_11, ... ,F_44 ")')
write(*,'(" chop float[np] Degrees chopped off forward scattering ")')
write(*,'(" --------------------------------------------------------------------------------------------")')
write(*,'(" plot() method Plot the cross sections and matrix elements ")')
write(*,'(" --------------------------------------------------------------------------------------------")')
write(*,'(" computemean() method Compute Planck/Rosseland mean opacities ")')
write(*,'(" tmin,tmax,ntemp float,float,int Temperature grid for mean opacities ")')
write(*,'(" temp float[ntemp] Temperatures used for mean opacities ")')
write(*,'(" kplanck,kross float[np,ntemp] Mean opacities, after calling computemean() ")')
write(*,'(" --------------------------------------------------------------------------------------------")')
write(*,'(" norm string Current scattering matrix normalization ")')
write(*,'(" scatnorm() method Check/change scat. matrix normalization ")')
write(*,'(" --------------------------------------------------------------------------------------------")')
write(*,'(" sizedist() method Sum opacities over a size distribution ")')
write(*,'("")')
write(*,'(" Applying the plot() method to a particle object like p.plot() ")')
write(*,'(" will show (see Fig 1):")')
write(*,'(" - a plot showing the opacities kappa_abs, kappa_sca, and kappa_ext as")')
write(*,'(" a function of wavelength, along with the asymmetry parameter /g/ (on")')
write(*,'(" a linear y-scale). Note that the blue /g/ curve does not have its")')
write(*,'(" own axis, imagine the full /y/ axis going from 0 to 1 for /g/.")')
write(*,'(" - a plot showing the scattering matrix elements as a function of")')
write(*,'(" scattering angle, with sliders to go through grain sizes and")')
write(*,'(" wavelengths. When interpreting the y axis, note that we plot the")')
write(*,'(" positive/negative $\log_{10}$ of positive/negative matrix elements,")')
write(*,'(" compressing the range from $10^{-2}$ to $10^2$ into a line (use the")')
write(*,'(" grey lines as a guide, ignore the y-axis labels). If you cannot see")')
write(*,'(" F_11, it is because it is equal to and hidden behind F_22. If you")')
write(*,'(" cannot see F_33, it is because it is equal to and hidden behind")')
write(*,'(" F_44.")')
write(*,'(" - If the computemean method has been called first, the mean")')
write(*,'(" opacities kappa_Planck and kappa_Ross are shown in a separate plot.")')
write(*,'(" The mean opacities are per unit of grain mass, so please apply a")')
write(*,'(" dust-to-gas mass ratio to obtain opacities for a gas-dust mixture.")')
write(*,'("")')
write(*,'(" <./maint/inspect.png>")')
write(*,'("")')
write(*,'(" The python module has a few more tricks up its sleeve (for details")')
write(*,'(" check the documentation inside the Python module file optool.py ):")')
write(*,'("")')
write(*,'(" - You can cache results of an expensive computation for quick")')
write(*,'(" reloading in a new python session. When running the following")')
write(*,'(" command twice, optool will be called only the first time. The")')
write(*,'(" cache will reside in the specified directory.")')
write(*,'(" ,----")')
write(*,'(" | sil = optool.particle(QQoptool -d ol-mg50 -na 100QQ,cache=QQsilcacheQQ)")')
write(*,'(" ----")')
write(*,'(" - You can read the results directory created by another run of optool")')
write(*,'(" by leaving the command parameter empty and specifying the directory")')
write(*,'(" as the second parameter.")')
write(*,'(" ,----")')
write(*,'(" | part = optool.particle(QQQQ,QQpath/to/directoryQQ)")')
write(*,'(" ----")')
write(*,'(" - A lnktable class to read, plot, modify and write lnk files.")')
write(*,'(" ,----")')
write(*,'(" | x = optool.lnktable(QQlnk_data/sio2-Kitamura2007.lnkQQ)")')
write(*,'(" | x.plot()")')
write(*,'(" ----")')
write(*,'("")')
write(*,'(" - Compute Planck and Rosseland mean opacities")')
write(*,'(" ,----")')
write(*,'(" | p = optool.particle(QQoptool pyr 0.87 c 0.13 -p 0.25QQ)")')
write(*,'(" | p.computemean(tmin=10.,tmax=1500.,ntemp=300)")')
write(*,'(" ----")')
write(*,'("")')
write(*,'(" - /Particle arithmetic/: multiplying optool.particle objects with")')
write(*,'(" factors and adding them, or applying size distributions to a")')
write(*,'(" pre-computed set of opacities. See the documentation of the optool")')
write(*,'(" module and [appendix C.1] for an example.")')
write(*,'("")')
write(*,'(" \appendix")')
write(*,'("")')
write(*,'("8 Size distribution")')
write(*,'("===================")')
write(*,'("")')
write(*,'(" optool implements powerlaw, log-normal, and normal size")')
write(*,'(" distributions. Each of these will be subject to a minimum and a")')
write(*,'(" maximum grain size. The grain size grid is logarithmic, so $da\propto")')
write(*,'(" a$. The logarithmic bins are then filled according to:")')
write(*,'(" powerlaw n(a) ~ a^{-p+1}")')
write(*,'(" log-normal distribution, triggered by sig>0 n(a) ~ exp[-0.5(\frac{ln (a/a_m)}{sig})^2]")')
write(*,'(" normal distribution[11], triggered by sig<0 n(a) ~ exp[-0.5(\frac{a-a_m}{sig})^2]")')
write(*,'("")')
write(*,'(" Other size distributions can be constructed using the [python")')
write(*,'(" interface]. Finally, optool can also read a size distribution from")')
write(*,'(" a file, and this is also the way to provide an arbitrary size")')
write(*,'(" grid. The first data line in the file gives the number of grain size")')
write(*,'(" bins, followed by lines with two numbers each: grain size in micron")')
write(*,'(" and number of grains in the corresponding bin. To get an example file,")')
write(*,'(" run optool with the option *-w*):")')
write(*,'("")')
write(*,'("9 Scattering Matrix: The fine print")')
write(*,'("===================================")')
write(*,'("")')
write(*,'("9.1 Phase function normalization")')
write(*,'("~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~")')
write(*,'("")')
write(*,'(" A number of different normalizations for the scattering matrix are")')
write(*,'(" being used in the literature and in computational tools. The")')
write(*,'(" differences are significant, and it is important to be aware of the")')
write(*,'(" choice. For optool we are using a convention ([Hovenier (2004)]) in")')
write(*,'(" which the average over all directions of the 1-1 element of the")')
write(*,'(" scattering matrix equals unity, i.e. the integral will be 4pi:")')
write(*,'("")')
write(*,'(" \begin{equation}")')
write(*,'(" \label{eq:1}")')
write(*,'(" \oint_{(4\pi)} F_{11}(\lambda,\Theta) d\Omega = ")')
write(*,'(" 2\pi \int_{-1}^{1} F_{11}(\lambda,\mu) {\rm d}\mu= 4\pi \quad ,")')
write(*,'(" \end{equation}")')
write(*,'("")')
write(*,'(" with $\mu=\cos\Theta$. optool can also produce output for [RADMC-3D]")')
write(*,'(" which uses instead")')
write(*,'("")')
write(*,'(" \begin{equation}")')
write(*,'(" \label{eq:2}")')
write(*,'(" \oint_{(4\pi)} Z_{11}(\lambda,\Theta) d\Omega =")')
write(*,'(" 2\pi \int_{-1}^{1} Z_{11}(\lambda,\mu) {\rm d}\mu =")')
write(*,'(" \kappa_{\rm sca}(\lambda) \quad .")')
write(*,'(" \end{equation}")')
write(*,'("")')
write(*,'(" The books by Bohren & Huffman and by Mishchenko use different")')
write(*,'(" normalizations again. You can change the normalization of the")')
write(*,'(" scattering matrix in the python interface with the scatnorm() ")')
write(*,'(" method. By default, that method checks the current normalization.")')
write(*,'(" Using an argument r , b , m , or h will modify the")')
write(*,'(" normalization.")')
write(*,'("")')
write(*,'("9.2 Forward-scattering peak")')
write(*,'("~~~~~~~~~~~~~~~~~~~~~~~~~~~")')
write(*,'("")')
write(*,'(" Particles that are much larger than the wavelength of the considered")')
write(*,'(" radiation can show extreme forward scattering, where much of the")')
write(*,'(" /scattered/ radiation is sent into just a few degrees around the")')
write(*,'(" forward direction. This can be difficult to handle for radiative")')
write(*,'(" transfer codes which have limited angular resolution or limited")')
write(*,'(" sampling. [MCMax3D] has the nspike keyword to deal with this")')
write(*,'(" issue. Other tools (e.g. RADMC-3D) require this to be taken care of by")')
write(*,'(" the process that creates the opacity files. The -chop switch")')
write(*,'(" specifies a number of degrees around the forward scattering")')
write(*,'(" direction. Inside that cone, the scattering matrix gets limited to the")')
write(*,'(" value at the edge of the cone. To compensate and ensure energy")')
write(*,'(" conservation, the scattering cross section will be reduced")')
write(*,'(" accordingly. As a result, the radiation that would be /scattered/ into")')
write(*,'(" this narrow range of angles will be treated as if it did have /no")')
write(*,'(" interaction at all/ with the grain.")')
write(*,'("")')
write(*,'("9.3 Angular grid")')
write(*,'("~~~~~~~~~~~~~~~~")')
write(*,'("")')
write(*,'(" optool uses an angular grid in one-degree steps from 0 to 180")')
write(*,'(" degrees. The full degrees are the cell /interfaces/ of that")')
write(*,'(" grid. optool computes the scattering matrix at the cell /midpoints/,")')
write(*,'(" i.e. at 0.5degree, 1.5degree etc to 179.5degree, for a total of 180")')
write(*,'(" values. The scattering matrix is normalized in this way, so that a")')
write(*,'(" numerical integral gives the correct result.")')
write(*,'("")')
write(*,'(" RADMC-3D requires the values of the scattering matrix on the cell")')
write(*,'(" /boundaries/, so at 0degree, 1degree etc to 180degree, for a total of")')
write(*,'(" 181 values. For the input files for RADMC-3D, we interpolate and")')
write(*,'(" extend the computed values to the cell boundaries.")')
write(*,'("")')
write(*,'("10 More on optical properties")')
write(*,'("=============================")')
write(*,'("")')
write(*,'("10.1 Crystalline materials")')
write(*,'("~~~~~~~~~~~~~~~~~~~~~~~~~~")')
write(*,'("")')
write(*,'(" Crystalline materials have optical properties that are dependent on")')
write(*,'(" the relative orientation of the electric field of the incoming")')
write(*,'(" radiation and the crystallographic axes of the material. A fully")')
write(*,'(" correct treatment would require the use of the refractive index")')
write(*,'(" /tensor/ which is not implemented in optool (DDA can in principle do")')
write(*,'(" that). There are two approximations that can be used. For the")')
write(*,'(" materials built into optool we have assumed that the material")')
write(*,'(" consists of /many small crystalline areas that are randomly oriented")')
write(*,'(" within each grain/. For this, the refractive index data for different")')
write(*,'(" orientations have been combined using the Bruggeman effective medium")')
write(*,'(" rule. It results in a single refractive index data set. However, we")')
write(*,'(" can also think of a situation where /each dust grain is a single")')
write(*,'(" crystal in a cloud of randomly oriented grains/. In that case, we need")')
write(*,'(" to compute opacities for individual orientations, and then average the")')
write(*,'(" opacities. Axis-specific data is (when available) also included in the")')
write(*,'(" optool distribution, in the lnk_data/ad directory. You could use the")')
write(*,'(" [python interface] to do this kind of mixing in the following way:")')
write(*,'(" ,----")')
write(*,'(" | import optool")')
write(*,'(" | px = optool.particle(QQoptool -a 1.5 lnk_data/ad/c-gra-x-Draine2003.lnkQQ)")')
write(*,'(" | py = optool.particle(QQoptool -a 1.5 lnk_data/ad/c-gra-y-Draine2003.lnkQQ)")')
write(*,'(" | pz = optool.particle(QQoptool -a 1.5 lnk_data/ad/c-gra-z-Draine2003.lnkQQ)")')
write(*,'(" | pmix = (px+py+pz)/3.")')
write(*,'(" ----")')
write(*,'("")')
write(*,'("10.2 How to ingest refractive index data for another material")')
write(*,'("~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~")')
write(*,'("")')
write(*,'(" Using external refractive index data means that you have to keep track")')
write(*,'(" of where those files are. It can be convenient to compile your")')
write(*,'(" favorite materials into optool , so that accessing them will be as")')
write(*,'(" simple as using the [built-in materials]. Here is how to do that:")')
write(*,'("")')
write(*,'(" 1. Give your lnk file a name exactly like")')
write(*,'(" pyr-mg70-Dorschner1995.lnk , where the start of the name")')
write(*,'(" ( pyr-mg70 ) is the key to access the material and Dorschner1995 ")')
write(*,'(" (the text after the final - ) is the reference.")')
write(*,'(" 2. Put this file into the lnk_data directory.")')
write(*,'(" 3. Optionally edit lnk_data/lnk-help.txt , so that [ optool -c ] will")')
write(*,'(" list the new material. Note that, in order to define /generic")')
write(*,'(" keys/, optool looks for pairs that look like genkey -> fullkey in")')
write(*,'(" this file.")')
write(*,'(" 4. Run make ingest to update optool_refind.f90 , now with your new")')
write(*,'(" material.")')
write(*,'(" 5. Recompile and install the code.")')
write(*,'("")')
write(*,'("10.3 Overview of optical properties")')
write(*,'("~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~")')
write(*,'("")')
write(*,'(" The grid plot in Fig. 2 shows the imaginary parts of all built-in")')
write(*,'(" materials, in the wavelength range from 0.05 to 300 mum. Some if the")')
write(*,'(" ices have only data in a small range, where the vibrational")')
write(*,'(" transitions lie. However, these materials can be used over a much")')
write(*,'(" broader wavelength range, because the extrapolation becomes")')
write(*,'(" problematic only in the UV where electronic transitions kick in.")')
write(*,'("")')
write(*,'(" <./maint/all_k.pdf>")')
write(*,'("")')
write(*,'("11 Internals")')
write(*,'("============")')
write(*,'("")')
write(*,'(" This appendix describes some key aspects of the internal workings of")')
write(*,'(" the code.")')
write(*,'("")')
write(*,'(" Refractive index data")')
write(*,'(" Measured refractive index data is obtained from data compiled")')
write(*,'(" into the code, or read-in from a file. That data is then")')
write(*,'(" interpolated and extrapolated onto the wavelength grid requested")')
write(*,'(" for the computation. Extrapolation toward short wavelengths is")')
write(*,'(" done keeping the refractive indices constant. Extrapolation")')
write(*,'(" toward long wavelengths assumes that the last two measured data")')
write(*,'(" points define a powerlaw. Interpolation in the measured grid is")')
write(*,'(" done using double-logarithmic interpolation.")')
write(*,'("")')
write(*,'(" Mixing")')
write(*,'(" Once the refractive index for all involved materials is")')
write(*,'(" available, the core and the mantle mixtures are created")')
write(*,'(" independently, using the Bruggeman rule. Mass fractions are")')
write(*,'(" converted into volume fractions, and porosity is implemented")')
write(*,'(" using vacuum as an additional material. The subroutine doing")')
write(*,'(" the mixing uses an iterative procedure that is very stable, also")')
write(*,'(" for a large number of components.")')
write(*,'(" If there is a mantle, the Maxwell-Garnett rule is applied with")')
write(*,'(" the core being treated as an inclusion inside a mantle matrix.")')
write(*,'(" With a -w switch, optool will write result of mixing into")')
write(*,'(" the file optool_mix.lnk .")')
write(*,'("")')
write(*,'(" DHS")')
write(*,'(" In order to simulate irregularities in grains (irregular shapes,")')
write(*,'(" or the properties of low-porosity aggregates), optool averages")')
write(*,'(" the opacities of grains with an inner empty region, over a range")')
write(*,'(" of volume fractions of this inner region between 0 and $f_{\rm")')
write(*,'(" max}$. The subroutine used to compute the opacities and")')
write(*,'(" scattering matrix elements for these structures is DMiLay ")')
write(*,'(" (Toon & Ackerman 1981). For speed, you can use -xlim 1e3 or so")')
write(*,'(" to set a limit for the size parameter ($x=2\pi a/\lambda$) where")')
write(*,'(" optool switches from DHS to Mie[12]. In rare cases, DMiLay ")')
write(*,'(" might not converge and return an error. optool then falls back")')
write(*,'(" to a Mie computation, using the routine MeerhoffMie (de Rooij+")')
write(*,'(" 1984), with a size parameter limit to ensure stability. In a")')
write(*,'(" situation where the imaginary part of the refractive index is")')
write(*,'(" extremely small, numerical inaccuracies may lead to an")')
write(*,'(" unphysical result with $q_{\rm sca}>q_{\rm ext}$, implying a")')
write(*,'(" small negative $q_{\rm abs}$. optool therefore enforces")')
write(*,'(" $q_{\rm abs}>q_{\rm ext}/10^4$.")')
write(*,'("")')
write(*,'(" MMF")')
write(*,'(" To construct fluffy/fractal aggregates, optool needs the")')
write(*,'(" number of monomers $N$, the fractal dimension $D_{\rm f}$, and a")')
write(*,'(" scaling factor $k_{\rm f}$ which are related to the radius of")')
write(*,'(" gyration $R_{\rm g}$ of the aggregate by $N=k_{\rm f}(R_{\rm")')
write(*,'(" g}/a_0)^{D_{\rm f}}$. The size $a$ of the particles as")')
write(*,'(" specified by the *-a* switch is interpreted as the /compact/[13]")')
write(*,'(" size of all material in the aggregate, so that $N=a^3/a_0^3$,")')
write(*,'(" where $a_0$ is the monomer radius. The average volume filling")')
write(*,'(" factor $f$ can be expressed by $f=N\cdot(\sqrt{3/5}\,a_0/R_{\rm")')
write(*,'(" g})^3$. To determine the structure of the aggregates, the user")')
write(*,'(" can specify a structure parameter. If that parameter is >1, it")')
write(*,'(" is interpreted as the /fractal dimension/ $D_{\rm f}$. Using a")')
write(*,'(" fixed fractal dimension means that the volume filling factor")')
write(*,'(" will decrease with aggregate size. If the parameter is <1, it")')
write(*,'(" is interpreted as a fixed /volume filling factor/ $f$ that")')
write(*,'(" applies to all aggregate sizes - with the implication that then")')
write(*,'(" the fractal dimension increases as a function of size. The")')
write(*,'(" fractal prefactor $k_{\rm f}$ is chosen automatically so that")')
write(*,'(" the asymptotic density of small aggregates is the monomer")')
write(*,'(" material density. To force another value for the prefactor, it")')
write(*,'(" can be given explicitly as the third value of the mmf ")')
write(*,'(" option. The following table summarizes the relevant equations.")')
write(*,'("")')
write(*,'(" | -mmf A0 DF | -mmf A0 FILL ")')
write(*,'(" -------------+-----------------------+-------------- --------------")')
write(*,'(" $f$ | $N^(D_\rm f -3)/3 $ | given by use 3/D_\rm f N^")')
write(*,'(" $D_\rm f $ | given by user | $3\ln N\,/\,\ ")')
write(*,'(" $k_\rm f $ | $(5/3)^D_\rm f /2 $ | $(5/3)^D_\r ")')
write(*,'("")')
write(*,'(" With the structure defined, optool then applies the formalism")')
write(*,'(" from Tazaki & Tanaka (2018) and Tazaki (2021) to compute cross")')
write(*,'(" sections and the scattering matrix. optool also computes the")')
write(*,'(" phase shift $\Delta\phi$ to check the validity of the scattering")')
write(*,'(" matrix. If the condition $\Delta\phi<1$ for accurate scattering")')
write(*,'(" matrix results is violated, a warning will be issued and the")')
write(*,'(" scattering matrix will be set to zero at the relevant")')
write(*,'(" wavelengths. However, the opacities will remain applicable. You")')
write(*,'(" can request to get the result computed under the assumption of")')
write(*,'(" single scattering at wavelengths where the phase shift is too")')
write(*,'(" large. This may be usable for absorbing materials, but we do not")')
write(*,'(" have a clear criterion on when it will be accurate. For this")')
write(*,'(" result, use *-mmfss* instead of *-mmf*. optool will then also")')
write(*,'(" print the wavelength below which the scattering matrix needs to")')
write(*,'(" be used with caution.")')
write(*,'("")')
write(*,'(" CDE")')
write(*,'(" CDE (Continuous Distribution of Ellipsoids) is an analytical")')
write(*,'(" formalism by Bohren & Huffman (1998) to compute the opacity of a")')
write(*,'(" very broad shape distribution. This method is only applicable")')
write(*,'(" in the Rayleigh limit $x=2\pi a\ll\lambda$ and")')
write(*,'(" $|mx|\ll1$. optool will issue a warning if the computation")')
write(*,'(" leaves the bounds of this condition. The scattering matrix will")')
write(*,'(" be computed from a single sphere in the Rayleigh limit.")')
write(*,'("")')
write(*,'("12 Troubleshooting")')
write(*,'("==================")')
write(*,'("")')
write(*,'(" 1. If you get a compilation error about the intrinsic module")')
write(*,'(" ISO_FORTRAN_ENV , compile with make oldio=true .")')
write(*,'(" 2. If you get oscillations in the opacities, in particular at long")')
write(*,'(" wavelengths, the grain size resolution is not sufficient. Use more")')
write(*,'(" grain sizes (*-a*, *-na* and *-d* switches).")')
write(*,'(" 3. If you do not remember how to reproduce a specific run, check the")')
write(*,'(" output file header. It contains the exact command that was used to")')
write(*,'(" produce the file.")')
write(*,'(" 4. If the optool command is not found by your shell, make sure the")')
write(*,'(" optool executable is on your binary search path. Or run it by")')
write(*,'(" giving the full path, like ./optool .")')
write(*,'("")')
write(*,'("13 Acknowledgments")')
write(*,'("==================")')
write(*,'("")')
write(*,'(" We are indebted to")')
write(*,'(" - the [Jena Database of Optical Constants] and the [Aerosol Refractive")')
write(*,'(" Index Archive] for their invaluable collections of refractive index")')
write(*,'(" datasets.")')
write(*,'(" - Rens Waters, Thomas Henning, Xander Tielens, Elisabetta Palumbo,")')
write(*,'(" Laurent Pilon, Jeroen Bouwman, and Melissa McClure for discussions")')
write(*,'(" around optical properties of cosmic dust analogues.")')
write(*,'(" - Charlène Lefèvre for [SIGMA], which inspired me to add grain")')
write(*,'(" mantles.")')
write(*,'(" - Kees Dullemond for discussions about the [RADMC-3D] input format and")')
write(*,'(" the scattering matrix, for the idea to write optool2tex and for")')