diff --git a/.gitignore b/.gitignore
index f27607245..9d509693c 100644
--- a/.gitignore
+++ b/.gitignore
@@ -48,6 +48,7 @@ html
### ignore OS generated files
*~
._*
+.vscode
.DS_Store
.DS_Store?
.Spotlight-V100
diff --git a/HEN_HOUSE/doc/src/pirs509a-beamnrc/figures/dynvmlcd.fig b/HEN_HOUSE/doc/src/pirs509a-beamnrc/figures/dynvmlcd.fig
index 41bf23726..dd309ff4f 100644
--- a/HEN_HOUSE/doc/src/pirs509a-beamnrc/figures/dynvmlcd.fig
+++ b/HEN_HOUSE/doc/src/pirs509a-beamnrc/figures/dynvmlcd.fig
@@ -1,41 +1,37 @@
-#FIG 3.2
+#FIG 3.2 Produced by xfig version 3.2.7a
Portrait
Center
Inches
A4
-100.00
+70.00
Single
-2
1200 2
5 1 0 2 -1 6 0 0 20 0.000 0 1 0 0 3960.674 13702.884 3421 13252 3259 13656 3402 14130
5 1 0 2 -1 6 0 0 20 0.000 0 0 0 0 1230.716 13710.728 1770 13260 1932 13664 1788 14139
6 4770 10260 4950 10485
-4 0 -1 0 0 18 16 0.0000 4 165 135 4809 10462 4\001
+4 0 -1 0 0 18 16 0.0000 4 195 150 4809 10462 4\001
-6
6 5445 10260 5625 10485
-4 0 -1 0 0 18 16 0.0000 4 165 135 5454 10447 5\001
+4 0 -1 0 0 18 16 0.0000 4 195 150 5454 10447 5\001
-6
6 6075 10260 6255 10485
-4 0 -1 0 0 18 16 0.0000 4 165 135 6114 10447 6\001
+4 0 -1 0 0 18 16 0.0000 4 195 150 6114 10447 6\001
-6
6 6795 10260 6975 10485
-4 0 -1 0 0 18 16 0.0000 4 165 135 6834 10447 7\001
+4 0 -1 0 0 18 16 0.0000 4 195 150 6834 10447 7\001
-6
6 7515 10260 7695 10485
-4 0 -1 0 0 18 16 0.0000 4 165 135 7524 10447 8\001
+4 0 -1 0 0 18 16 0.0000 4 195 150 7524 10447 8\001
-6
1 3 0 1 -1 0 0 0 20 0.000 1 0.0000 728 13682 46 46 728 13682 774 13728
1 3 0 1 0 0 50 0 20 0.000 1 0.0000 5091 6016 51 51 5091 6016 5142 6016
-2 3 0 2 -1 6 2 0 20 0.000 0 0 -1 0 0 6
- 6013 13274 6013 14134 7121 14144 7501 13264 6023 13264 6013 13274
2 3 0 2 -1 6 2 0 20 0.000 0 0 -1 0 0 6
9970 13254 8820 13254 9480 14154 9980 14154 9962 13254 9970 13254
2 1 0 2 -1 7 0 0 -1 0.000 0 0 -1 1 0 2
1 1 1.00 60.00 120.00
8006 12094 8006 14554
-2 1 1 2 -1 7 0 0 -1 4.000 0 0 -1 0 0 2
- 7996 12264 7516 13264
-2 1 1 2 -1 7 0 0 -1 4.000 0 0 -1 0 0 2
+2 1 1 2 -1 7 12 0 -1 4.000 0 0 -1 0 0 2
8025 12284 8845 13264
2 2 2 2 -1 7 0 0 -1 3.000 0 0 -1 0 0 5
558 13041 4533 13041 4533 14135 558 14135 558 13041
@@ -55,9 +51,9 @@ Single
2 1 1 2 -1 7 0 0 -1 6.000 0 0 -1 1 0 2
1 1 2.00 60.00 120.00
2547 14276 4509 14276
-2 1 0 2 -1 7 0 0 -1 0.000 0 0 -1 0 0 4
+2 1 0 2 -1 7 4 0 -1 0.000 0 0 -1 0 0 4
4677 13256 4802 13256 4802 13681 4677 13681
-2 2 0 0 -1 7 0 0 20 0.000 0 0 -1 0 0 5
+2 2 0 0 -1 7 1 0 20 0.000 0 0 -1 0 0 5
4589 13318 5777 13318 5777 13593 4589 13593 4589 13318
2 1 0 2 -1 7 0 0 -1 0.000 0 0 -1 1 0 2
1 1 2.00 60.00 120.00
@@ -78,18 +74,12 @@ Single
2 1 0 2 -1 7 0 0 -1 0.000 0 0 -1 1 0 2
1 1 2.00 60.00 120.00
7964 13259 7574 13259
-2 2 0 0 -1 7 0 0 20 0.000 0 0 -1 0 0 5
- 7434 12949 7974 12949 7974 13209 7434 13209 7434 12949
-2 2 0 0 -1 7 0 0 20 0.000 0 0 -1 0 0 5
- 8093 12949 8633 12949 8633 13229 8093 13229 8093 12949
2 1 0 2 -1 7 1 0 -1 0.000 0 0 -1 0 0 4
3421 13252 4516 13252 4520 14126 3407 14131
2 1 1 2 -1 7 0 0 -1 10.000 0 0 -1 0 0 2
565 13693 4527 13693
2 2 0 0 -1 6 50 0 20 0.000 0 0 -1 0 0 5
563 13241 1793 13241 1793 14132 563 14132 563 13241
-2 1 1 2 -1 7 1 0 -1 5.000 0 0 -1 0 0 2
- 6032 13262 10008 13262
2 1 0 2 -1 7 0 0 -1 0.000 0 0 -1 1 1 2
1 1 2.00 60.00 120.00
1 1 2.00 60.00 120.00
@@ -101,10 +91,6 @@ Single
8701 9235 9776 9235 9776 9468 8701 9468 8701 9235
2 2 0 0 -1 7 0 0 20 0.000 0 0 -1 0 0 5
9277 8562 10014 8562 10014 8786 9277 8786 9277 8562
-2 2 2 2 -1 7 0 0 -1 3.000 0 0 -1 0 0 5
- 6017 13056 9992 13056 9992 14150 6017 14150 6017 13056
-2 2 0 0 -1 7 0 0 20 0.000 0 0 -1 0 0 5
- 5998 13120 6735 13120 6735 13395 5998 13395 5998 13120
2 2 0 2 -1 6 0 0 10 0.000 0 0 -1 0 0 5
12726 8773 12791 8773 12791 12018 12726 12018 12726 8773
2 1 0 2 -1 0 0 0 -1 0.000 0 0 -1 0 0 2
@@ -129,26 +115,16 @@ Single
12149 10259 12149 12503
2 2 0 2 -1 6 1 0 10 0.000 0 0 -1 0 0 5
13016 8749 13065 8749 13065 12002 13016 12002 13016 8749
-2 2 0 0 -1 7 0 0 20 0.000 0 0 -1 0 0 5
- 7300 12953 8037 12953 8037 13228 7300 13228 7300 12953
-2 2 0 0 -1 7 0 0 20 0.000 0 0 -1 0 0 5
- 8012 12938 8749 12938 8749 13213 8012 13213 8012 12938
2 1 1 1 0 7 50 0 -1 4.000 0 0 -1 0 0 2
2493 408 2493 5656
2 1 1 1 0 7 50 0 -1 4.000 0 0 -1 0 0 2
2768 408 2768 5687
-2 1 1 1 0 7 50 0 -1 4.000 0 0 -1 0 0 2
- 11848 5435 1060 5442
2 1 1 1 0 7 50 0 -1 4.000 0 0 -1 0 0 2
3256 5199 1060 5199
-2 1 1 1 0 7 50 0 -1 4.000 0 0 -1 0 0 2
- 3256 4589 1060 4589
2 1 1 1 0 7 50 0 -1 4.000 0 0 -1 0 0 2
3317 2575 1120 2575
2 1 1 1 0 7 50 0 -1 4.000 0 0 -1 0 0 2
3347 744 1152 744
-2 1 1 1 0 6 50 0 20 4.000 0 0 -1 0 0 2
- 1486 408 1486 5656
2 1 1 1 0 7 55 0 -1 4.000 0 0 -1 0 0 2
11678 495 1152 500
2 1 1 1 0 7 50 0 -1 4.000 0 0 -1 0 0 2
@@ -159,8 +135,6 @@ Single
1274 774 1274 2575 1486 2575 1486 4589 2493 4589 2493 5199
2280 5199 2280 5442 3013 5442 3013 2819 2768 2819 2768 1049
2067 1049 2067 744 1274 744 1274 806 1274 774
-2 1 1 1 0 7 50 0 -1 4.000 0 0 -1 0 0 2
- 3371 3813 1022 3813
2 1 1 1 0 7 50 0 -1 4.000 0 0 -1 0 0 2
3279 3558 1084 3558
2 1 1 1 0 7 52 0 -1 4.000 0 0 -1 0 0 2
@@ -180,14 +154,6 @@ Single
2 1 0 2 0 7 50 0 -1 0.000 0 0 -1 1 0 2
1 1 1.00 96.00 96.00
3296 2302 3288 2834
-2 1 0 2 0 7 50 0 -1 0.000 0 0 -1 1 1 2
- 1 1 1.00 94.50 94.50
- 1 1 1.00 94.50 94.50
- 1298 876 2036 876
-2 1 1 1 0 7 50 0 -1 4.000 0 0 -1 0 0 2
- 3317 1049 1120 1049
-2 2 0 0 0 7 50 0 20 0.000 0 0 -1 0 0 5
- 1415 945 1999 945 1999 1146 1415 1146 1415 945
2 1 0 2 0 7 50 0 -1 0.000 0 0 -1 1 1 2
1 1 1.00 66.00 36.00
1 1 1.00 66.00 36.00
@@ -211,65 +177,37 @@ Single
1 1 1.00 66.00 36.00
1 1 1.00 66.00 36.00
2297 5137 2491 5137
-2 1 0 2 0 7 50 0 -1 0.000 0 0 -1 1 1 2
- 1 1 1.00 96.00 96.00
- 1 1 1.00 96.00 96.00
- 1500 3975 2992 3975
-2 2 0 0 0 7 50 0 20 0.000 0 0 -1 0 0 5
- 2152 4875 2617 4875 2617 5085 2152 5085 2152 4875
2 1 0 2 0 7 50 0 -1 0.000 0 0 -1 1 0 2
1 1 1.00 96.00 96.00
3141 4162 3143 4585
2 1 0 2 0 7 50 0 -1 0.000 0 0 -1 1 0 2
1 1 1.00 96.00 96.00
3214 4831 3221 5208
-2 2 0 0 0 7 50 0 20 0.000 0 0 -1 0 0 5
- 1762 4320 2227 4320 2227 4530 1762 4530 1762 4320
-2 2 0 0 0 7 50 0 20 0.000 0 0 -1 0 0 5
- 1955 3877 2293 3877 2293 4095 1955 4095 1955 3877
2 2 0 0 -1 7 0 0 20 0.000 0 0 -1 0 0 5
1650 1911 2400 1911 2400 2175 1650 2175 1650 1911
-2 1 1 1 0 7 50 0 -1 4.000 0 0 -1 0 0 2
- 2280 408 2280 5656
2 2 0 2 0 14 51 0 20 0.000 0 0 -1 0 0 5
1496 3551 3011 3551 3011 3813 1496 3813 1496 3551
-2 1 1 1 0 7 51 0 -1 4.000 0 0 -1 0 0 2
- 2067 408 2067 5656
2 1 1 1 0 7 50 0 -1 4.000 0 0 -1 0 0 2
6245 5190 4053 5190
2 1 0 1 0 7 50 0 -1 0.000 0 0 -1 0 0 2
4265 1838 6032 1838
2 1 0 1 0 7 50 0 -1 0.000 0 0 -1 0 0 2
4234 1564 6002 1564
-2 3 0 2 0 6 52 0 20 0.000 0 0 -1 0 0 15
- 5666 5190 5665 2368 6010 2370 6002 1137 4905 1137 4905 802
- 5879 802 5879 498 4234 498 4218 2213 4475 2218 4473 4978
- 5300 4976 5300 5190 5666 5190
-2 1 1 1 0 7 50 0 -1 4.000 0 0 -1 0 0 2
- 6245 802 4053 802
2 1 1 1 0 7 50 0 -1 4.000 0 0 -1 0 0 2
6274 1564 4084 1564
-2 1 1 1 0 7 50 0 -1 4.000 0 0 -1 0 0 2
- 5666 345 5666 5585
2 1 1 1 0 7 50 0 -1 4.000 0 0 -1 0 0 2
5879 345 5879 5616
2 1 1 1 0 7 50 0 -1 4.000 0 0 -1 0 0 2
6002 345 6002 5616
-2 1 1 1 0 7 50 0 -1 4.000 0 0 -1 0 0 2
- 4905 316 4905 5556
-2 1 1 1 0 7 50 0 -1 4.000 0 0 -1 0 0 2
- 5300 316 5300 5556
2 1 1 1 0 7 55 0 -1 4.000 0 0 -1 0 0 2
6245 498 4053 498
-2 1 1 1 0 7 50 0 -1 4.000 0 0 -1 0 0 2
+2 1 1 1 0 7 49 0 -1 4.000 0 0 -1 0 0 2
6245 1137 4053 1137
2 1 1 1 0 7 50 0 -1 4.000 0 0 -1 0 0 2
6274 1838 4084 1838
-2 1 1 1 0 7 52 0 -1 4.000 0 0 -1 0 0 2
- 4232 314 4232 5554
2 2 0 2 0 14 51 0 20 0.000 0 0 -1 0 0 5
4227 1564 6002 1564 6002 1824 4227 1824 4227 1564
-2 1 1 1 0 7 50 0 -1 4.000 0 0 -1 0 0 2
+2 1 1 1 0 7 51 0 -1 4.000 0 0 -1 0 0 2
4482 316 4482 5556
2 1 0 2 0 7 50 0 -1 0.000 0 0 -1 1 0 2
1 1 1.00 96.00 96.00
@@ -289,29 +227,9 @@ Single
2 1 0 2 0 7 50 0 -1 0.000 0 0 -1 1 0 2
1 1 1.00 96.00 96.00
6086 4622 6086 5177
-2 1 0 2 0 7 50 0 -1 0.000 0 0 -1 1 1 2
- 1 1 1.00 94.50 94.50
- 1 1 1.00 94.50 94.50
- 4916 904 5878 904
-2 2 0 0 0 7 50 0 20 0.000 0 0 -1 0 0 5
- 5131 787 5596 787 5596 997 5131 997 5131 787
-2 1 0 2 0 7 50 0 -1 0.000 0 0 -1 1 1 2
- 1 1 1.00 66.00 36.00
- 1 1 1.00 66.00 36.00
- 5307 5294 5656 5294
-2 2 0 0 0 7 50 0 20 0.000 0 0 -1 0 0 5
- 5182 5333 5766 5333 5766 5534 5182 5534 5182 5333
-2 1 0 2 0 7 50 0 -1 0.000 0 0 -1 1 1 2
- 1 1 1.00 94.50 94.50
- 1 1 1.00 94.50 94.50
- 4926 1257 5976 1257
-2 2 0 0 0 7 50 0 20 0.000 0 0 -1 0 0 5
- 5221 1147 5686 1147 5686 1357 5221 1357 5221 1147
2 1 1 1 0 7 50 0 -1 4.000 0 0 -1 0 0 2
8550 393 8550 5672
-2 1 1 1 0 7 50 0 -1 4.000 0 0 -1 0 0 2
- 9040 393 9040 5703
-2 1 1 1 0 7 50 0 -1 4.000 0 0 -1 0 0 2
+2 1 1 1 0 7 51 0 -1 4.000 0 0 -1 0 0 2
9285 393 9285 5703
2 1 1 1 0 7 50 0 -1 4.000 0 0 -1 0 0 2
9530 5457 7325 5457
@@ -329,7 +247,7 @@ Single
7510 3798 9285 3798
2 1 0 1 0 7 50 0 -1 0.000 0 0 -1 0 0 2
7510 4109 9285 4109
-2 3 0 2 0 6 51 0 20 0.000 0 0 -1 0 0 15
+2 3 0 2 0 6 53 0 20 0.000 0 0 -1 0 0 15
7999 731 7999 2574 7627 2573 7627 4688 8550 4691 8550 5211
7754 5211 7754 5457 9285 5457 9285 2817 9040 2817 9040 1038
8337 1038 8337 731 7999 731
@@ -389,10 +307,6 @@ Single
1 1 1.00 94.50 94.50
1 1 1.00 94.50 94.50
7751 5138 8531 5138
-2 1 0 2 0 7 50 0 -1 0.000 0 0 -1 1 1 2
- 1 1 1.00 94.50 94.50
- 1 1 1.00 94.50 94.50
- 7650 4598 8533 4598
2 2 0 0 -1 7 0 0 20 0.000 0 0 -1 0 0 5
8093 3815 8843 3815 8843 4079 8093 4079 8093 3815
2 1 1 1 0 7 50 0 -1 4.000 0 0 -1 0 0 2
@@ -575,124 +489,210 @@ Single
1 1 1.00 60.00 90.00
1 1 1.00 60.00 90.00
6815 14424 7527 14424
-2 2 0 0 0 7 50 0 20 0.000 0 0 -1 0 0 5
- 7902 4895 8367 4895 8367 5105 7902 5105 7902 4895
-2 2 0 0 0 7 50 0 20 0.000 0 0 -1 0 0 5
- 7827 4445 8292 4445 8292 4655 7827 4655 7827 4445
2 1 0 2 0 7 50 0 -1 0.000 0 0 -1 0 0 2
8993 4033 10058 4273
2 1 0 2 0 7 50 0 -1 0.000 0 0 -1 1 1 2
1 1 1.00 96.00 96.00
1 1 1.00 90.00 90.00
11819 521 11827 5426
-2 2 0 0 -1 7 0 0 20 0.000 0 0 -1 0 0 5
+2 2 0 0 -1 7 31 0 20 0.000 0 0 -1 0 0 5
11509 2784 12259 2784 12259 3048 11509 3048 11509 2784
-2 1 0 2 0 7 50 0 -1 0.000 0 0 -1 1 1 2
- 1 1 1.00 96.00 96.00
- 1 1 1.00 96.00 96.00
- 4500 3657 6000 3675
2 1 1 1 0 7 50 0 -1 4.000 0 0 -1 0 0 2
6258 2205 4001 2213
-2 2 0 0 0 7 50 0 20 0.000 0 0 -1 0 0 5
- 4972 3551 5310 3551 5310 3769 4972 3769 4972 3551
2 1 1 1 0 7 50 0 -1 4.000 0 0 -1 0 0 2
7999 393 7999 5672
+2 1 1 1 0 7 50 0 -1 4.000 0 0 -1 0 0 2
+ 8337 393 8337 5672
+2 1 0 2 0 7 50 0 -1 0.000 0 0 -1 1 1 2
+ 1 1 1.00 66.00 36.00
+ 1 1 1.00 66.00 36.00
+ 4918 5266 5655 5265
+2 1 1 1 0 7 50 0 -1 4.000 0 0 -1 0 0 2
+ 4905 310 4905 5550
+2 1 1 1 0 7 50 0 -1 4.000 0 0 -1 0 0 2
+ 5025 318 5025 5558
+2 1 1 1 0 7 50 0 -1 4.000 0 0 -1 0 0 2
+ 5648 295 5648 5535
+2 1 1 1 0 7 50 0 -1 4.000 0 0 -1 0 0 2
+ 3256 4589 1060 4589
+2 1 1 1 0 7 50 0 -1 4.000 0 0 -1 0 0 2
+ 3408 3810 1059 3810
+2 1 1 1 0 7 49 0 -1 4.000 0 0 -1 0 0 2
+ 2067 408 2067 5656
+2 2 0 0 0 7 48 0 20 0.000 0 0 -1 0 0 5
+ 2025 3869 2363 3869 2363 4087 2025 4087 2025 3869
+2 1 1 1 0 7 50 0 -1 4.000 0 0 -1 0 0 2
+ 2280 408 2280 5656
+2 2 0 0 0 7 48 0 20 0.000 0 0 -1 0 0 5
+ 2145 4830 2610 4830 2610 5040 2145 5040 2145 4830
+2 2 0 0 0 7 49 0 20 0.000 0 0 -1 0 0 5
+ 5288 1155 5753 1155 5753 1365 5288 1365 5288 1155
+2 2 0 0 0 7 48 0 20 0.000 0 0 -1 0 0 5
+ 7842 4448 8307 4448 8307 4658 7842 4658 7842 4448
+2 1 0 2 0 7 50 0 -1 0.000 0 0 -1 1 1 2
+ 1 1 1.00 94.50 94.50
+ 1 1 1.00 94.50 94.50
+ 1298 876 2036 876
+2 1 1 1 0 7 51 0 -1 4.000 0 0 -1 0 0 2
+ 3334 1055 1137 1055
+2 2 0 0 0 7 49 0 20 0.000 0 0 -1 0 0 5
+ 1422 946 2006 946 2006 1147 1422 1147 1422 946
+2 1 1 1 0 6 50 0 20 4.000 0 0 -1 0 0 2
+ 1486 408 1486 5656
+2 3 0 2 0 6 52 0 20 0.000 0 0 -1 0 0 15
+ 5666 5190 5665 2368 6010 2370 6002 1137 5025 1133 5025 803
+ 5879 802 5879 498 4234 498 4218 2213 4475 2218 4473 4978
+ 4905 4980 4905 5190 5666 5190
+2 1 1 1 0 7 50 0 -1 4.000 0 0 -1 0 0 2
+ 6190 800 3998 800
+2 1 0 2 0 7 50 0 -1 0.000 0 0 -1 1 1 2
+ 1 1 1.00 94.50 94.50
+ 1 1 1.00 94.50 94.50
+ 5038 867 5883 864
+2 1 1 1 0 7 52 0 -1 4.000 0 0 -1 0 0 2
+ 4232 314 4232 5554
+2 2 0 0 0 7 49 0 20 0.000 0 0 -1 0 0 5
+ 5235 770 5700 770 5700 980 5235 980 5235 770
+2 1 0 2 0 7 50 0 -1 0.000 0 0 -1 1 1 2
+ 1 1 1.00 94.50 94.50
+ 1 1 1.00 94.50 94.50
+ 5047 1215 5998 1219
2 1 0 2 0 7 50 0 -1 0.000 0 0 -1 1 1 2
1 1 1.00 96.00 96.00
1 1 1.00 96.00 96.00
- 7995 1721 9291 1721
-2 2 0 0 0 7 50 0 20 0.000 0 0 -1 0 0 5
- 8332 1593 8670 1593 8670 1811 8332 1811 8332 1593
+ 4485 3772 5995 3770
+2 2 0 0 0 7 49 0 20 0.000 0 0 -1 0 0 5
+ 4932 3635 5270 3635 5270 3853 4932 3853 4932 3635
+2 2 0 0 0 7 49 0 20 0.000 0 0 -1 0 0 5
+ 5006 5304 5590 5304 5590 5505 5006 5505 5006 5304
+2 1 0 2 0 7 50 0 -1 0.000 0 0 -1 1 1 2
+ 1 1 1.00 96.00 96.00
+ 1 1 1.00 96.00 96.00
+ 1524 3950 3016 3950
+2 2 0 0 0 7 48 0 20 0.000 0 0 -1 0 0 5
+ 1770 4340 2235 4340 2235 4550 1770 4550 1770 4340
+2 1 0 2 0 7 50 0 -1 0.000 0 0 -1 1 1 2
+ 1 1 1.00 96.00 96.00
+ 1 1 1.00 96.00 96.00
+ 7995 1715 9291 1715
+2 2 0 0 0 7 48 0 20 0.000 0 0 -1 0 0 5
+ 8352 1607 8690 1607 8690 1825 8352 1825 8352 1607
+2 2 0 0 0 7 48 0 20 0.000 0 0 -1 0 0 5
+ 9062 2022 9400 2022 9400 2240 9062 2240 9062 2022
+2 2 0 0 0 7 48 0 20 0.000 0 0 -1 0 0 5
+ 7507 1382 7845 1382 7845 1600 7507 1600 7507 1382
+2 2 0 0 0 7 48 0 20 0.000 0 0 -1 0 0 5
+ 7913 4883 8378 4883 8378 5093 7913 5093 7913 4883
+2 1 0 2 0 7 50 0 -1 0.000 0 0 -1 1 1 2
+ 1 1 1.00 94.50 94.50
+ 1 1 1.00 94.50 94.50
+ 7650 4598 8533 4598
2 1 1 1 0 7 50 0 -1 4.000 0 0 -1 0 0 2
- 8337 393 8337 5672
-4 0 -1 0 0 18 16 0.0000 4 180 150 8160 14739 Z\001
-4 0 -1 0 0 18 16 0.0000 4 180 150 2687 14859 Z\001
-4 0 -1 0 0 18 16 0.0000 4 225 1230 2851 14528 RMAX_CM\001
-4 0 -1 0 0 18 16 0.0000 4 180 1095 4638 13562 ZTHICK/2\001
-4 0 -1 0 0 18 16 0.0000 4 180 525 1950 13607 NEG\001
-4 0 -1 0 0 18 16 0.0000 4 180 510 2595 13607 POS\001
-4 0 -1 0 0 18 16 0.0000 4 225 1290 8081 12373 ZFOCUS(2)\001
-4 0 -1 0 0 18 16 0.0000 4 180 525 7424 13179 NEG\001
-4 0 -1 0 0 18 16 0.0000 4 180 510 8083 13189 POS\001
-4 0 -1 0 0 14 30 0.0000 4 285 2430 1500 11947 ENDTYPE=0\001
-4 0 -1 0 0 14 30 0.0000 4 300 2430 7171 11882 ENDTYPE=1\001
-4 0 -1 0 0 18 16 0.0000 4 180 1590 511 12377 LEAFRADIUS\001
-4 0 -1 0 0 18 16 0.0000 4 180 165 10095 12384 Y\001
-4 0 -1 0 0 18 16 0.0000 4 180 165 4351 12549 Y\001
-4 0 -1 0 0 19 30 0.0000 4 345 3360 3251 11433 Y cross-section\001
-4 0 -1 0 0 18 16 0.0000 4 225 1230 6560 10775 RMAX_CM\001
-4 0 -1 0 0 18 16 0.0000 4 180 900 8838 9422 ZTHICK\001
-4 0 -1 0 0 18 16 0.0000 4 180 585 9294 8746 ZMIN\001
-4 0 -1 0 0 18 16 0.0000 4 180 585 6060 13353 ZMIN\001
-4 0 -1 0 0 18 16 0.0000 4 180 825 11076 12480 START\001
-4 0 -1 0 0 18 16 0.0000 4 225 1230 12506 12465 RMAX_CM\001
-4 0 -1 0 0 18 16 0.0000 4 180 165 13893 10176 X\001
-4 0 -1 0 0 18 16 0.0000 4 180 165 12075 12771 Y\001
-4 0 0 50 0 18 16 0.0000 4 240 405 1504 61 ztip\001
-4 0 0 50 0 18 16 0.0000 4 180 180 2352 187 zl\001
-4 0 0 50 0 18 16 0.0000 4 180 420 3306 3451 zbh\001
-4 0 0 50 0 18 16 0.0000 4 195 270 3210 2259 zg\001
-4 0 0 50 0 18 16 0.0000 4 240 465 1455 1116 wtip\001
-4 0 0 50 0 18 16 0.0000 4 180 195 1006 1958 zt\001
-4 0 0 50 0 18 16 0.0000 4 180 255 928 3287 wt\001
-4 0 0 50 0 18 16 0.0000 4 195 330 3111 1786 wg\001
-4 0 0 50 0 18 16 0.0000 4 180 345 3078 3145 zth\001
-4 0 0 50 0 18 16 0.0000 4 180 465 2155 5056 wbs\001
-4 0 0 50 0 18 16 0.0000 4 180 330 3083 4125 zts\001
-4 0 0 50 0 18 16 0.0000 4 180 405 3194 4808 zbs\001
-4 0 0 50 0 18 16 0.0000 4 180 390 1796 4493 wts\001
-4 0 0 50 0 18 16 0.0000 4 180 240 2009 4064 wl\001
-4 0 -1 0 0 18 16 0.0000 4 180 675 1695 2130 MED2\001
-4 0 0 50 0 18 16 0.0000 4 180 330 5985 180 zts\001
-4 0 0 50 0 18 16 0.0000 4 180 405 6223 691 zbs\001
-4 0 0 50 0 18 16 0.0000 4 180 345 6162 1327 zth\001
-4 0 0 50 0 18 16 0.0000 4 180 420 6327 1590 zbh\001
-4 0 0 50 0 18 16 0.0000 4 180 180 5714 4380 zl\001
-4 0 0 50 0 18 16 0.0000 4 240 405 6041 4605 ztip\001
-4 0 0 50 0 18 16 0.0000 4 180 390 5187 962 wts\001
-4 0 0 50 0 18 16 0.0000 4 180 465 5226 1325 wbs\001
-4 0 0 50 0 18 16 0.0000 4 240 465 5236 5515 wtip\001
-4 0 0 50 0 18 16 0.0000 4 240 405 8020 111 ztip\001
-4 0 0 50 0 18 16 0.0000 4 180 180 8668 437 zl\001
-4 0 0 50 0 18 16 0.0000 4 180 195 7765 1959 zt\001
-4 0 0 50 0 18 16 0.0000 4 195 270 9062 2197 zg\001
-4 0 0 50 0 18 16 0.0000 4 180 345 9367 3418 zth\001
-4 0 0 50 0 18 16 0.0000 4 180 420 9592 3771 zbh\001
-4 0 0 50 0 18 16 0.0000 4 180 330 9340 4361 zts\001
-4 0 0 50 0 18 16 0.0000 4 180 405 9528 4872 zbs\001
-4 0 0 50 0 18 16 0.0000 4 240 465 7330 1581 wtip\001
-4 0 0 50 0 18 16 0.0000 4 180 255 7074 3367 wt\001
-4 0 0 50 0 18 16 0.0000 4 195 330 9787 2226 wg\001
-4 0 0 50 0 18 16 0.0000 4 180 390 7862 4633 wts\001
-4 0 -1 0 0 18 16 0.0000 4 180 675 8126 4030 MED3\001
-4 0 0 50 0 18 16 0.0000 4 180 195 3925 1869 zt\001
-4 0 0 50 0 18 16 0.0000 4 195 270 6297 2091 zg\001
-4 0 0 50 0 18 16 0.0000 4 180 255 3912 2933 wt\001
-4 0 0 50 0 18 16 0.0000 4 195 330 6087 3017 wg\001
-4 0 -1 0 0 18 16 0.0000 4 180 165 9194 6387 X\001
-4 0 -1 0 0 18 16 0.0000 4 195 600 2335 10457 leaf 1\001
-4 0 -1 0 0 18 16 0.0000 4 165 135 3490 10472 2\001
-4 0 -1 0 0 18 16 0.0000 4 165 135 4195 10475 3\001
-4 0 -1 0 0 18 16 0.0000 4 225 1290 3620 6108 ZFOCUS(1)\001
-4 0 -1 0 0 18 16 0.0000 4 180 825 3695 8433 START\001
-4 0 -1 0 0 18 18 0.0000 4 195 1110 4610 -267 TARGET\001
-4 0 -1 0 0 18 18 0.0000 4 195 690 1880 -252 FULL\001
-4 0 -1 0 0 18 18 0.0000 4 195 1575 7910 -237 ISOCENTER\001
-4 0 -1 0 0 18 16 0.0000 4 180 1170 6530 7638 LEAFGAP\001
-4 0 -1 0 0 18 16 0.0000 4 180 585 10593 572 ZMIN\001
-4 0 -1 0 0 18 16 0.0000 4 180 900 11486 3009 ZTHICK\001
-4 0 -1 0 0 18 16 0.0000 4 180 150 5219 10857 Z\001
-4 0 -1 0 0 18 16 0.0000 4 180 675 7527 9119 MED1\001
-4 0 -1 0 0 18 16 0.0000 4 180 1170 10813 8253 LEAFGAP\001
-4 0 0 50 0 18 16 0.0000 4 180 420 -163 13662 zbh\001
-4 0 0 50 0 18 16 0.0000 4 180 345 103 13388 zth\001
-4 0 -1 0 0 18 16 0.0000 4 180 1185 1094 14588 HOLEPOS\001
-4 0 -1 0 0 18 16 0.0000 4 180 1185 6605 14725 HOLEPOS\001
-4 0 0 50 0 18 16 0.0000 4 180 465 7911 5084 wbs\001
-4 0 -1 0 0 19 30 0.0000 4 435 1950 11306 7685 Top View\001
-4 0 0 50 0 18 16 0.0000 4 240 1530 10162 4266 driving screw\001
-4 0 0 50 0 18 16 0.0000 4 180 495 10162 4506 hole\001
-4 0 -1 0 0 19 40 0.0000 4 435 2970 3660 -1050 DYNVMLC\001
-4 0 -1 0 0 19 30 0.0000 4 345 3615 69 6705 X cross-sections\001
-4 0 0 50 0 18 16 0.0000 4 180 240 5010 3743 wl\001
-4 0 0 50 0 18 16 0.0000 4 180 240 8380 1784 wl\001
+ 11853 5449 1065 5456
+2 1 1 1 0 7 50 0 -1 4.000 0 0 -1 0 0 2
+ 9040 393 9040 5703
+2 1 1 2 -1 7 1 0 -1 5.000 0 0 -1 0 0 2
+ 5552 12603 9528 12603
+2 3 0 2 -1 6 39 0 20 0.000 0 0 -1 0 0 6
+ 6021 13273 6021 14133 7129 14143 7509 13263 6031 13263 6021 13273
+2 2 0 0 -1 7 3 0 20 0.000 0 0 -1 0 0 5
+ 5960 13085 6697 13085 6697 13360 5960 13360 5960 13085
+2 2 0 0 -1 7 4 0 20 0.000 0 0 -1 0 0 5
+ 7226 12883 7963 12883 7963 13158 7226 13158 7226 12883
+2 2 0 0 -1 7 4 0 20 0.000 0 0 -1 0 0 5
+ 8135 12915 8675 12915 8675 13195 8135 13195 8135 12915
+2 1 1 2 -1 7 8 0 -1 4.000 0 0 -1 0 0 2
+ 7996 12264 7516 13264
+2 2 2 2 -1 7 7 0 -1 3.000 0 0 -1 0 0 5
+ 6017 13045 9992 13045 9992 14139 6017 14139 6017 13045
+4 0 -1 0 0 18 16 0.0000 4 210 165 8160 14739 Z\001
+4 0 -1 0 0 18 16 0.0000 4 210 165 2687 14859 Z\001
+4 0 -1 0 0 18 16 0.0000 4 270 1410 2851 14528 RMAX_CM\001
+4 0 -1 0 0 18 16 0.0000 4 210 600 1950 13607 NEG\001
+4 0 -1 0 0 18 16 0.0000 4 210 600 2595 13607 POS\001
+4 0 -1 0 0 18 16 0.0000 4 270 1530 8081 12373 ZFOCUS(2)\001
+4 0 -1 0 0 14 30 0.0000 4 315 2700 1500 11947 ENDTYPE=0\001
+4 0 -1 0 0 14 30 0.0000 4 315 2700 7171 11882 ENDTYPE=1\001
+4 0 -1 0 0 18 16 0.0000 4 210 1770 511 12377 LEAFRADIUS\001
+4 0 -1 0 0 18 16 0.0000 4 210 195 10095 12384 Y\001
+4 0 -1 0 0 18 16 0.0000 4 210 195 4351 12549 Y\001
+4 0 -1 0 0 19 30 0.0000 4 375 3645 3251 11433 Y cross-section\001
+4 0 -1 0 0 18 16 0.0000 4 270 1410 6560 10775 RMAX_CM\001
+4 0 -1 0 0 18 16 0.0000 4 210 1020 8838 9422 ZTHICK\001
+4 0 -1 0 0 18 16 0.0000 4 210 690 9294 8746 ZMIN\001
+4 0 -1 0 0 18 16 0.0000 4 210 900 11076 12480 START\001
+4 0 -1 0 0 18 16 0.0000 4 270 1410 12506 12465 RMAX_CM\001
+4 0 -1 0 0 18 16 0.0000 4 210 180 13893 10176 X\001
+4 0 -1 0 0 18 16 0.0000 4 210 195 12075 12771 Y\001
+4 0 0 50 0 18 16 0.0000 4 270 495 1504 61 ztip\001
+4 0 0 50 0 18 16 0.0000 4 210 225 2352 187 zl\001
+4 0 0 50 0 18 16 0.0000 4 210 495 3306 3451 zbh\001
+4 0 0 50 0 18 16 0.0000 4 210 330 3210 2259 zg\001
+4 0 0 50 0 18 16 0.0000 4 195 240 1006 1958 zt\001
+4 0 0 50 0 18 16 0.0000 4 195 285 928 3287 wt\001
+4 0 0 50 0 18 16 0.0000 4 210 375 3111 1786 wg\001
+4 0 0 50 0 18 16 0.0000 4 210 405 3078 3145 zth\001
+4 0 0 50 0 18 16 0.0000 4 195 390 3083 4125 zts\001
+4 0 0 50 0 18 16 0.0000 4 210 480 3194 4808 zbs\001
+4 0 -1 0 0 18 16 0.0000 4 210 780 1695 2130 MED2\001
+4 0 0 50 0 18 16 0.0000 4 195 390 5985 180 zts\001
+4 0 0 50 0 18 16 0.0000 4 210 480 6223 691 zbs\001
+4 0 0 50 0 18 16 0.0000 4 210 405 6162 1327 zth\001
+4 0 0 50 0 18 16 0.0000 4 210 495 6327 1590 zbh\001
+4 0 0 50 0 18 16 0.0000 4 210 225 5714 4380 zl\001
+4 0 0 50 0 18 16 0.0000 4 270 495 6041 4605 ztip\001
+4 0 0 50 0 18 16 0.0000 4 270 495 8020 111 ztip\001
+4 0 0 50 0 18 16 0.0000 4 210 225 8668 437 zl\001
+4 0 0 50 0 18 16 0.0000 4 195 240 7765 1959 zt\001
+4 0 0 45 0 18 16 0.0000 4 210 330 9062 2197 zg\001
+4 0 0 50 0 18 16 0.0000 4 210 405 9367 3418 zth\001
+4 0 0 50 0 18 16 0.0000 4 210 495 9592 3771 zbh\001
+4 0 0 50 0 18 16 0.0000 4 195 390 9340 4361 zts\001
+4 0 0 50 0 18 16 0.0000 4 210 480 9528 4872 zbs\001
+4 0 0 45 0 18 16 0.0000 4 270 540 7330 1581 wtip\001
+4 0 0 50 0 18 16 0.0000 4 195 285 7074 3367 wt\001
+4 0 0 50 0 18 16 0.0000 4 210 375 9787 2226 wg\001
+4 0 0 50 0 18 16 0.0000 4 195 240 3925 1869 zt\001
+4 0 0 50 0 18 16 0.0000 4 210 330 6297 2091 zg\001
+4 0 0 50 0 18 16 0.0000 4 195 285 3912 2933 wt\001
+4 0 0 50 0 18 16 0.0000 4 210 375 6087 3017 wg\001
+4 0 -1 0 0 18 16 0.0000 4 210 180 9194 6387 X\001
+4 0 -1 0 0 18 16 0.0000 4 210 720 2335 10457 leaf 1\001
+4 0 -1 0 0 18 16 0.0000 4 195 150 3490 10472 2\001
+4 0 -1 0 0 18 16 0.0000 4 195 150 4195 10475 3\001
+4 0 -1 0 0 18 16 0.0000 4 270 1530 3620 6108 ZFOCUS(1)\001
+4 0 -1 0 0 18 16 0.0000 4 210 900 3695 8433 START\001
+4 0 -1 0 0 18 18 0.0000 4 210 1170 4610 -267 TARGET\001
+4 0 -1 0 0 18 18 0.0000 4 210 720 1880 -252 FULL\001
+4 0 -1 0 0 18 18 0.0000 4 210 1680 7910 -237 ISOCENTER\001
+4 0 -1 0 0 18 16 0.0000 4 210 1305 6530 7638 LEAFGAP\001
+4 0 -1 0 0 18 16 0.0000 4 210 690 10593 572 ZMIN\001
+4 0 -1 18 0 18 16 0.0000 4 210 1020 11486 3009 ZTHICK\001
+4 0 -1 0 0 18 16 0.0000 4 210 165 5219 10857 Z\001
+4 0 -1 0 0 18 16 0.0000 4 210 780 7527 9119 MED1\001
+4 0 -1 0 0 18 16 0.0000 4 210 1305 10813 8253 LEAFGAP\001
+4 0 0 50 0 18 16 0.0000 4 210 495 -163 13662 zbh\001
+4 0 0 50 0 18 16 0.0000 4 210 405 103 13388 zth\001
+4 0 -1 0 0 18 16 0.0000 4 210 1380 1094 14588 HOLEPOS\001
+4 0 -1 0 0 18 16 0.0000 4 210 1380 6605 14725 HOLEPOS\001
+4 0 -1 0 0 19 30 0.0000 4 480 2145 11306 7685 Top View\001
+4 0 0 50 0 18 16 0.0000 4 270 1755 10162 4266 driving screw\001
+4 0 0 50 0 18 16 0.0000 4 210 585 10162 4506 hole\001
+4 0 -1 0 0 19 40 0.0000 4 495 3255 3660 -1050 DYNVMLC\001
+4 0 -1 0 0 19 30 0.0000 4 375 3915 69 6705 X cross-sections\001
+4 0 0 45 0 18 16 0.0000 4 195 435 7869 4633 wts\001
+4 0 -1 0 0 18 16 0.0000 4 210 780 8111 4053 MED3\001
+4 0 0 47 0 18 16 0.0000 4 270 540 1476 1118 wtip\001
+4 0 0 47 0 18 16 0.0000 4 195 435 5260 942 wts\001
+4 0 0 47 0 18 16 0.0000 4 210 525 5279 1352 wbs\001
+4 0 0 47 0 18 16 0.0000 4 210 270 4955 3838 wl\001
+4 0 0 47 0 18 16 0.0000 4 270 540 5031 5492 wtip\001
+4 0 0 45 0 18 16 0.0000 4 210 270 2064 4077 wl\001
+4 0 0 46 0 18 16 0.0000 4 195 435 1786 4532 wts\001
+4 0 0 45 0 18 16 0.0000 4 210 525 2089 5019 wbs\001
+4 0 0 45 0 18 16 0.0000 4 210 270 8395 1809 wl\001
+4 0 0 45 0 18 16 0.0000 4 210 525 7911 5073 wbs\001
+4 0 -1 0 0 18 16 0.0000 4 210 690 6022 13323 ZMIN\001
+4 0 -1 0 0 18 16 0.0000 4 210 1245 4623 13578 ZTHICK/2\001
+4 0 -1 0 0 18 16 0.0000 4 210 600 8136 13152 POS\001
+4 0 -1 0 0 18 16 0.0000 4 210 600 7349 13142 NEG\001
diff --git a/HEN_HOUSE/doc/src/pirs509a-beamnrc/figures/dynvmlcd.pdf b/HEN_HOUSE/doc/src/pirs509a-beamnrc/figures/dynvmlcd.pdf
index 59122bedb..1cd088d48 100644
Binary files a/HEN_HOUSE/doc/src/pirs509a-beamnrc/figures/dynvmlcd.pdf and b/HEN_HOUSE/doc/src/pirs509a-beamnrc/figures/dynvmlcd.pdf differ
diff --git a/HEN_HOUSE/doc/src/pirs509a-beamnrc/inputformats/DYNVMLC.inp b/HEN_HOUSE/doc/src/pirs509a-beamnrc/inputformats/DYNVMLC.inp
index ade48d77a..9339c579a 100644
--- a/HEN_HOUSE/doc/src/pirs509a-beamnrc/inputformats/DYNVMLC.inp
+++ b/HEN_HOUSE/doc/src/pirs509a-beamnrc/inputformats/DYNVMLC.inp
@@ -206,7 +206,7 @@
If MODE_$DYNVMLC=1 or 2 (dynamic delivery or step-and-shoot delivery):
- 14b mlc_file (A80)
+ 14b mlc_file (A256)
mlc_file: The full name of the file containing leaf opening
data. The format of the file contents is as follows:
@@ -229,7 +229,7 @@
MLC_TITLE: A title line
NFIELDS_$DYNVMLC: Total number of fields
INDEX_$DYNVMLC(I): Index of field I. 0 <= INDEX_$DYNVMLC(I) <= 1 and
- INDEX_$DYNVMLC(I) >= INDEX_$DYNVMLC(I-1). This
+ INDEX_$DYNVMLC(I) > INDEX_$DYNVMLC(I-1). This
number is compared to a random number on (0,1) at
the start of each history; if the random number is
<= INDEX_$DYNVMLC(I), then field I is used.
@@ -325,7 +325,7 @@
6.7, ZTHICK
0.5, 0.04, 0.04, 0.1354, 0.3252, 0.1227, 48.25, 48.533, 51.524, 51.732,
52.98, 53.28, 2, 54.5474, 54.812, FULL leaf
- 0.25, 0.04, 0.04, 0.0929, 0.1371, 0.1371, 48.345, 48.6096, 49.5277,
+ 0.25, 0.04, 0.04, 0.0929, 0.132, 0.132, 48.345, 48.6096, 49.5277,
49.8277, 2, 51.625, 51.627, 54.7, 54.746, TARGET leaf
0.25, 0.04, 0.04, 0.0354, 0.1285, 0.1235, 48.412, 48.531, 51.631, 51.732,
53.3293, 53.6293, 2, 54.5474, 54.812, ISOCENTER leaf
diff --git a/HEN_HOUSE/doc/src/pirs509a-beamnrc/inputformats/SYNCHDMLC.inp b/HEN_HOUSE/doc/src/pirs509a-beamnrc/inputformats/SYNCHDMLC.inp
index 07e9348d1..5dec45469 100644
--- a/HEN_HOUSE/doc/src/pirs509a-beamnrc/inputformats/SYNCHDMLC.inp
+++ b/HEN_HOUSE/doc/src/pirs509a-beamnrc/inputformats/SYNCHDMLC.inp
@@ -14,9 +14,8 @@
NGROUP_$SYNCHDMLC = number of groups of adjacent leaves where
all leaves in a group are:
1. FULL leaves
- 2. HALF TARGET/ISOCENTER pairs with TARGET leaf
- on the -X (ORIENT=0) or -Y (ORIENT=1) side
- 3. QUARTER TARGET/ISOCENTER pairs
+ 2. HALF TARGET/ISOCENTER leaves
+ 3. QUARTER TARGET/ISOCENTER leaves
NGROUP_$SYNCHDMLC defaults to 3 if set <=0
MODE_$SYNCHCMLC = 0 for single setting of leaf openings (static
field)
@@ -135,7 +134,7 @@
leaf on +X [ORIENT=0] or +Y [ORIENT=1] side of ISOCENTER
leaf only) ZTONGUE_$SYNCHDMLC(1)<=ZGROOVE_$SYNCHDMLC(3)
- Note: If there are no HALF leaf pairs in your model, you can input
+ Note: If there are no HALF leaves in your model, you can input
blanks, zeroes, or arbitrary floating-point numbers for these
cross section dimensions.
@@ -206,15 +205,22 @@
NUM_LEAF_$SYNCHDMLC(I): Number of adjacent leaves in group I
LEAFTYPE: Type of leaf in group I.
Set to: 1 for FULL leaves
- 2 for HALF TARGET/ISOCENTER pair with
- TARGET leaf on the -X (ORIENT=0)
- or -Y (ORIENT=1) side
- 3 for QUARTER TARGET/ISOCENTER pair with
- TARGET leaf on the -X (ORIENT=0)
- or -Y (ORIENT=1) side
-
- Note: If LEAFTYPE is 2 or 3, then you must have an even number
- of leaves in the group.
+ 2 for HALF TARGET/ISOCENTER leaves.
+ Starting leaf on -X (ORIENT=0) or -Y
+ (ORIENT=1) side depends on adjacent leaf:
+ If adjacecent leaf is a HALF or QUARTER
+ TARGET, then it will be HALF ISOCENTER.
+ Otherwise, it will be HALF TARGET.
+ 3 for QUARTER TARGET/ISOCENTER leaves.
+ Starting leaf on -X (ORIENT=0) or -Y
+ (ORIENT=1) side will depend on the
+ adjacent leaf: If it is a HALF or
+ QUARTER TARGET leaf, then the starting
+ leaf will be QUARTER ISOCENTER.
+ Otherwise, it will be QUARTER TARGET.
+
+ Restrictions: 1. Cannot have a HALF or QUARTER TARGET leaf on the -ve
+ side of a FULL leaf.
9 START_$SYNCHDMLC (F15.0) : the start position (cm) wrt the CAX of
leaf 1 as projected to ZMIN_$SYNCHDMLC.
@@ -379,7 +385,7 @@
The collimator starts at Z=47.0 cm, is 7 cm thick and has 80 tungsten
leaves opening in the X direction. Leaves 1-10 and 71-80 are FULL
(type 1), leaves 11-24 and 57-70 are HALF TARGET/HALF ISOCENTER
- pairs, and leaves 25-56 are QUARTER TARGET/QUARTER ISOCENTER pairs.
+ leaves, and leaves 25-56 are QUARTER TARGET/QUARTER ISOCENTER leaves.
The Z focus of the leaf sides is at Z=0 cm which is the position of the
source. The leaf ends are straight and also focused at Z=0 cm.
In this example, the mlc is being run in step-and-shoot mode
@@ -402,9 +408,9 @@
0.11, 0.04, 0.04, 0.035, 0.08, 0.065, 47.2, 47.3, 50.53, 50.52, 52.9, 53.1,
1.7, 53.6, 53.9, QUARTER ISOCENTER cross-section
10, 1, FULL leaves
- 14, 2, HALF TARGET/ISOCENTER pairs
- 32, 3, QUARTER TARGET/ISOCENTER pairs
- 14, 2, HALF TARGET/ISOCENTER pairs
+ 14, 2, HALF TARGET/ISOCENTER leaves
+ 32, 3, QUARTER TARGET/ISOCENTER leaves
+ 14, 2, HALF TARGET/ISOCENTER leaves
10, 1, FULL leaves
-10.5, START
0.005, LEAFGAP
diff --git a/HEN_HOUSE/doc/src/pirs509a-beamnrc/pirs509a-beamnrc.tex b/HEN_HOUSE/doc/src/pirs509a-beamnrc/pirs509a-beamnrc.tex
index 27ec49abd..e2d59e58e 100644
--- a/HEN_HOUSE/doc/src/pirs509a-beamnrc/pirs509a-beamnrc.tex
+++ b/HEN_HOUSE/doc/src/pirs509a-beamnrc/pirs509a-beamnrc.tex
@@ -6842,8 +6842,8 @@ \subsubsubsection{SYNCMLCE}
\subsubsubsection{SYNCHDMLC}
SYNCHDMLC is a version of SYNCVMLC optimized for modeling the High-definition Multileaf
Collimator (HDMLC 120) available on TrueBeam and Novalis linacs. SYNCHDMLC features five
-possible leaf cross sections: FULL, HALF TARGET/HALF ISOCENTER pairs, and QUARTER TARGET/QUARTER ISOCENTER
-pairs. FULL leaves are similar to FULL leaves in DYNVMLC/SYNCVMLC. HALF TARGET/HALF ISOCENTER leaves are
+possible leaf cross sections: FULL, HALF TARGET/HALF ISOCENTER leaves, and QUARTER TARGET/QUARTER ISOCENTER
+leaves. FULL leaves are similar to FULL leaves in DYNVMLC/SYNCVMLC. HALF TARGET/HALF ISOCENTER leaves are
equivalent to the TARGET/ISOCENTER leaves in DYNVMLC/SYNCVMLC. QUARTER TARGET/QUARTER ISOCENTER have the
same shaped cross sections as the HALF TARGET/HALF ISOCENTER leaves, but the cross section perpendicular to
the opening direction is generally thinner. Similar to other synchronized CMs, leaf opening
@@ -8463,8 +8463,9 @@ \subsubsection{SYNCHDMLC}
SYNCHDMLC is a CM coded by Lobo \& Popescu and based on an earlier code by Borges et al\cite{Bo12} for
modeling the high-definition micro mlc (HD120) available on TrueBeam and Novalis linacs. In overall
-design it is similar to DYNVMLC/SYNCVMLC but features two different types of TARGET/ISOCENTER leaf pairs:
-HALF TARGET/HALF ISOCENTER and QUARTER TARGET/QUARTER ISOCENTER. In terms of leaf cross-section
+design it is similar to DYNVMLC/SYNCVMLC but features two different types of TARGET/ISOCENTER leaves:
+HALF TARGET/HALF ISOCENTER and QUARTER TARGET/QUARTER ISOCENTER. In addition, the last TARGET or ISOCENTER
+leaf in a group need not have a matching ISOCENTER or TARGET counterpart. In terms of leaf cross-section
(perpendicular to the opening direction), these are
similar to the TARGET/ISOCENTER leaves in DYNVMLC, with HALF TARGET/HALF ISOCENTER leaves being
wider than the QUARTER TARGET/QUARTER ISOCENTER leaves.
@@ -8490,10 +8491,9 @@ \subsubsection{SYNCHDMLC}
user must input (see Table~\ref{dynvmlc_tab} for label key).
QUARTER leaf cross-sections have similar Z-dimensions
to HALF leaf cross-sections but are thinner in X.
-TARGET/ISOCENTER leaves must always occur in pairs.
In this illustration leaves 1, 2, 11, 12 are FULL; leaves 3, 4, 9, 10
-are HALF TARGET/HALF ISOCENTER pairs; leaves 5--8 are QUARTER TARGET/QUARTER ISOCENTER
-pairs.}
+are HALF TARGET/HALF ISOCENTER leaves; leaves 5--8 are QUARTER TARGET/QUARTER ISOCENTER
+leaves. Unlike DYNVMLC and SYNCVMLC, an odd number of TARGET/ISOCENTER leaves may be specified in a group.}
\label{synchdmlc_fig}
\end{figure}
@@ -8505,7 +8505,7 @@ \subsubsection{SYNCHDMLC}
Table~\ref{dynvmlc_tab}). In general, QUARTER leaves
are used closer to middle of the mlc, and, while their Z-dimensions are similar to those for their HALF counterparts, they
are significantly thinner in the X- ({\tt ORIENT}=0) or Y- ({\tt ORIENT}=1)
-dimension. HALF TARGET/HALF ISOCENTER leaves must occur in pairs, as must QUARTER TARGET/QUARTER ISOCENTER leaves.
+dimension.
\index{SYNCHDMLC!restrictions on leaf dimensions}
Similar to DYNVMLC, leaf cross-section dimensions must be defined so that the Z and X ({\tt ORIENT}=0) or Y ({\tt ORIENT}=1) grid lines
@@ -8523,14 +8523,20 @@ \subsubsection{SYNCHDMLC}
QUARTER ISOCENTER/HALF TARGET ({\tt zg}$_{QUARTER~ISOCENTER}$ $\geq$ {\tt zt}$_{HALF~TARGET}$)\\
and HALF ISOCENTER/FULL (ie {\tt zg}$_{HALF~ISOCENTER}$ $\geq$ {\tt zt}$_{FULL}$).
-As in DYNVMLC, rather than specifying the type of each leaf, the user is able to specify groups
-of adjacent leaves all having the same type. Since HALF TARGET/HALF ISOCENTER leaves and
-QUARTER TARGET/QUARTER ISOCENTER leaves must always occur in pairs, there are really
-only three types to choose from. This also means there must be an even number of
-HALF or QUARTER TARGET/ISOCENTER leaves. Also note that, even though the user
-must specify cross-sections for all leaf types, they need not use all types in a
-given MLC simulation. Thus, there may be no FULL leaves, HALF TARGET/ISOCENTER
-pairs, or QUARTER TARGET/ISOCENTER pairs in the simulation.
+As in DYNVMLC and SYNCVMLC, rather than specifying the type of each leaf in the MLC, the user is able to specify groups
+of adjacent leaves of the same type. Since TARGET leaves always alternate with ISOCENTER leaves within a group,
+for the purpose of defining a leaf group, HALF or QUARTER TARGET/ISOCENTER leaves are considered one type.
+Unlike DYNVMLC and SYNCVMLC, however, the number of TARGET/ISOCENTER leaves in a group does not have to be even. Thus,
+the final TARGET or ISOCENTER leaf in a group need not be matched by its ISOCENTER or TARGET (respectively)
+counterpart. If a group of TARGET/ISOCENTER is defined adjacent to another group of TARGET/ISOCENTER leaves, the
+first leaf in the group will automatically be chosen to match the last leaf of the previous group. Thus, if
+the last leaf of the previous group was a QUARTER or HALF TARGET leaf, then the first leaf of the current
+group will be a QUARTER or HALF ISOCENTER leaf, and vice versa. In terms of leaf groups, the only restriction
+is that the numbers and types of leaves cannot be specified so that a QUARTER or HALF TARGET leaf ends up
+immediately adjacent to a FULL leaf on the -X (ORIENT=0) or -Y (ORIENT=1) side.
+Finally, note that the user need not use all leaf types in a simulated MLC. Thus, there may be no FULL leaves, HALF TARGET/ISOCENTER
+leaves, or QUARTER TARGET/ISOCENTER leaves in the simulation. In this case, inputs for the
+cross-section parameters of the unused leaves can be left as blanks or zeros.
Similar to DYNVMLC and SYNCVMLC, leaf ends can be straight, focused ({\tt ENDTYPE}=1) or cylindrical
({\tt ENDTYPE}=0). For straight, focused ends, the user must specify the Z-position of the
diff --git a/HEN_HOUSE/doc/src/pirs701-egsnrc/inputs/chapter2_3.tex b/HEN_HOUSE/doc/src/pirs701-egsnrc/inputs/chapter2_3.tex
index fecb55f5f..4fd0a7fde 100644
--- a/HEN_HOUSE/doc/src/pirs701-egsnrc/inputs/chapter2_3.tex
+++ b/HEN_HOUSE/doc/src/pirs701-egsnrc/inputs/chapter2_3.tex
@@ -1432,12 +1432,12 @@ \subsubsection{Collision stopping power}
\left[ \ln {T^2 \over I^2} + \ln(1 + \tau/2) + G^\pm(\tau) - \delta \right]
\end{equation}
where $\tau = T/m$ and the functions $G^\pm$ are different for
-electrons and positrons due differences in the M{\o}ller and
+electrons and positrons due to differences in the M{\o}ller and
Bhabha cross sections and are given by
\index{Bhabha scattering}
\begin{equation}
\begin{split}
-G^-(\tau) & = -1 - \beta^2 + \ln\left[ \eta (1 - \eta) \right] +
+G^-(\tau) & = -1 - \beta^2 + \ln\left[4 \eta (1 - \eta) \right] +
{1 \over 1 - \eta} + (1 - \beta^2) \left[ {\tau^2 \eta^2 \over 2} +
(2 \tau + 1) \ln(1 - \eta) \right] \\
G^+(\tau) & = \ln(4 \eta) - \beta^2 \left[ 1 + (2 - y^2) \eta -
diff --git a/HEN_HOUSE/doc/src/pirs701-egsnrc/inputs/new_um_app2.tex b/HEN_HOUSE/doc/src/pirs701-egsnrc/inputs/new_um_app2.tex
index 5d1fc225b..095c721cc 100644
--- a/HEN_HOUSE/doc/src/pirs701-egsnrc/inputs/new_um_app2.tex
+++ b/HEN_HOUSE/doc/src/pirs701-egsnrc/inputs/new_um_app2.tex
@@ -522,7 +522,7 @@ \subsection{ The COMMON Blocks}
&AE &Array(\$MXMED) containing PEGS lower charged
particle cutoff energy for each
medium in MeV.\\
- &UE &Array(\$MXMED) containing PEGS lower charged
+ &UE &Array(\$MXMED) containing PEGS upper charged
particle cutoff energy for each
medium in MeV.\\
&TE &Same as AE except kinetic energy
diff --git a/HEN_HOUSE/doc/src/pirs898-egs++/Doxyfile b/HEN_HOUSE/doc/src/pirs898-egs++/Doxyfile
index 0775afb85..362433333 100644
--- a/HEN_HOUSE/doc/src/pirs898-egs++/Doxyfile
+++ b/HEN_HOUSE/doc/src/pirs898-egs++/Doxyfile
@@ -748,6 +748,7 @@ INPUT = main.doxy \
egs_cbct.doxy \
egs_chamber.doxy \
egs_fac.doxy \
+ egs_kerma.doxy \
../../../egs++ \
../../../interface/egs_interface2.h
diff --git a/HEN_HOUSE/doc/src/pirs898-egs++/egs_kerma.doxy b/HEN_HOUSE/doc/src/pirs898-egs++/egs_kerma.doxy
new file mode 100644
index 000000000..9ed239ff3
--- /dev/null
+++ b/HEN_HOUSE/doc/src/pirs898-egs++/egs_kerma.doxy
@@ -0,0 +1,574 @@
+/*
+###############################################################################
+#
+# EGSnrc egs++ egs_kerma application documentation
+# Copyright (C) 2019 National Research Council Canada
+#
+# This file is part of EGSnrc.
+#
+# EGSnrc is free software: you can redistribute it and/or modify it under
+# the terms of the GNU Affero General Public License as published by the
+# Free Software Foundation, either version 3 of the License, or (at your
+# option) any later version.
+#
+# EGSnrc is distributed in the hope that it will be useful, but WITHOUT ANY
+# WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+# FOR A PARTICULAR PURPOSE. See the GNU Affero General Public License for
+# more details.
+#
+# You should have received a copy of the GNU Affero General Public License
+# along with EGSnrc. If not, see .
+#
+###############################################################################
+#
+# Author: Ernesto Mainegra-Hing, 2019
+#
+# Contributors:
+#
+###############################################################################
+*/
+
+/*! \file egs_kerma.doxy
+ * \brief Documents the egs_kerma application
+ * \EM
+ */
+
+/*! \page egs_kerma egs_kerma: Kerma calculations in a volume
+
+
+- \ref egs_kerma_intro
+- \ref egs_kerma_geometry
+- \ref egs_kerma_sources
+- \ref egs_kerma_options
+
+- \ref egs_kerma_calculation
+ - \ref egs_kerma_FD
+ - \ref egs_kerma_fluence
+
+- \ref egs_kerma_output
+- \ref common_mc
+- \ref egs_kerma_usage
+- \ref egs_kerma_example
+
+
+\section egs_kerma_intro Introduction
+
+The C++ application \c egs_kerma is an advanced EGSnrc application whose
+main goal is the estimation of kerma in a volume defined by one or
+more geometrical regions by means of a track-length estimator.
+The efficiency of the calculation can be increased by several orders
+of magnitude by using the variance reduction technique (VRT) known as
+\ref egs_kerma_FD "forced detection" (FD). It can also be used for estimating the total and
+differential photon fluence in the scoring volume.
+
+This code is by no means in its final stage. It is still very experimental
+and it has been cleaned up for release in the hope it will be useful
+to EGSnrc users. A number of collaborators have been using it for some time and
+have helped with their comments and bug reporting.
+
+\section egs_kerma_geometry Geometries
+
+Geometries are specified in an input file as explained in
+the \ref Geometry "geometry module". Besides the simulation geometry, a geometry
+should be defined that encompasses all scoring regions to be able to take advantage
+of the \ref egs_kerma_FD "FD" VRT.
+
+\section egs_kerma_sources Particle source definition
+
+Any source from the \c egs++ \link Sources particle sources module\endlink can be used. Here is an example for
+an isotropic point source of 40 keV photons placed at 100 cm on the negative z-axis:
+
+\verbatim
+:start source:
+ name = isotropic
+ library = egs_isotropic_source
+ charge = 0
+
+ ### source shape
+ :start shape:
+ type = point
+ position = 0, 0, -100
+ :stop shape:
+
+ ### source spectrum
+ :start spectrum:
+ type = monoenergetic
+ energy = 0.04 # MeV
+ :stop spectrum:
+
+:stop source:
+\endverbatim
+
+BEWARE: No check is made for the charge of the initial particles.
+Although kerma is only defined for neutral particles such as photons and neutrons,
+one could potentially want to calculate kerma for secondary neutral particles
+generated by charged particles.
+
+\section egs_kerma_options Scoring options input block
+
+Kerma in the scoring volume can be calculated
+for several geometries at once. Hence a correlated scoring
+of kerma ratios can be requested by linking any combination of
+two calculation geometries. In cases of small differences between
+these geometries, the kerma ratios will be strongly correlated,
+resulting in a significant reduction of the variance.
+Correlated scoring is requested by assigning two geometry names to
+the correlated geometries
key.
+A pseudo input block of such an input is shown below:
+
+\verbatim
+:start scoring options:
+
+ :start calculation geometry:
+ geometry name = name_1
+ cavity regions = list_of_cavity_region_indices
+ cavity mass = total cavity mass in grams
+ excluded regions = list_of_regions # exclude particles touching these regions
+ :stop calculation geometry:
+ :start calculation geometry:
+ geometry name = name_1
+ cavity regions = list_of_cavity_region_indices
+ cavity mass = total cavity mass in grams
+ excluded regions = list_of_regions
+ :stop calculation geometry:
+ :start calculation geometry:
+ geometry name = name_2
+ ...
+ :stop calculation geometry:
+ ...
+ :start calculation geometry:
+ geometry name = name_n
+ ...
+ :stop calculation geometry:
+
+ correlated geometries = name_i name_j
+ ...
+ correlated geometries = name_k name_l
+
+ ### fluence scoring requested (common to all calculation geometries)
+ :start fluence scoring:
+ minimum energy = Emin
+ maximum energy = Emax
+ number of bins = N
+ scale = linear # linear or logarithmic
+ :stop fluence scoring:
+
+ ### E*muen file (could also be E*mutr): absolute or relative file path
+ emuen file = absolute_file_name
+
+ ### geometry for forced-detection (if omitted, an analog scoring is used)
+ cavity geometry = cavity # For forced detection
+
+:stop scoring options:
+\endverbatim
+
+The meaning of the various keys defining a calculation geometry should
+be self explanatory. The reason one must specify the cavity mass is that
+due to the generality of the geometry package it is not possible to
+automatically compute the cavity volume and therefore the user must supply
+this information for proper normalization. Geometrical regions spanning the
+scoring volume (cavity) are provided via the
+cavity regions
key. If no valid region is entered or
+this input is ignored, no kerma to the cavity of the calculation geometry
+is calculated.
+
+\section egs_kerma_calculation Kerma scoring
+
+Kerma calculations require a file containing pairs of energy \f$ E \f$ and
+\f$ E \times \frac{\mu_{en}}{\rho} \f$ values for the medium of interest.
+This allows a faster evaluation of the
+summation over all contributing particles by saving the multiplication operation.
+The full name (absolute path) of this file must be provided using the
+emuen file
key as shown in the above example.
+
+Kerma for medium \f$m\f$, \f$K_m\f$ is computed by summing up the individual contributions
+of N photons crossing the scoring volume \f$V\f$ using the expression
+
+\f[
+K_m=\sum_{i=1}^{N} w_i \cdot \frac{l_i}{V} \cdot E_i \cdot
+\left(\frac{\mu_\mathrm{en}}{\rho}\right)_i^m
+\f]
+
+The example file emuen_icru90_1.5MeV.data contains \f$ E \times \frac{\mu_{en}}{\rho} \f$
+values for air on a logarithmic energy scale from 1 keV up to 1.5 MeV. These values
+were calculated with the g application using the following (non-default)
+transport parameters:
+
+\verbatim
+##############################
+:start mc transport parameter:
+
+ Global ECUT = 0.512
+ Global PCUT = 0.001
+
+ Photon cross sections = mcdf-xcom # XCOM with renormalized PE xsections
+ Pair cross sections = nrc
+ Triplet production = On
+ Radiative Compton corrections = On
+
+ Brems cross sections = nrc
+ Electron Impact Ionization = penelope # Could be also ik
+
+:stop mc transport parameter:
+##############################
+ \endverbatim
+
+This file is distributed with egs_kerma in the same folder.
+
+\section egs_kerma_FD Variance reduction: Forced Detection
+
+If the user provides the name of a defined geometry encompassing all scoring regions,
+then a variance reduction technique known as forced detection can be used
+which increases the calculation efficiency by about 3 orders of magnitude compared
+to an analog calculation. The name of such a geometry should be provided to the
+application via the key cavity geometry
in the
+scoring options
input block.
+If no such geometry is provided, kerma is calculated in an analog manner,
+\em i.e. by scoring every time a photon crosses the scoring volume.
+
+The usual approach to compute kerma \f$K_m\f$ for medium \f$m\f$ is to score the individual
+contributions from photons crossing the scoring region. This approach is already several
+times more efficient than using the kerma approximation where a full simulation of
+all photon interactions is performed assuming electrons deposit their energy on the spot.
+One can go one step further and estimate the contribution to kerma from any photon
+that could potentially cross the scoring region due to its direction. An algorithm
+to determine whether the direction of the photon intercepts with the scoring region
+ needs to be implemented, as well as a ray-tracing algorithm to account for the attenuation
+through the geometry. This way, one does not have to wait until the photon reaches
+the scoring region to estimate its contribution to the kerma and a more efficient
+sampling of all possible contributions is achieved. This approach is equivalent to
+running N separate MC simulations where in the \f$ k^{th} \f$ simulation, only the contribution
+to the kerma from particles interacting \f$k-1\f$ times is scored. After running all these
+simulations, they can be added together to obtain the total kerma.
+
+Kerma to medium \f$m\f$ is computed by adding up the contribution
+from all photons aimed at the cavity geometry
using the expression
+
+\f[
+K_m=\sum_{i=1}^{N} w_i \cdot \frac{l_i}{V} \cdot E_i \cdot
+\left(\frac{\mu_\mathrm{en}}{\rho}\right)_i^m \cdot
+e^{- \sum_{j=1}^{N_j} \mu_{ij} \cdot l_{ij}}
+\f]
+
+with \f$N_j\f$ the number of regions crossed by the \f$ i^\mathrm{th}\f$ photon,
+\f$\mu_{ij}\f$ the attenuation coefficient in region \f$j\f$ and \f$l_{ij}\f$
+the path across region \f$j\f$.
+
+\section egs_kerma_fluence Fluence scoring
+
+A photon fluence calculation can be requested using the following input block:
+
+\verbatim
+:start scoring options:
+ ...
+ :start fluence scoring:
+ minimum energy = Emin
+ maximum energy = Emax
+ number of bins = N
+ scale = linear, logarithmic
+ :stop fluence scoring:
+:stop scoring options:
+\endverbatim
+
+Differential fluence \f$ \phi_j \f$ is calculated by scoring the volume-averaged
+track length for particles with energies E between \f$E_j\f$ and \f$E_{j+1}\f$
+using Chilton's fluence concept for any sampling volume V:
+
+\f[
+\phi_j=\frac{\sum_{i=1}^{N_j} w_i \cdot l_i}{V}
+\f]
+
+where \f$\displaystyle l_i\f$ is the length of the \f$\displaystyle i^{th} \f$
+particle track. The width of the energy bins is defined by \f$ E_{min} \f$,
+\f$ E_{max} \f$ and the number of bins N. Total fluence \f$ \Phi \f$
+is calculated as the sum over all contributions to the fluence
+
+\f[
+\Phi=\sum_{j=0}^{N-1} \phi_j
+\f]
+
+The fluence scoring
input block can be used
+to request the photon fluence spectrum in a volume \f$V\f$
+defined by the cavity regions
.
+
+\section egs_kerma_output Output
+
+Kerma to medium in the volume of interest for each geometry as well as kerma
+ratios requested by the user are output
+to an *.egslog file (batch execution) or to the screen (interactive execution).
+If a fluence calculation is requested, total and differential photon fluences
+are also output for each calculation geometry. Additionally,
+a Grace plot of
+the differential fluence for all geometries is generated in a file (*.agr).
+
+\section egs_kerma_usage Usage
+
+As any other EGSnrc application, \c egs_kerma can be started from the command line
+using
+\verbatim
+egs_kerma -i input_file [-p pegs_file] [-o output_file] [-b] [-s] [-P N -j i]
+\endverbatim
+where the arguments in square brackets are optional. With the \c -o option
+one can change the name of the output files (by default \c input_file.xxx
+is used, where \c xxx is \c .egslog for the log file, \c .egsdat for the
+data file, etc.). With \c -b one specifies a "batch" run, \em i.e. the
+output is redirected to \c output_file.egslog. With \c -s one can force
+\c egs_kerma to use a simple RCO instead of a JCF-RCO in parallel runs
+specified with -P N -j i
, where \c N is the number of parallel
+jobs and \c i the job index. Note that one Unix-type systems it is easier
+to use the \c exb command to submit parallel jobs
+\verbatim
+exb egs_kerma input_file pegs_file [p=N] [batch=pbs]
+\endverbatim
+where the \c batch option specifies the queuing system to be used.
+If using a \em pegsless input file, then pegs_file
must be substituted with
+the keyword pegsless.
+The EGSnrc GUI can be also used, see
+see PIRS-877 for more details on running EGSnrc applications.
+
+
+\section egs_kerma_example An input example: Kerma in 5 cm air sphere for 40 keV
+photons
+
+\verbatim
+###############################################################################
+#
+# EGSnrc egs++ egs_kerma application sample input file
+# Copyright (C) 2016 National Research Council Canada
+#
+# This file is part of EGSnrc.
+#
+# EGSnrc is free software: you can redistribute it and/or modify it under
+# the terms of the GNU Affero General Public License as published by the
+# Free Software Foundation, either version 3 of the License, or (at your
+# option) any later version.
+#
+# EGSnrc is distributed in the hope that it will be useful, but WITHOUT ANY
+# WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+# FOR A PARTICULAR PURPOSE. See the GNU Affero General Public License for
+# more details.
+#
+# You should have received a copy of the GNU Affero General Public License
+# along with EGSnrc. If not, see .
+#
+###############################################################################
+#
+# Author: Ernesto Mainegra-Hing, 2016
+#
+# Contributors:
+#
+###############################################################################
+#
+# A simple example input file for the egs_kerma C++ application.
+#
+# Simulates a 40 keV isotropic photon source in a room with concrete walls.
+# Kerma and fluence inside a 5 cm radius air sphere are calculated using a
+# forced-detection (FD) technique and a track-length estimator. The air sphere
+# is in the centre of the room and the source is 1 m from the sphere.
+#
+# NOTE 1: Material data is generated in pegsless mode. If source energy
+# changed, make sure to adjust energy cut-offs (ae, ue, ap, up) below.
+#
+# NOTE 2: Two geometries are used here for illustration purposes. Most likely
+# only one geometry is needed for kerma and fluence calculations.
+#
+###############################################################################
+
+
+###############################################################################
+### Geometry
+###############################################################################
+:start geometry definition:
+
+ ### air cavity, spherical, 5 cm radius
+ :start geometry:
+ name = cavity
+ library = egs_spheres
+ midpoint = 0 0 0
+ radii = 5.0
+ :start media input:
+ media = air
+ :stop media input:
+ :stop geometry:
+
+ ### air box (8 m x 8 m x 8 m)
+ :start geometry:
+ name = air
+ library = egs_box
+ box size = 800
+ :start media input:
+ media = air
+ :stop media input:
+ :stop geometry:
+
+ ### room with 1 m thick concrete walls
+ :start geometry:
+ name = walls
+ library = egs_box
+ box size = 900
+ :start media input:
+ media = concrete
+ :stop media input:
+ :stop geometry:
+
+ ### room with concrete walls
+ :start geometry:
+ name = room
+ library = egs_genvelope
+ base geometry = walls
+ inscribed geometries = air
+ :stop geometry:
+
+ ###########################################################################
+ #
+ # The two geometries below are identical
+ #
+ # The purpose is to account for wall contributions to the air sphere. The
+ # first geometry does NOT include contributions from the wall during the
+ # calculation, while the second one does. See the 'scoring options' block
+ # for more detail.
+ #
+ # There are several ways of accomplishing this. One could just have used
+ # the same geometry for both calculations, with and without the sensitive
+ # regions, or have one geometry with and another without the walls.
+ #
+ ###########################################################################
+
+ ### air sphere in room with concrete walls (wall contribution NOT included)
+ :start geometry:
+ name = cavity_in_room_no_wall
+ library = egs_genvelope
+ base geometry = room
+ inscribed geometries = cavity
+ :stop geometry:
+
+ ### air sphere in room with concrete walls (wall contribution included)
+ :start geometry:
+ name = cavity_in_room_all
+ library = egs_genvelope
+ base geometry = room
+ inscribed geometries = cavity
+ :stop geometry:
+
+ ### simulation geometry
+ simulation geometry = cavity_in_room_no_wall
+
+:stop geometry definition:
+
+
+###############################################################################
+### Media
+###############################################################################
+:start media definition:
+
+ ### energy cutoffs
+ ae = 0.512
+ ue = 0.555
+ ap = 0.001
+ up = 0.045
+
+ ### air
+ :start air:
+ density correction file = air_dry_nearsealevel
+ :stop air:
+
+ ### concrete
+ :start concrete:
+ density correction file = concrete_ordinary
+ :stop concrete:
+
+:stop media definition:
+
+
+###############################################################################
+### Source
+###############################################################################
+:start source definition:
+
+ ### isotropic
+ :start source:
+ name = isotropic
+ library = egs_isotropic_source
+ charge = 0
+
+ ### source shape
+ :start shape:
+ type = point
+ position = 0, 0, -100
+ :stop shape:
+
+ ### source spectrum
+ :start spectrum:
+ type = monoenergetic
+ energy = 0.04 # MeV
+ :stop spectrum:
+
+ :stop source:
+
+ ### simulation source
+ simulation source = isotropic
+
+:stop source definition:
+
+
+###############################################################################
+### Scoring options
+###############################################################################
+:start scoring options:
+
+ ### use the same geometry under two different names, for easier bookeeping
+ :start calculation geometry:
+ geometry name = cavity_in_room_no_wall
+ cavity regions = 2
+ excluded regions = 0 # exclude particles passing through these regions
+ cavity mass = 0.630831804841 # 5 cm radius air sphere
+ :stop calculation geometry:
+
+ :start calculation geometry:
+ geometry name = cavity_in_room_all
+ cavity regions = 2
+ #excluded regions = 0 # exclude particles passing through these regions
+ cavity mass = 0.630831804841 # 5 cm radius air sphere
+ :stop calculation geometry:
+
+ ### ratio estimates wall contribution to air sphere
+ correlated geometries = cavity_in_room_all cavity_in_room_no_wall
+
+ ### fluence scoring requested (common to all calculation geometries)
+ :start fluence scoring:
+ minimum energy = 0.0
+ maximum energy = 0.040
+ number of bins = 400
+ scale = linear
+ :stop fluence scoring:
+
+ ### E*muen file (could also be E*mutr): absolute or relative file path
+ emuen file = emuen_icru90_1.5MeV.data
+
+ ### geometry for forced-detection (if omitted, an analog scoring is used)
+ cavity geometry = cavity
+
+:stop scoring options:
+
+
+###############################################################################
+### Transport parameters
+###############################################################################
+:start MC transport parameter:
+
+ ### you can include here any of the EGSnrc transport parameters
+
+ Global ECUT = 2000. # Turn-off electron transport
+ Photon cross sections = mcdf-xcom # XCOM with renormalized PE xsections
+
+:stop MC transport parameter:
+
+
+###############################################################################
+### Run control
+###############################################################################
+:start run control:
+ ncase = 1000000
+:stop run control:
+\endverbatim
+*/
diff --git a/HEN_HOUSE/doc/src/pirs898-egs++/main.doxy b/HEN_HOUSE/doc/src/pirs898-egs++/main.doxy
index 1301f70db..8a13adb1b 100644
--- a/HEN_HOUSE/doc/src/pirs898-egs++/main.doxy
+++ b/HEN_HOUSE/doc/src/pirs898-egs++/main.doxy
@@ -104,6 +104,9 @@
\ref egs_cbct "egs_cbct" An application for CBCT/CT scanner simulations |
+
+ \ref egs_kerma "egs_kerma" An application for volumetric kerma calculations |
+
How to write an egs++ application |
diff --git a/HEN_HOUSE/egs++/Makefile b/HEN_HOUSE/egs++/Makefile
index 316d61fdd..9f29384d9 100644
--- a/HEN_HOUSE/egs++/Makefile
+++ b/HEN_HOUSE/egs++/Makefile
@@ -51,16 +51,18 @@ geometry_libs = egs_planes egs_cd_geometry egs_gtransformed egs_nd_geometry \
egs_box egs_genvelope egs_spheres egs_cylinders egs_iplanes \
egs_cones egs_gstack egs_prism egs_union egs_pyramid egs_conez\
egs_space egs_elliptic_cylinders egs_smart_envelope \
- egs_vhp_geometry egs_octree egs_roundrect_cylinders
+ egs_vhp_geometry egs_octree egs_roundrect_cylinders \
+ egs_rz egs_autoenvelope egs_glib
source_libs = egs_collimated_source egs_isotropic_source egs_parallel_beam \
egs_point_source egs_source_collection egs_transformed_source \
egs_beam_source egs_phsp_source egs_angular_spread \
- iaea_phsp_source egs_radionuclide_source egs_dynamic_source
+ iaea_phsp_source egs_radionuclide_source egs_dynamic_source \
+ egs_fano_source
shape_libs = egs_circle egs_ellipse egs_extended_shape egs_gaussian_shape \
egs_line_shape egs_polygon_shape egs_rectangle egs_shape_collection \
- egs_voxelized_shape
+ egs_voxelized_shape egs_spherical_shell egs_conical_shell
aobject_libs = egs_track_scoring egs_dose_scoring egs_radiative_splitting egs_phsp_scoring
diff --git a/HEN_HOUSE/egs++/ausgab_objects/egs_phsp_scoring/egs_phsp_scoring.cpp b/HEN_HOUSE/egs++/ausgab_objects/egs_phsp_scoring/egs_phsp_scoring.cpp
index 79059f4ee..958a8b71c 100644
--- a/HEN_HOUSE/egs++/ausgab_objects/egs_phsp_scoring/egs_phsp_scoring.cpp
+++ b/HEN_HOUSE/egs++/ausgab_objects/egs_phsp_scoring/egs_phsp_scoring.cpp
@@ -322,7 +322,6 @@ void EGS_PhspScoring::openPhspFile() const {
}
}
else if (oformat == 1) { //IAEA format
- int rwmode;
int iaea_iostat;
iaea_id = 1; //numerical index indicating this is the 1st file associated with this object scored
//hard coded to 1
@@ -525,7 +524,7 @@ extern "C" {
int stype = 0; //default is to use scoring geom
int phspouttype;
int ptype;
- int sdir;
+ int sdir=0;
int imuscore = 0;
float xyzconst[3];
bool xyzisconst[3] = {false, false, false};
diff --git a/HEN_HOUSE/egs++/ausgab_objects/egs_phsp_scoring/egs_phsp_scoring.h b/HEN_HOUSE/egs++/ausgab_objects/egs_phsp_scoring/egs_phsp_scoring.h
index c6966f88a..9fd6caddb 100644
--- a/HEN_HOUSE/egs++/ausgab_objects/egs_phsp_scoring/egs_phsp_scoring.h
+++ b/HEN_HOUSE/egs++/ausgab_objects/egs_phsp_scoring/egs_phsp_scoring.h
@@ -484,11 +484,7 @@ class EGS_PHSP_SCORING_EXPORT EGS_PhspScoring : public EGS_AusgabObject {
wt = p.w >= 0 ? p.wt : -p.wt;
};
};
- mutable fstream phsp_file; //output file -- mutable so we can write to it during storeState
- EGS_I64 count; //total no. of particles in file
- EGS_I64 countg; //no. of photons in file
- float emax; //max. k.e. of particles in phsp file
- float emin; //min. k.e. of charged particles in file
+
bool score_mc; //set to true to score multiple crossers and their descendents
//variables specific to IAEA format
@@ -514,8 +510,13 @@ class EGS_PHSP_SCORING_EXPORT EGS_PhspScoring : public EGS_AusgabObject {
return (1 << 31);
}
- int store_max; //max. no. of particles to store in p_stack
mutable int phsp_index; //index in p_stack array -- mutable so we can change it in storeState
+ int store_max; //max. no. of particles to store in p_stack
+ mutable fstream phsp_file; //output file -- mutable so we can write to it during storeState
+ EGS_I64 count; //total no. of particles in file
+ EGS_I64 countg; //no. of photons in file
+ float emin; //min. k.e. of charged particles in file
+ float emax; //max. k.e. of particles in phsp file
mutable bool first_flush; //first time writing to file in this run -- mutable so we can change it in storeState
bool is_restart; //true if this is a restart
diff --git a/HEN_HOUSE/egs++/egs_alias_table.cpp b/HEN_HOUSE/egs++/egs_alias_table.cpp
index 207887d44..689db96c6 100644
--- a/HEN_HOUSE/egs++/egs_alias_table.cpp
+++ b/HEN_HOUSE/egs++/egs_alias_table.cpp
@@ -128,12 +128,10 @@ void EGS_AliasTable::make() {
else {
fcum[i] = 0.5*(fi[i]+fi[i+1])*(xi[i+1]-xi[i]);
}
- //egsWarning("i=%d fcum=%g x=%g f=%g\n",i,fcum[i],xi[i],fi[i]);
sum += fcum[i];
wi[i] = 1;
bin[i] = 0;
not_done[i] = true;
- wi[i] = 1;
if (type == 0) {
sum1 += fcum[i]*xi[i];
}
@@ -144,6 +142,7 @@ void EGS_AliasTable::make() {
fi[i+1]*(xi[i]+2*xi[i+1]))/(3*(fi[i]+fi[i+1]));
}
average = sum1/sum;
+
for (i=0; i sum) {
+ high_bin = jh;
break;
}
}
- for (jl=0; jl n && fcum[jh] > sum) {
- break;
- }
+
+ // normalize distribution
+ if (sum <= 0) {
+ egsFatal("Error: %s, line %d: degenerate distribution, histogram sum <= 0", __FILE__, __LINE__);
+ }
+ else {
+ for (i=0; i big_list; // bins above average
+ vector small_list; // bins below average
+ for (i=0; i n && fcum[jl] < sum) {
- break;
- }
- }
- jl_old = jl;
+ big_list.push_back(i);
}
- double aux = sum - fcum[jl];
- fcum[jh] -= aux;
- wi[jl] = fcum[jl]/sum;
- bins[jl] = jh;
- //egsInformation("i=%d jh=%d jl=%d w=%g\n",i,jh,jl,wi[jl]);
- jh_old = jh; //jl_old = jl;
}
- for (i=0; i n) {
- bins[i] = i;
- wi[i] = 1;
+
+ // alias
+ int loopCount=0;
+ while (big_list.size() > 0 && small_list.size() > 0 && loopCount++ <= loopMax) {
+
+ // get a pair of big and small bins
+ int big = big_list.back();
+ int small = small_list.back();
+
+ // alias: fill small bin
+ bins[small] = big; // alias small to big
+ wi[big] -= (1.0 - wi[small]); // remove aliased portion from big bin
+ small_list.pop_back(); // small bin is now filled
+
+ // check if big bin is now small
+ if (wi[big] < 1.0 + epsilon) {
+ big_list.pop_back();
+ small_list.push_back(big);
}
- //egsInformation("alias table: %d %d %g\n",i,bins[i],wi[i]);
}
- delete [] fcum;
+ if (big_list.size() > 0) {
+ egsWarning("Warning: %s, line %d: table aliasing may be incomplete", __FILE__, __LINE__);
+ }
+
+ // delete local arrays
+ delete [] p;
}
+
EGS_SimpleAliasTable::~EGS_SimpleAliasTable() {
if (n > 0) {
delete [] wi;
delete [] bins;
}
}
-
diff --git a/HEN_HOUSE/egs++/egs_alias_table.h b/HEN_HOUSE/egs++/egs_alias_table.h
index e7bf2aade..922fbcbe1 100644
--- a/HEN_HOUSE/egs++/egs_alias_table.h
+++ b/HEN_HOUSE/egs++/egs_alias_table.h
@@ -39,6 +39,7 @@
#include "egs_libconfig.h"
#include "egs_rndm.h"
+#include
typedef EGS_Float(*EGS_AtFunction)(EGS_Float,void *);
diff --git a/HEN_HOUSE/egs++/egs_application.cpp b/HEN_HOUSE/egs++/egs_application.cpp
index b7c06edf4..c4159c3de 100644
--- a/HEN_HOUSE/egs++/egs_application.cpp
+++ b/HEN_HOUSE/egs++/egs_application.cpp
@@ -245,6 +245,12 @@ EGS_Application::EGS_Application(int argc, char **argv) : input(0), geometry(0),
//
if (!getArgument(argc,argv,"-a","--application",app_name)) {
app_name = egsStripPath(argv[0]);
+
+ // In Windows PowerShell, we need to remove a .exe extension
+ size_t exePos = app_name.rfind(".exe");
+ if (exePos != string::npos) {
+ app_name = app_name.substr(0, exePos);
+ }
}
if (!app_name.size()) egsFatal("%s\n failed to determine application "
"name from %s or command line arguments\n",__egs_app_msg1,argv[0]);
diff --git a/HEN_HOUSE/egs++/egs_polygon.h b/HEN_HOUSE/egs++/egs_polygon.h
index 4d85a8a1d..6518b7e9a 100644
--- a/HEN_HOUSE/egs++/egs_polygon.h
+++ b/HEN_HOUSE/egs++/egs_polygon.h
@@ -251,7 +251,12 @@ class EGS_EXPORT EGS_2DPolygon {
if (tt <= t+epsilon) {
EGS_Float lam = uj[j]*(x-p[j]+u*tt);
if (lam >= 0 && lam < uj[j].length2()) {
- t = tt;
+ if (tt < 0) {
+ t = 0;
+ }
+ else {
+ t = tt;
+ }
res = true;
jhit = j;
//if( normal ) *normal = a[j]*(-1);
@@ -437,7 +442,7 @@ class EGS_EXPORT EGS_PolygonT {
/*! \brief Will the line defined by position \a x and direction \a u
intersect the polygon plane at a position inside the polygon ?
- The interprtation of the return value and the value of \a in and
+ The interpretation of the return value and the value of \a in and
\a t is the same as in EGS_2DPolygon::howfar()
*/
inline bool howfar(bool in, const EGS_Vector &x, const EGS_Vector &u,
diff --git a/HEN_HOUSE/egs++/egs_transformations.h b/HEN_HOUSE/egs++/egs_transformations.h
index f48537c22..dac540615 100644
--- a/HEN_HOUSE/egs++/egs_transformations.h
+++ b/HEN_HOUSE/egs++/egs_transformations.h
@@ -215,11 +215,23 @@ class EGS_EXPORT EGS_RotationMatrix {
(rzx == m.rxz) && (rzy == m.rzy) && (rzz == m.rzz)) ? true : false;
};
- /*! \brief Returns \c true, if this object is a unity matrix, \a false
+ /*! \brief Returns \c true, if this object is approximately the unit matrix, \a false
otherwise */
bool isI() const {
- return (rxx == 1 && rxy == 0 && rxz == 0 && ryx == 0 && ryy == 1 &&
- ryz == 0 && rzx == 0 && rzy == 0 && rzz == 1);
+
+ // compare with unit matrix components
+ if (fabs(rxy) < epsilon &&
+ fabs(rxz) < epsilon &&
+ fabs(ryx) < epsilon &&
+ fabs(ryz) < epsilon &&
+ fabs(rzx) < epsilon &&
+ fabs(rzy) < epsilon &&
+ fabs(rxx-1) < epsilon &&
+ fabs(ryy-1) < epsilon &&
+ fabs(rzz-1) < epsilon) {
+ return true;
+ }
+ return false;
};
/*! \brief Returns the rotated a vector \f$ R \cdot \vec{v}\f$ */
diff --git a/HEN_HOUSE/egs++/geometry/egs_autoenvelope/Makefile b/HEN_HOUSE/egs++/geometry/egs_autoenvelope/Makefile
new file mode 100644
index 000000000..4e5c51a88
--- /dev/null
+++ b/HEN_HOUSE/egs++/geometry/egs_autoenvelope/Makefile
@@ -0,0 +1,90 @@
+
+###############################################################################
+#
+# EGSnrc egs++ makefile to build auto envelope geometry
+# Copyright (C) 2016 Randle E. P. Taylor, Rowan M. Thomson,
+# Marc J. P. Chamberland, D. W. O. Rogers
+#
+# This file is part of EGSnrc.
+#
+# EGSnrc is free software: you can redistribute it and/or modify it under
+# the terms of the GNU Affero General Public License as published by the
+# Free Software Foundation, either version 3 of the License, or (at your
+# option) any later version.
+#
+# EGSnrc is distributed in the hope that it will be useful, but WITHOUT ANY
+# WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+# FOR A PARTICULAR PURPOSE. See the GNU Affero General Public License for
+# more details.
+#
+# You should have received a copy of the GNU Affero General Public License
+# along with EGSnrc. If not, see .
+#
+###############################################################################
+#
+# Author: Randle Taylor, 2016
+#
+# Contributors: Marc Chamberland
+# Rowan Thomson
+# Dave Rogers
+#
+###############################################################################
+
+
+include $(EGS_CONFIG)
+include $(SPEC_DIR)egspp.spec
+include $(SPEC_DIR)egspp_$(my_machine).conf
+
+SOBOL =
+SOBOL_DEF =
+SOBOL_H =
+ifneq ("$(wildcard egs_sobol.cpp)","")
+SOBOL = sobol egs_sobol
+SOBOL_DEF = -DBUILD_SOBOL_DLL -DBUILD_AE_SOBOL -DHAS_SOBOL
+SOBOL_H = sobol.h egs_sobol.h
+endif
+
+GZSTREAM =
+GZSTREAM_DEF =
+GZSTREAM_H =
+GZSTREAM_LIB =
+ifneq ("$(wildcard gzstream.cpp)","")
+GZSTREAM = gzstream
+GZSTREAM_DEF = -DBUILD_GZSTREAM -DHAS_GZSTREAM
+GZSTREAM_H = gzstream.h
+GZSTREAM_LIB = z
+endif
+
+DEFS = $(DEF1) -DBUILD_AENVELOPE_DLL $(SOBOL_DEF) $(GZSTREAM_DEF)
+
+
+library = egs_autoenvelope
+lib_files = egs_autoenvelope $(SOBOL) $(GZSTREAM)
+
+link2_libs = egs_gtransformed egspp $(GZSTREAM_LIB)
+
+
+include $(SPEC_DIR)egspp_libs.spec
+
+$(make_depend)
+
+
+$(DSO2)autoenvelope.$(obje): volcor.h $(GZSTREAM_H) \
+ ..$(DSEP)egs_gtransformed$(DSEP)egs_gtransformed.h
+
+$(ABS_DSO)$(libpre)egs_autoenvelope$(libext): volcor.h $(GZSTREAM_H) \
+ ..$(DSEP)egs_gtransformed$(DSEP)egs_gtransformed.h
+
+clean:
+ $(REMOVE) $(ABS_DSO)$(libpre)egs_autoenvelope$(libext) $(ABS_DSO)$(libpre)egs_sobol$(libext) $(ABS_DSO)$(libpre)sobol$(libext) \
+ $(ABS_DSO)$(libpre)gzstream$(libext) $(ABS_DSO)$(libpre)gzstream$(libext)
+ $(REMOVE) $(DSO2)autoenvelope.$(obje) $(DSO2)sobol.$(obje) $(DSO2)egs_sobol.$(obje) \
+ $(DSO2)gzstream.$(obje) $(DSO2)gzstream.$(obje) \
+
+touchup:
+ touch egs_autoenvelope.cpp
+
+fresh: clean touchup $(lib_objects)
+
+
+.PHONY: clean touchup fresh
diff --git a/HEN_HOUSE/egs++/geometry/egs_autoenvelope/egs_autoenvelope.cpp b/HEN_HOUSE/egs++/geometry/egs_autoenvelope/egs_autoenvelope.cpp
new file mode 100644
index 000000000..5f16882c2
--- /dev/null
+++ b/HEN_HOUSE/egs++/geometry/egs_autoenvelope/egs_autoenvelope.cpp
@@ -0,0 +1,1061 @@
+/*
+###############################################################################
+#
+# EGSnrc egs++ auto envelope geometry
+# Copyright (C) 2016 Randle E. P. Taylor, Rowan M. Thomson,
+# Marc J. P. Chamberland, D. W. O. Rogers
+#
+# This file is part of EGSnrc.
+#
+# EGSnrc is free software: you can redistribute it and/or modify it under
+# the terms of the GNU Affero General Public License as published by the
+# Free Software Foundation, either version 3 of the License, or (at your
+# option) any later version.
+#
+# EGSnrc is distributed in the hope that it will be useful, but WITHOUT ANY
+# WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+# FOR A PARTICULAR PURPOSE. See the GNU Affero General Public License for
+# more details.
+#
+# You should have received a copy of the GNU Affero General Public License
+# along with EGSnrc. If not, see .
+#
+###############################################################################
+#
+# Author: Randle Taylor, 2016
+#
+# Contributors: Marc Chamberland
+# Rowan Thomson
+# Dave Rogers
+#
+###############################################################################
+#
+# egs_autoenvelope was developed for the Carleton Laboratory for
+# Radiotherapy Physics.
+#
+###############################################################################
+*/
+
+
+/*! \file egs_autoenvelope.cpp
+ * \brief A fast envelope geometry (based on EGS_FastEnvelope) with automatic region detection
+ * \author Randle Taylor (randle.taylor@gmail.com)
+ */
+
+#include "egs_input.h"
+#include "egs_autoenvelope.h"
+#include "../egs_gtransformed/egs_gtransformed.h"
+
+#include
+#include
+#include
+#include
+#include
+
+#ifdef HAS_SOBOL
+ #include "sobol.h"
+#endif
+
+
+string EGS_AENVELOPE_LOCAL EGS_AEnvelope::type = "EGS_AEnvelope";
+string EGS_AENVELOPE_LOCAL EGS_ASwitchedEnvelope::type = "EGS_ASwitchedEnvelope";
+/* only geometries that support getting regional volume/mass can be used as phantoms */
+const string EGS_AENVELOPE_LOCAL EGS_AEnvelope::allowed_base_geom_types[] = {"EGS_cSpheres", "EGS_cSphericalShell", "EGS_XYZGeometry", "EGS_RZ"};
+
+static char EGS_AENVELOPE_LOCAL geom_class_msg[] = "createGeometry(AEnvelope): %s\n";
+static char EGS_AENVELOPE_LOCAL base_geom_keyword[] = "base geometry";
+static char EGS_AENVELOPE_LOCAL inscribed_geom_keyword[] = "inscribed geometry";
+static char EGS_AENVELOPE_LOCAL inscribed_geom_name_keyword[] = "inscribed geometry name";
+static char EGS_AENVELOPE_LOCAL transformations_keyword[] = "transformations";
+static char EGS_AENVELOPE_LOCAL type_keyword[] = "type";
+static char EGS_AENVELOPE_LOCAL transformation_keyword[] = "transformation";
+
+
+EGS_AEnvelope::EGS_AEnvelope(EGS_BaseGeometry *base_geom,
+ const vector inscribed, const string &Name, bool debug, string output_vc_file) :
+ EGS_BaseGeometry(Name), base_geom(base_geom), debug_info(debug), output_vc(output_vc_file) {
+
+ bool volcor_available = allowedBaseGeomType(base_geom->getType());
+ bool volcor_requested = inscribed.size() > 0 && inscribed[0].vcopts->mode == CORRECT_VOLUME;
+
+ if (!volcor_available && volcor_requested) {
+
+ string msg(
+ "EGS_AEnvelope:: Volume correction is not available for geometry type '%s (%s)'. "
+ "Geometry types must implement getMass. Valid choices are:\n\t"
+ );
+
+ int end = (int)(sizeof(allowed_base_geom_types)/sizeof(string));
+ for (int i=0; i < end; i++) {
+ msg += allowed_base_geom_types[i] + " ";
+ }
+ msg += "\nPlease set `action = discover -or- correct and zero volume` or use a different base geometry type.\n";
+ egsFatal(msg.c_str(), base_geom->getType().c_str(), base_geom->getName().c_str());
+ }
+
+
+ base_geom->ref();
+ nregbase = base_geom->regions();
+ is_convex = base_geom->isConvex();
+ has_rho_scaling = getHasRhoScaling();
+
+ // initialize uncorrected/corrected masses
+ for (int ir = 0; ir < nregbase; ir++) {
+ uncorrected_mass.push_back(base_geom->getMass(ir));
+ corrected_mass.push_back(base_geom->getMass(ir));
+ }
+
+ ninscribed = inscribed.size();
+ if (ninscribed == 0) {
+ egsFatal("EGS_AEnvelope: no inscribed geometries!\n");
+ }
+
+ // nreg is total regions in geometry = nregbase + sum(nreginscribed_j)
+ nreg = nregbase;
+
+ // create a transformed copy of the inscribed geometry
+ int global = nreg;
+ for (int geom_idx=0; geom_idx < ninscribed; geom_idx++) {
+
+ EGS_BaseGeometry *inscribed_geom = inscribed[geom_idx].geom;
+ EGS_AffineTransform *transform = inscribed[geom_idx].transform;
+ VCOptions *vcopt = inscribed[geom_idx].vcopts;
+ EGS_TransformedGeometry *transformed = new EGS_TransformedGeometry(inscribed_geom, *transform);
+
+ inscribed_geoms.push_back(transformed);
+ transforms.push_back(transform);
+ opts.push_back(vcopt);
+
+ int ninscribed_reg= transformed->regions();
+ nreg += ninscribed_reg;
+
+ for (int lreg = 0; lreg < ninscribed_reg; lreg++) {
+ GeomRegPairT local(transformed, lreg);
+ local_to_global_reg[local] = global;
+ global_reg_to_local[global] = local;
+ global++;
+ }
+ }
+
+ geoms_in_region = new vector[nregbase];
+
+ if (inscribed[0].vcopts->vc_file != "") {
+ vc_results = loadFileResults(inscribed[0].vcopts, base_geom, inscribed_geoms, transforms);
+ egsInformation("loaded from %s\n", inscribed[0].vcopts->vc_file.c_str());
+ }
+ else {
+ // now run volume correction and figure out which regions have inscribed geometries
+ vc_results = findRegionsWithInscribed(inscribed[0].vcopts, base_geom, inscribed_geoms, transforms);
+ }
+
+
+ nreg_with_inscribed = 0;
+ for (volcor::RegionGeomSetT::iterator it=vc_results.regions_with_inscribed.begin();
+ it!=vc_results.regions_with_inscribed.end(); it++) {
+ nreg_with_inscribed++;
+ copy(it->second.begin(), it->second.end(), back_inserter(geoms_in_region[it->first]));
+ }
+
+ // set mass correction based on volume correction
+ for (int ir = 0; ir < nregbase; ir++) {
+ double volume_ratio = vc_results.corrected_volumes[ir]/vc_results.uncorrected_volumes[ir];
+ corrected_mass[ir] = uncorrected_mass[ir]*volume_ratio;
+ }
+
+ if (getNRegWithInscribed() == 0) {
+ egsFatal("EGS_AEnvelope:: Failed to find any regions with inscribed geometries\n");
+ }
+
+ if (output_vc=="yes" || output_vc=="text" || output_vc=="gzip") {
+ writeVolumeCorrection();
+ }
+
+ if (debug_info) {
+ printInfo();
+ }
+
+}
+
+
+EGS_AEnvelope::~EGS_AEnvelope() {
+ if (!base_geom->deref()) {
+ delete base_geom;
+ }
+}
+
+bool EGS_AEnvelope::allowedBaseGeomType(const string &geom_type) {
+ // Check if the input geometry is one that aenvelope can handle
+
+ int end = (int)(sizeof(allowed_base_geom_types)/sizeof(string));
+
+ for (int i=0; i::const_iterator glob_it = local_to_global_reg.find(local);
+ if (glob_it != local_to_global_reg.end()) {
+ return (*glob_it).second;
+ }
+ return -1;
+}
+
+volcor::GeomRegPairT EGS_AEnvelope::getLocalFromGlobalReg(int ireg) const {
+
+ map::const_iterator glob_it = global_reg_to_local.find(ireg);
+ if (glob_it != global_reg_to_local.end()) {
+ return (*glob_it).second;
+ }
+ return volcor::GeomRegPairT();
+}
+
+int EGS_AEnvelope::getNRegWithInscribed() const {
+ return nreg_with_inscribed;
+}
+
+bool EGS_AEnvelope::isRealRegion(int ireg) const {
+
+ bool is_outside = ireg < 0 || ireg >= nreg;
+ if (is_outside) {
+ return false;
+ }
+
+ bool in_base_geom = ireg < nregbase;
+ if (in_base_geom) {
+ return base_geom->isRealRegion(ireg);
+ }
+
+ volcor::GeomRegPairT local;
+ local = getLocalFromGlobalReg(ireg);
+ return local.first->isRealRegion(local.second);
+
+};
+
+
+bool EGS_AEnvelope::isInside(const EGS_Vector &x) {
+ return base_geom->isInside(x);
+};
+
+
+int EGS_AEnvelope::isWhere(const EGS_Vector &x) {
+
+ int base_reg = base_geom->isWhere(x);
+
+ // which inscribed geometries are in this region
+ vector geoms_in_reg = getGeomsInRegion(base_reg);
+
+ if (base_reg < 0 || geoms_in_reg.size()==0) {
+ return base_reg;
+ }
+
+ // loop over all inscribed geometries in this region and check
+ // if position is inside any of them
+ for (vector::iterator geom=geoms_in_reg.begin(); geom!= geoms_in_reg.end(); ++geom) {
+ int inscribed_reg = (*geom)->isWhere(x);
+ if (inscribed_reg >= 0) {
+ return getGlobalRegFromLocalReg(*geom, inscribed_reg);
+ }
+ }
+
+ // not in any of the inscribed geometries so must be in envelope region
+ return base_reg;
+};
+
+int EGS_AEnvelope::inside(const EGS_Vector &x) {
+ return isWhere(x);
+};
+
+int EGS_AEnvelope::medium(int ireg) const {
+
+ if (ireg < nregbase) {
+ return base_geom->medium(ireg);
+ }
+
+ volcor::GeomRegPairT local = getLocalFromGlobalReg(ireg);
+
+ return local.first->medium(local.second);
+};
+
+vector EGS_AEnvelope::getGeomsInRegion(int ireg) {
+ if (ireg < 0 || ireg >= nregbase) {
+ vector empty;
+ return empty;
+ }
+ return geoms_in_region[ireg];
+};
+
+int EGS_AEnvelope::computeIntersections(int ireg, int n, const EGS_Vector &X,
+ const EGS_Vector &u, EGS_GeometryIntersections *isections) {
+
+ // For a given position x, direction u and region number ireg, this
+ // method returns the number of intersections with the geometry and
+ // distances, medium indeces and relative mass densities for all
+ // intersections, or, if this number is larger than the size n of
+ // isections, -1 (but the n intersections are still put in isections).
+ // If the position is outside, the method checks if the trajectory
+ // intersects the geometry and if yes, puts in the first element of
+ // isections the distance to the entry point and then finds all other
+ // intersections as in the case of x inside.
+
+ if (n < 1) {
+ return -1;
+ }
+
+ int isec_idx = 0; // current intersection index
+ EGS_Float t; // distance to next boundary
+ EGS_Float ttot = 0; // total distance along path
+ EGS_Vector x(X); // current position
+
+ int imed;
+
+ bool outside_envelope = ireg < 0;
+ bool in_base_geom = ireg >= 0 && ireg < nregbase;
+
+ if (outside_envelope) {
+ t = 1e30;
+
+ // region we would hit
+ ireg = howfar(ireg,x,u,t,&imed);
+ bool would_not_intersect = ireg < 0;
+
+ if (would_not_intersect) {
+ return 0;
+ }
+
+ isections[0].t = t; // distance to entry point
+ isections[0].rhof = 1;
+ isections[0].ireg = -1;
+ isections[0].imed = -1;
+ ttot = t;
+
+ x += u*t;
+
+ isec_idx = 1;
+ }
+ else {
+ // we're already inside the envelope
+ imed = medium(ireg);
+ }
+
+ // keep looping until we exit the geometry!
+ while (1) {
+ //egsInformation("in loop: j=%d ireg=%d imed=%d x=(%g,%g,%g)\n",
+ // j,ireg,imed,x.x,x.y,x.z);
+ //
+
+ // why are we setting these here?
+ // on first pass if we were outside envelope then:
+ // ireg is set to entry region
+ // imed is set to med of entry region (set in howfar above)
+ // else
+ // ireg is region of input ireg (current region?)
+ // imed is medium of input ireg (current region?)
+ isections[isec_idx].imed = imed; //
+ isections[isec_idx].ireg = ireg;
+ isections[isec_idx].rhof = getRelativeRho(ireg);
+
+ if (in_base_geom) { // in one of the regions of the base geometry
+
+ t = 1e30;
+
+ // next region will be in base geom by default
+ ireg = base_geom->howfar(ireg,x,u,t,&imed);
+
+ // if there's inscribed geoms in this region see if we will be intersecting one
+ vector geoms_in_reg = getGeomsInRegion(ireg);
+ for (vector::iterator geom=geoms_in_reg.begin(); geom!= geoms_in_reg.end(); ++geom) {
+
+ int inscribed_reg = (*geom)->howfar(-1, x, u, t, &imed);
+
+ bool hits_inscribed = inscribed_reg >= 0;
+
+ if (hits_inscribed) {
+ ireg = getGlobalRegFromLocalReg(*geom, inscribed_reg);
+ }
+ }
+
+ ttot += t;
+ isections[isec_idx++].t = ttot;
+
+ // left the geometry (implies that no inscribed geoms hit
+ if (ireg < 0) {
+ return isec_idx;
+ }
+
+
+ // over limit of number of intersections
+ if (isec_idx >= n) {
+ return -1;
+ }
+
+ // move to next boundary
+ x += u*t;
+
+ }
+ else {
+
+ // in inscribed geometry
+ volcor::GeomRegPairT local = getLocalFromGlobalReg(ireg);
+
+ int max_isec = n - isec_idx;
+
+ int ninter_sec = local.first->computeIntersections(local.second, max_isec,
+ x, u, &isections[isec_idx]);
+
+ // convert intersection region indexes to global indexex
+ int nmax = ninter_sec >= 0 ? ninter_sec + isec_idx : n;
+ for (int i=isec_idx; i < nmax; i++) {
+ isections[i].ireg = getGlobalRegFromLocalReg(local.first, isections[i].ireg);
+ isections[i].t += ttot;
+ }
+
+ //egsInformation("last intersection: %g\n",isections[nm-1].t);
+ if (ninter_sec < 0) {
+ return ninter_sec;
+ }
+
+ isec_idx += ninter_sec;
+
+ if (isec_idx >= n) {
+ return -1;
+ }
+
+ t = isections[isec_idx-1].t - ttot;
+ x += u*t;
+ ttot = isections[isec_idx-1].t;
+ ireg = base_geom->isWhere(x);
+ //egsInformation("new region: %d\n",ireg);
+
+ if (ireg < 0) {
+ return isec_idx;
+ }
+
+ imed = base_geom->medium(ireg);
+ }
+ }
+ return -1;
+}
+
+EGS_Float EGS_AEnvelope::howfarToOutside(int ireg, const EGS_Vector &x, const EGS_Vector &u) {
+
+ if (ireg < 0) {
+ return 0;
+ }
+
+ EGS_Float d;
+
+ bool in_base_geom = ireg < nregbase;
+ if (in_base_geom) {
+ d = base_geom->howfarToOutside(ireg, x, u);
+ }
+ else if (base_geom->regions() == 1) {
+ d = base_geom->howfarToOutside(0, x, u);
+ }
+ else {
+ int ir = base_geom->isWhere(x);
+ d = base_geom->howfarToOutside(ir,x,u);
+ }
+ return d;
+};
+
+
+int EGS_AEnvelope::howfar(int ireg, const EGS_Vector &x, const EGS_Vector &u,
+ EGS_Float &t, int *newmed, EGS_Vector *normal) {
+
+ bool inside_geom = ireg >= 0;
+ if (inside_geom) {
+
+ bool inside_base_geom = ireg < nregbase;
+
+ if (inside_base_geom) {
+
+ int base_reg = base_geom->howfar(ireg,x,u,t,newmed,normal);
+ vector geoms_in_reg = getGeomsInRegion(ireg);
+
+ bool no_inscribed_in_reg = geoms_in_reg.size() == 0;
+ if (no_inscribed_in_reg) {
+ return base_reg;
+ }
+
+ int inscribed_global = -1;
+ for (vector::iterator geom=geoms_in_reg.begin(); geom!= geoms_in_reg.end(); ++geom) {
+ int local_reg = (*geom)->howfar(-1, x, u, t, newmed, normal);
+ bool hits_inscribed = local_reg >= 0;
+ if (hits_inscribed) {
+ inscribed_global = getGlobalRegFromLocalReg(*geom, local_reg);
+ }
+ }
+
+ bool hit_inscribed_first = inscribed_global >= 0;
+
+ if (hit_inscribed_first) {
+ return inscribed_global;
+ }
+
+ return base_reg;
+ }
+ else {
+
+ // if here, we are in an inscribed geometry.
+ volcor::GeomRegPairT local = getLocalFromGlobalReg(ireg);
+
+ // and then check if we will hit a boundary in this geometry.
+ int new_local = local.first->howfar(local.second, x, u, t, newmed, normal);
+ bool hit_boundary_in_inscribed = new_local >=0;
+ if (hit_boundary_in_inscribed) {
+ return getGlobalRegFromLocalReg(local.first, new_local);
+ }
+
+ // new_local < 0 implies that we have exited the inscribed geometry
+ // => check to see in which base geometry region we are.
+ int inew = base_geom->isWhere(x + u*t);
+ if (inew >= 0 && newmed) {
+ *newmed = base_geom->medium(inew);
+ }
+ return inew;
+ }
+ }
+
+ // if here, we are outside the base geometry.
+ // check to see if we will enter.
+ int new_base_reg = base_geom->howfar(ireg,x,u,t,newmed,normal);
+ vector geoms_in_reg = getGeomsInRegion(new_base_reg);
+ bool has_geoms_in_reg = geoms_in_reg.size() > 0;
+
+ //cout << "new reg "<= 0 && has_geoms_in_reg) {
+
+ // check if we will be inside an inscribed geometry when we enter
+ for (vector::iterator geom=geoms_in_reg.begin(); geom!= geoms_in_reg.end(); ++geom) {
+ int new_local_reg = (*geom)->isWhere(x+u*t);
+ bool landed_in_inscribed = new_local_reg >= 0;
+ if (landed_in_inscribed) {
+ if (newmed) {
+ *newmed = (*geom)->medium(new_local_reg);
+ }
+ return getGlobalRegFromLocalReg(*geom, new_local_reg);
+ }
+
+ }
+ }
+
+ return new_base_reg;
+};
+
+
+EGS_Float EGS_AEnvelope::hownear(int ireg, const EGS_Vector &x) {
+ bool inside_envelope = ireg >= 0;
+
+ if (!inside_envelope) {
+ return base_geom->hownear(ireg, x);
+ }
+
+ bool in_base_geom = ireg < nregbase ;
+ if (!in_base_geom) {
+ volcor::GeomRegPairT local = getLocalFromGlobalReg(ireg);
+ return local.first->hownear(local.second, x);
+ }
+
+ EGS_Float tmin = base_geom->hownear(ireg, x);
+ if (tmin <= 0) {
+ return tmin;
+ }
+ vector geoms_in_reg = getGeomsInRegion(ireg);
+ for (vector::iterator geom=geoms_in_reg.begin(); geom!= geoms_in_reg.end(); ++geom) {
+ EGS_Float local_t = (*geom)->hownear(-1, x);
+ if (local_t < tmin) {
+ tmin = local_t;
+ if (tmin < 0) {
+ return tmin;
+ }
+ }
+ }
+ return tmin;
+};
+
+bool EGS_AEnvelope::hasBooleanProperty(int ireg, EGS_BPType prop) const {
+ if (ireg >= 0 && ireg < nreg) {
+ if (ireg < nregbase) {
+ return base_geom->hasBooleanProperty(ireg,prop);
+ }
+ volcor::GeomRegPairT local = getLocalFromGlobalReg(ireg);
+ return local.first->hasBooleanProperty(local.second, prop);
+ }
+ return false;
+};
+
+void EGS_AEnvelope::setBooleanProperty(EGS_BPType prop) {
+ setPropertyError("setBooleanProperty()");
+};
+
+void EGS_AEnvelope::addBooleanProperty(int bit) {
+ setPropertyError("addBooleanProperty()");
+};
+
+void EGS_AEnvelope::setBooleanProperty(EGS_BPType prop, int start, int end,int step) {
+ setPropertyError("setBooleanProperty()");
+};
+
+void EGS_AEnvelope::addBooleanProperty(int bit,int start,int end,int step) {
+ setPropertyError("addBooleanProperty()");
+};
+
+int EGS_AEnvelope::getMaxStep() const {
+ int nstep = base_geom->getMaxStep();
+ for (size_t j=0; j< inscribed_geoms.size(); ++j) {
+ nstep += inscribed_geoms[j]->getMaxStep();
+ }
+ return nstep + inscribed_geoms.size();
+};
+
+EGS_Float EGS_AEnvelope::getMass(int ireg) {
+
+ if (ireg < 0) {
+ return -1;
+ }
+ else if (ireg < nregbase) {
+ return corrected_mass[ireg];
+ }
+
+ volcor::GeomRegPairT local = global_reg_to_local[ireg];
+
+ return local.first->getMass(local.second);
+};
+
+EGS_Float EGS_AEnvelope::getMassCorrectionRatio(int ireg) {
+
+ if (0 <= ireg && ireg < nregbase) {
+ return corrected_mass[ireg]/uncorrected_mass[ireg];
+ }
+
+ return 1;
+};
+
+
+/* Print information about the geometry. If `print debug info = yes` is present
+ * in the geometry, extra information about which regions of the base geometry
+ * contain inscribed geometries and how their volumes were corrected.
+*/
+void EGS_AEnvelope::printInfo() const {
+
+ EGS_BaseGeometry::printInfo();
+
+ if (debug_info) {
+
+ for (int ir=0; ir < nregbase; ir++) {
+ if (fabs(uncorrected_mass[ir] - corrected_mass[ir]) > 1E-8) {
+ egsInformation(" mass of region %d was corrected from %.5E g to %.5E g\n",
+ ir, uncorrected_mass[ir], corrected_mass[ir]);
+ }
+ }
+
+ for (int ir=0; ir < nregbase; ir++) {
+ if (geoms_in_region[ir].size() > 0) {
+ egsInformation(" region %d has %d insribed geometries\n", ir, geoms_in_region[ir].size());
+ }
+ }
+ }
+
+ vc_results.outputResults();
+
+ egsInformation(" base geometry = %s (type %s) nreg=%d\n",base_geom->getName().c_str(),
+ base_geom->getType().c_str(), nregbase, base_geom->regions());
+ egsInformation(" # of inscribed geometries= %d in %d regions\n",inscribed_geoms.size(), getNRegWithInscribed());
+
+
+ egsInformation("=======================================================\n");
+}
+
+void EGS_AEnvelope::writeVCToFile(ostream &out) {
+
+ vector to_write;
+
+ for (int i=0; i < base_geom->regions(); i++) {
+ EGS_Float cor = vc_results.corrected_volumes[i];
+ EGS_Float uncor = vc_results.uncorrected_volumes[i];
+ bool has_correction = fabs(cor-uncor) > 1E-8;
+ if (has_correction) {
+ to_write.push_back(i);
+ }
+ }
+
+ size_t nrecords = to_write.size();
+ out << nrecords << endl;
+
+ for (size_t i=0; i geoms_in_reg = getGeomsInRegion(ir);
+ vector geom_idxs_in_reg;
+ for (size_t i=0; igetName();
+ string fname = gname+".autoenv.volcor";
+ fname += (output_vc == "gzip" ? ".gz" : "");
+ string mode = (output_vc == "gzip" ? "gzip" : "text");
+
+ egsInformation(
+ "Writing volume correction file for %s to %s file %s\n",
+ gname.c_str(), mode.c_str(), fname.c_str());
+
+ if (output_vc == "gzip") {
+#ifdef HAS_GZSTREAM
+ ogzstream outgz(fname.c_str());
+ writeVCToFile(outgz);
+ outgz.close();
+#endif
+ }
+ else {
+ ofstream out(fname.c_str());
+ writeVCToFile(out);
+ out.close();
+ }
+
+}
+
+
+/* check if the base or any of inscribed geometries have density scaling on */
+bool EGS_AEnvelope::getHasRhoScaling() {
+
+ has_rho_scaling = base_geom->hasRhoScaling();
+ if (has_rho_scaling) {
+ return true;
+ }
+
+ for (size_t geom_idx=0; geom_idx < inscribed_geoms.size(); geom_idx++) {
+ if (inscribed_geoms[geom_idx]->hasRhoScaling()) {
+ return true;
+ }
+ }
+ return false;
+
+}
+
+
+void EGS_AEnvelope::setMedia(EGS_Input *,int,const int *) {
+ egsWarning("EGS_AEnvelope::setMedia: don't use this method. Use the\n"
+ " setMedia() methods of the geometry objects that make up this geometry\n");
+}
+
+
+void EGS_AEnvelope::setRelativeRho(int start, int end, EGS_Float rho) {
+ setRelativeRho(0);
+}
+
+
+void EGS_AEnvelope::setRelativeRho(EGS_Input *) {
+ egsWarning("EGS_AEnvelope::setRelativeRho(): don't use this method."
+ " Use the\n setRelativeRho methods of the geometry objects that make up"
+ " this geometry\n");
+}
+
+
+EGS_Float EGS_AEnvelope::getRelativeRho(int ireg) const {
+ if (ireg < 0 || ireg >= nreg) {
+ return 1;
+ }
+ if (ireg < nregbase) {
+ EGS_Float v = base_geom->getRelativeRho(ireg);
+ return v;
+ }
+ volcor::GeomRegPairT local = getLocalFromGlobalReg(ireg);
+ if (local.first) {
+ return local.first->getRelativeRho(local.second);
+ }
+ return 1;
+};
+
+
+/*************************************************************************/
+/* Switched envelope *****************************************************/
+/*************************************************************************/
+
+EGS_ASwitchedEnvelope::EGS_ASwitchedEnvelope(EGS_BaseGeometry *base_geom,
+ const vector inscribed, const string &Name, bool debug, string output_vc_file):
+ EGS_AEnvelope(base_geom, inscribed, Name, debug, output_vc_file) {
+
+ // activate a single geometry
+ cur_ptr = 0;
+ active_inscribed.push_back(inscribed_geoms[cur_ptr]);
+
+};
+
+
+//TODO: this gets called a lot and is probably quite slow. Instead fo doing a
+//set intersection on every call we can probably do it once when activated
+//geometries change and cache it
+vector EGS_ASwitchedEnvelope::getGeomsInRegion(int ireg) {
+ // find which geometries are present AND active in this region
+ vector active_in_reg;
+ if ((0 <= ireg) && (ireg < nregbase)) {
+ set_intersection(
+ geoms_in_region[ireg].begin(), geoms_in_region[ireg].end(),
+ active_inscribed.begin(), active_inscribed.end(),
+ back_inserter(active_in_reg)
+ );
+ }
+
+ return active_in_reg;
+
+}
+
+
+void EGS_ASwitchedEnvelope::setActiveGeometries(vector geoms) {
+ active_inscribed.clear();
+ active_inscribed = geoms;
+ cur_ptr = -1;
+ if (geoms.size() <= 0) {
+ return;
+ }
+
+ for (size_t i=0; i < inscribed_geoms.size(); i++) {
+ if (inscribed_geoms[i] == geoms[0]) {
+ cur_ptr = i;
+ return;
+ }
+ }
+}
+
+
+void EGS_ASwitchedEnvelope::setActiveGeometries(vector geom_indexes) {
+
+ active_inscribed.clear();
+
+ for (size_t i = 0; i < geom_indexes.size(); i++) {
+ int idx = geom_indexes[i];
+ if (idx < 0 || idx >= ninscribed) {
+ egsFatal("EGS_ASwitchedEnvelope:: %d is not a valid geometry index\n", idx);
+ }
+ active_inscribed.push_back(inscribed_geoms[i]);
+ }
+
+ cur_ptr = geom_indexes[0];
+
+}
+
+
+bool EGS_ASwitchedEnvelope::hasActiveGeom(int ireg) {
+
+ if ((0 <= ireg) && (ireg < nregbase)) {
+ size_t nactive = getGeomsInRegion(ireg).size();
+ return nactive > 0;
+ }
+
+ return false;
+}
+
+
+bool EGS_ASwitchedEnvelope::hasInactiveGeom(int ireg) {
+
+ if ((0 <= ireg) && (ireg < nregbase)) {
+ size_t nactive = getGeomsInRegion(ireg).size();
+ size_t ntotal = EGS_AEnvelope::getGeomsInRegion(ireg).size();
+ return nactive < ntotal;
+ }
+
+ return false;
+
+}
+
+void EGS_ASwitchedEnvelope::setActiveByIndex(int inscribed_index) {
+ vector to_activate;
+ to_activate.push_back(inscribed_geoms[inscribed_index]);
+ setActiveGeometries(to_activate);
+ cur_ptr = inscribed_index;
+}
+
+void EGS_ASwitchedEnvelope::activateByIndex(int inscribed_index) {
+ active_inscribed.push_back(inscribed_geoms[inscribed_index]);
+}
+
+void EGS_ASwitchedEnvelope::deactivateByIndex(int inscribed_index) {
+ vector::iterator loc = find(
+ active_inscribed.begin(), active_inscribed.end(),
+ inscribed_geoms[inscribed_index]
+ );
+ if (loc != active_inscribed.end()) {
+ active_inscribed.erase(loc);
+ }
+}
+
+void EGS_ASwitchedEnvelope::cycleActive() {
+ cur_ptr = cur_ptr == ninscribed -1 ? 0 : cur_ptr + 1;
+ vector to_activate;
+ to_activate.push_back(inscribed_geoms[cur_ptr]);
+ setActiveGeometries(to_activate);
+}
+
+
+vector EGS_AEnvelope::createTransforms(EGS_Input *input) {
+
+ vector transforms;
+ if (input) {
+ EGS_Input *i;
+
+ while ((i = input->takeInputItem(transformation_keyword))) {
+ EGS_AffineTransform *transform = EGS_AffineTransform::getTransformation(i);
+ if (!transform) {
+ egsWarning("Invalid transform input given\n");
+
+ }
+ transforms.push_back(transform);
+ delete i;
+
+ }
+ }
+
+ return transforms;
+
+}
+
+
+extern "C" {
+
+ EGS_AENVELOPE_EXPORT EGS_BaseGeometry *createGeometry(EGS_Input *input) {
+
+ if (!input) {
+ egsWarning(geom_class_msg, "null input");
+ return 0;
+ }
+
+ bool debug;
+ vector debug_choices;
+ debug_choices.push_back("no");
+ debug_choices.push_back("yes");
+ debug = input->getInput("print debug info", debug_choices, 0);
+
+ int output_vc_file_choice;
+ vector vc_file_choices;
+ vc_file_choices.push_back("no");
+ vc_file_choices.push_back("yes");
+ vc_file_choices.push_back("text");
+ vc_file_choices.push_back("gzip");
+ output_vc_file_choice = input->getInput("output volume correction file", vc_file_choices, 0);
+ string output_vc_file = vc_file_choices[output_vc_file_choice];
+
+#ifndef HAS_GZSTREAM
+ if (output_vc_file == "gzip") {
+ egsWarning(
+ "GZip file output requested but not compiled with gzstream.\n"
+ "Please recompile with gzstream support.\n"
+ );
+ return 0;
+ }
+#endif
+
+
+
+ string base_geom_name;
+ int err = input->getInput(base_geom_keyword, base_geom_name);
+ if (err) {
+ egsWarning(geom_class_msg, ("'"+string(base_geom_keyword)+"' input not found").c_str());
+ return 0;
+ }
+
+ string type;
+ err = input->getInput(type_keyword, type);
+ if (err) {
+ type = "AEnvelope";
+ }
+
+ EGS_BaseGeometry *base_geom = EGS_BaseGeometry::getGeometry(base_geom_name);
+ if (!base_geom) {
+ egsWarning(geom_class_msg, ("Unable to find geometry '"+base_geom_name+"'").c_str());
+ return 0;
+ }
+
+ EGS_Input *inscribed_input = input->takeInputItem(inscribed_geom_keyword);
+ if (!inscribed_input) {
+ egsWarning(geom_class_msg, ("Missing '"+string(inscribed_geom_keyword)+"' input item").c_str());
+ return 0;
+ }
+
+ string inscribed_geom_name;
+ err = inscribed_input->getInput(inscribed_geom_name_keyword, inscribed_geom_name);
+ if (err) {
+ egsWarning(geom_class_msg, ("'"+string(inscribed_geom_name_keyword)+"' input not found").c_str());
+ return 0;
+ }
+
+ EGS_BaseGeometry *inscribed_geom = EGS_BaseGeometry::getGeometry(inscribed_geom_name);
+ if (!inscribed_geom) {
+ egsWarning(geom_class_msg, ("Unable to find geometry '"+inscribed_geom_name+"'").c_str());
+ return 0;
+ }
+
+ vector transforms;
+ EGS_Input *trans_input = inscribed_input->takeInputItem(transformations_keyword);
+ if (trans_input) {
+ transforms = EGS_AEnvelope::createTransforms(trans_input);
+ }
+
+ delete trans_input;
+
+ EGS_Input *volcor_input = inscribed_input->takeInputItem("region discovery");
+ VCOptions *vcopts = new VCOptions(volcor_input);
+ if (!vcopts->valid) {
+ egsWarning(geom_class_msg, "Missing or invalid 'region discovery' input item");
+ return 0;
+ }
+
+ delete inscribed_input;
+
+
+ vector inscribed;
+ if (transforms.size()>0) {
+ for (size_t i=0; i < transforms.size(); i++) {
+ inscribed.push_back(AEnvelopeAux(inscribed_geom, transforms[i], vcopts));
+ }
+ }
+ else {
+ EGS_AffineTransform *unityt = new EGS_AffineTransform();
+ inscribed.push_back(AEnvelopeAux(inscribed_geom, unityt, vcopts));
+ }
+
+ EGS_BaseGeometry *result;
+ if (type == "EGS_ASwitchedEnvelope") {
+ result = new EGS_ASwitchedEnvelope(base_geom, inscribed, "", (bool)debug, output_vc_file);
+ }
+ else {
+ result = new EGS_AEnvelope(base_geom, inscribed, "", (bool)debug, output_vc_file);
+ }
+
+ result->setName(input);
+ return result;
+ }
+
+}
diff --git a/HEN_HOUSE/egs++/geometry/egs_autoenvelope/egs_autoenvelope.h b/HEN_HOUSE/egs++/geometry/egs_autoenvelope/egs_autoenvelope.h
new file mode 100644
index 000000000..915211e2a
--- /dev/null
+++ b/HEN_HOUSE/egs++/geometry/egs_autoenvelope/egs_autoenvelope.h
@@ -0,0 +1,603 @@
+/*
+###############################################################################
+#
+# EGSnrc egs++ auto envelope geometry headers
+# Copyright (C) 2016 Randle E. P. Taylor, Rowan M. Thomson,
+# Marc J. P. Chamberland, D. W. O. Rogers
+#
+# This file is part of EGSnrc.
+#
+# EGSnrc is free software: you can redistribute it and/or modify it under
+# the terms of the GNU Affero General Public License as published by the
+# Free Software Foundation, either version 3 of the License, or (at your
+# option) any later version.
+#
+# EGSnrc is distributed in the hope that it will be useful, but WITHOUT ANY
+# WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+# FOR A PARTICULAR PURPOSE. See the GNU Affero General Public License for
+# more details.
+#
+# You should have received a copy of the GNU Affero General Public License
+# along with EGSnrc. If not, see .
+#
+###############################################################################
+#
+# Author: Randle Taylor, 2016
+#
+# Contributors: Marc Chamberland
+# Rowan Thomson
+# Dave Rogers
+#
+###############################################################################
+#
+# egs_autoenvelope was developed for the Carleton Laboratory for
+# Radiotherapy Physics.
+#
+###############################################################################
+*/
+
+
+/*! \file egs_autoenvelope.h
+ * \brief An envelope geometry with automatic inscribed region detection (inspired by EGS_FastEnvelope)
+ * \author Randle Taylor
+ */
+
+#ifndef EGS_AENVELOPE_GEOMETRY_
+#define EGS_AENVELOPE_GEOMETRY_
+
+#ifdef WIN32
+
+ #ifdef BUILD_AENVELOPE_DLL
+ #define EGS_AENVELOPE_EXPORT __declspec(dllexport)
+ #else
+ #define EGS_AENVELOPE_EXPORT __declspec(dllimport)
+ #endif
+ #define EGS_AENVELOPE_LOCAL
+
+#else
+
+ #ifdef HAVE_VISIBILITY
+ #define EGS_AENVELOPE_EXPORT __attribute__ ((visibility ("default")))
+ #define EGS_AENVELOPE_LOCAL __attribute__ ((visibility ("hidden")))
+ #else
+ #define EGS_AENVELOPE_EXPORT
+ #define EGS_AENVELOPE_LOCAL
+ #endif
+
+#endif
+
+
+#include "egs_base_geometry.h"
+#include "egs_functions.h"
+#include "egs_transformations.h"
+#include "egs_shapes.h"
+#include "volcor.h"
+
+#include
+#include
+#include