diff --git a/BONDI_SOLUTION.dat b/BONDI_SOLUTION.dat deleted file mode 100644 index f5f7b52..0000000 --- a/BONDI_SOLUTION.dat +++ /dev/null @@ -1,1000 +0,0 @@ -2.18776e-07 1.74699e-08 2.18984e+10 4275.92 -2.23872e-07 1.74752e-08 2.1155e+10 4226.97 -2.29087e-07 1.74808e-08 2.04368e+10 4178.59 -2.34423e-07 1.74865e-08 1.9743e+10 4130.76 -2.39883e-07 1.74925e-08 1.90727e+10 4083.47 -2.45471e-07 1.74987e-08 1.84252e+10 4036.73 -2.51189e-07 1.75051e-08 1.77997e+10 3990.52 -2.57039e-07 1.75117e-08 1.71954e+10 3944.84 -2.63027e-07 1.75185e-08 1.66116e+10 3899.68 -2.69153e-07 1.75256e-08 1.60477e+10 3855.04 -2.75423e-07 1.75329e-08 1.55029e+10 3810.92 -2.81838e-07 1.75405e-08 1.49766e+10 3767.29 -2.88403e-07 1.75484e-08 1.44681e+10 3724.17 -2.95121e-07 1.75565e-08 1.3977e+10 3681.54 -3.01995e-07 1.75649e-08 1.35024e+10 3639.4 -3.09029e-07 1.75736e-08 1.30441e+10 3597.74 -3.16228e-07 1.75826e-08 1.26012e+10 3556.55 -3.23594e-07 1.7592e-08 1.21734e+10 3515.84 -3.31131e-07 1.76016e-08 1.17601e+10 3475.6 -3.38844e-07 1.76116e-08 1.13609e+10 3435.81 -3.46737e-07 1.7622e-08 1.09752e+10 3396.48 -3.54814e-07 1.76327e-08 1.06026e+10 3357.6 -3.63078e-07 1.76438e-08 1.02427e+10 3319.17 -3.71535e-07 1.76553e-08 9.89492e+09 3281.17 -3.80189e-07 1.76672e-08 9.55901e+09 3243.61 -3.89045e-07 1.76795e-08 9.23448e+09 3206.48 -3.98107e-07 1.76922e-08 8.92099e+09 3169.78 -4.0738e-07 1.77054e-08 8.61812e+09 3133.49 -4.16869e-07 1.77191e-08 8.32555e+09 3097.63 -4.2658e-07 1.77332e-08 8.0429e+09 3062.17 -4.36516e-07 1.77478e-08 7.76986e+09 3027.12 -4.46684e-07 1.7763e-08 7.50607e+09 2992.46 -4.57088e-07 1.77786e-08 7.25125e+09 2958.21 -4.67735e-07 1.77949e-08 7.00507e+09 2924.35 -4.7863e-07 1.78117e-08 6.76726e+09 2890.87 -4.89779e-07 1.7829e-08 6.53752e+09 2857.78 -5.01187e-07 1.7847e-08 6.31558e+09 2825.07 -5.12861e-07 1.78657e-08 6.10117e+09 2792.73 -5.24807e-07 1.7885e-08 5.89404e+09 2760.76 -5.37032e-07 1.79049e-08 5.69394e+09 2729.16 -5.49541e-07 1.79256e-08 5.50064e+09 2697.92 -5.62341e-07 1.7947e-08 5.3139e+09 2667.03 -5.7544e-07 1.79691e-08 5.13349e+09 2636.5 -5.88844e-07 1.7992e-08 4.95922e+09 2606.33 -6.0256e-07 1.80157e-08 4.79085e+09 2576.49 -6.16595e-07 1.80403e-08 4.62821e+09 2547 -6.30958e-07 1.80657e-08 4.47109e+09 2517.84 -6.45654e-07 1.8092e-08 4.3193e+09 2489.02 -6.60694e-07 1.81193e-08 4.17266e+09 2460.53 -6.76083e-07 1.81475e-08 4.03101e+09 2432.36 -6.91831e-07 1.81766e-08 3.89415e+09 2404.52 -7.07946e-07 1.82068e-08 3.76195e+09 2377 -7.24436e-07 1.82381e-08 3.63424e+09 2349.79 -7.4131e-07 1.82705e-08 3.51086e+09 2322.89 -7.58578e-07 1.8304e-08 3.39167e+09 2296.3 -7.76247e-07 1.83387e-08 3.27653e+09 2270.01 -7.94328e-07 1.83746e-08 3.16529e+09 2244.03 -8.1283e-07 1.84118e-08 3.05783e+09 2218.34 -8.31764e-07 1.84502e-08 2.95402e+09 2192.95 -8.51138e-07 1.849e-08 2.85374e+09 2167.84 -8.70964e-07 1.85313e-08 2.75685e+09 2143.03 -8.91251e-07 1.85739e-08 2.66326e+09 2118.5 -9.12011e-07 1.86181e-08 2.57285e+09 2094.25 -9.33254e-07 1.86638e-08 2.4855e+09 2070.27 -9.54993e-07 1.87112e-08 2.40112e+09 2046.58 -9.77237e-07 1.87602e-08 2.31961e+09 2023.15 -1e-06 1.88109e-08 2.24086e+09 1999.99 -1.02329e-06 1.88634e-08 2.16478e+09 1977.09 -1.04713e-06 1.89177e-08 2.09129e+09 1954.46 -1.07152e-06 1.8974e-08 2.02029e+09 1932.09 -1.09648e-06 1.90322e-08 1.95171e+09 1909.97 -1.12202e-06 1.90925e-08 1.88545e+09 1888.11 -1.14815e-06 1.91549e-08 1.82144e+09 1866.5 -1.1749e-06 1.92195e-08 1.7596e+09 1845.13 -1.20226e-06 1.92863e-08 1.69987e+09 1824.01 -1.23027e-06 1.93555e-08 1.64216e+09 1803.13 -1.25893e-06 1.94272e-08 1.58641e+09 1782.49 -1.28825e-06 1.95013e-08 1.53255e+09 1762.09 -1.31826e-06 1.95781e-08 1.48052e+09 1741.92 -1.34896e-06 1.96576e-08 1.43026e+09 1721.98 -1.38038e-06 1.97398e-08 1.3817e+09 1702.26 -1.41254e-06 1.9825e-08 1.3348e+09 1682.78 -1.44544e-06 1.99131e-08 1.28948e+09 1663.52 -1.47911e-06 2.00043e-08 1.2457e+09 1644.47 -1.51356e-06 2.00988e-08 1.20341e+09 1625.65 -1.54882e-06 2.01965e-08 1.16256e+09 1607.04 -1.58489e-06 2.02977e-08 1.12309e+09 1588.64 -1.62181e-06 2.04025e-08 1.08496e+09 1570.46 -1.65959e-06 2.05109e-08 1.04813e+09 1552.48 -1.69824e-06 2.06231e-08 1.01255e+09 1534.71 -1.7378e-06 2.07393e-08 9.78174e+08 1517.14 -1.77828e-06 2.08596e-08 9.44965e+08 1499.77 -1.8197e-06 2.09841e-08 9.12884e+08 1482.61 -1.86209e-06 2.1113e-08 8.81893e+08 1465.64 -1.90546e-06 2.12464e-08 8.51953e+08 1448.86 -1.94984e-06 2.13844e-08 8.23031e+08 1432.27 -1.99526e-06 2.15274e-08 7.95089e+08 1415.88 -2.04174e-06 2.16753e-08 7.68097e+08 1399.67 -2.0893e-06 2.18285e-08 7.42021e+08 1383.65 -2.13796e-06 2.19871e-08 7.16831e+08 1367.81 -2.18776e-06 2.21512e-08 6.92494e+08 1352.15 -2.23872e-06 2.23211e-08 6.68986e+08 1336.67 -2.29087e-06 2.24969e-08 6.46274e+08 1321.37 -2.34423e-06 2.26789e-08 6.24334e+08 1306.25 -2.39883e-06 2.28674e-08 6.03138e+08 1291.29 -2.45471e-06 2.30624e-08 5.82663e+08 1276.51 -2.51189e-06 2.32643e-08 5.62882e+08 1261.9 -2.57039e-06 2.34733e-08 5.43773e+08 1247.45 -2.63027e-06 2.36897e-08 5.25312e+08 1233.17 -2.69153e-06 2.39136e-08 5.07479e+08 1219.06 -2.75423e-06 2.41455e-08 4.9025e+08 1205.1 -2.81838e-06 2.43854e-08 4.73607e+08 1191.31 -2.88403e-06 2.46338e-08 4.57528e+08 1177.67 -2.95121e-06 2.4891e-08 4.41996e+08 1164.19 -3.01995e-06 2.51571e-08 4.2699e+08 1150.86 -3.09029e-06 2.54327e-08 4.12495e+08 1137.69 -3.16228e-06 2.57179e-08 3.98491e+08 1124.67 -3.23594e-06 2.60131e-08 3.84962e+08 1111.79 -3.31131e-06 2.63187e-08 3.71894e+08 1099.06 -3.38844e-06 2.6635e-08 3.59268e+08 1086.48 -3.46737e-06 2.69625e-08 3.47072e+08 1074.05 -3.54814e-06 2.73015e-08 3.35289e+08 1061.75 -3.63078e-06 2.76523e-08 3.23906e+08 1049.6 -3.71535e-06 2.80156e-08 3.1291e+08 1037.58 -3.80189e-06 2.83915e-08 3.02287e+08 1025.7 -3.89045e-06 2.87807e-08 2.92025e+08 1013.96 -3.98107e-06 2.91836e-08 2.82111e+08 1002.36 -4.0738e-06 2.96006e-08 2.72534e+08 990.881 -4.16869e-06 3.00323e-08 2.63282e+08 979.538 -4.2658e-06 3.04791e-08 2.54343e+08 968.325 -4.36516e-06 3.09417e-08 2.45709e+08 957.24 -4.46684e-06 3.14205e-08 2.37367e+08 946.282 -4.57088e-06 3.19161e-08 2.29309e+08 935.45 -4.67735e-06 3.24291e-08 2.21524e+08 924.741 -4.7863e-06 3.29602e-08 2.14004e+08 914.156 -4.89779e-06 3.351e-08 2.06739e+08 903.691 -5.01187e-06 3.4079e-08 1.9972e+08 893.346 -5.12861e-06 3.46681e-08 1.9294e+08 883.119 -5.24807e-06 3.52778e-08 1.8639e+08 873.01 -5.37032e-06 3.5909e-08 1.80062e+08 863.016 -5.49541e-06 3.65624e-08 1.7395e+08 853.137 -5.62341e-06 3.72387e-08 1.68044e+08 843.371 -5.7544e-06 3.79389e-08 1.62339e+08 833.716 -5.88844e-06 3.86635e-08 1.56828e+08 824.172 -6.0256e-06 3.94137e-08 1.51504e+08 814.737 -6.16595e-06 4.01903e-08 1.46361e+08 805.411 -6.30958e-06 4.09941e-08 1.41392e+08 796.191 -6.45654e-06 4.18262e-08 1.36592e+08 787.076 -6.60694e-06 4.26875e-08 1.31955e+08 778.066 -6.76083e-06 4.35791e-08 1.27475e+08 769.159 -6.91831e-06 4.4502e-08 1.23147e+08 760.354 -7.07946e-06 4.54573e-08 1.18967e+08 751.65 -7.24436e-06 4.64463e-08 1.14928e+08 743.045 -7.4131e-06 4.74699e-08 1.11027e+08 734.54 -7.58578e-06 4.85296e-08 1.07257e+08 726.131 -7.76247e-06 4.96265e-08 1.03616e+08 717.818 -7.94328e-06 5.07619e-08 1.00098e+08 709.601 -8.1283e-06 5.19372e-08 9.67004e+07 701.478 -8.31764e-06 5.31539e-08 9.34175e+07 693.447 -8.51138e-06 5.44133e-08 9.02462e+07 685.509 -8.70964e-06 5.5717e-08 8.71824e+07 677.661 -8.91251e-06 5.70664e-08 8.42228e+07 669.904 -9.12011e-06 5.84634e-08 8.13635e+07 662.235 -9.33254e-06 5.99093e-08 7.86014e+07 654.654 -9.54993e-06 6.14062e-08 7.5933e+07 647.159 -9.77237e-06 6.29556e-08 7.33553e+07 639.751 -1e-05 6.45595e-08 7.08649e+07 632.427 -1.02329e-05 6.62197e-08 6.84591e+07 625.187 -1.04713e-05 6.79383e-08 6.61351e+07 618.03 -1.07152e-05 6.97173e-08 6.38899e+07 610.955 -1.09648e-05 7.15587e-08 6.1721e+07 603.961 -1.12202e-05 7.3465e-08 5.96257e+07 597.046 -1.14815e-05 7.54381e-08 5.76015e+07 590.212 -1.1749e-05 7.74807e-08 5.5646e+07 583.455 -1.20226e-05 7.9595e-08 5.3757e+07 576.776 -1.23027e-05 8.17837e-08 5.1932e+07 570.172 -1.25893e-05 8.40492e-08 5.01691e+07 563.645 -1.28825e-05 8.63944e-08 4.84659e+07 557.192 -1.31826e-05 8.88219e-08 4.68206e+07 550.814 -1.34896e-05 9.13349e-08 4.52311e+07 544.508 -1.38038e-05 9.3936e-08 4.36957e+07 538.274 -1.41254e-05 9.66287e-08 4.22123e+07 532.112 -1.44544e-05 9.94159e-08 4.07793e+07 526.02 -1.47911e-05 1.02301e-07 3.93949e+07 519.998 -1.51356e-05 1.05288e-07 3.80575e+07 514.045 -1.54882e-05 1.08379e-07 3.67655e+07 508.16 -1.58489e-05 1.11579e-07 3.55175e+07 502.343 -1.62181e-05 1.14892e-07 3.43117e+07 496.592 -1.65959e-05 1.18321e-07 3.31469e+07 490.907 -1.69824e-05 1.21871e-07 3.20216e+07 485.286 -1.7378e-05 1.25545e-07 3.09346e+07 479.731 -1.77828e-05 1.29349e-07 2.98844e+07 474.238 -1.8197e-05 1.33286e-07 2.88699e+07 468.809 -1.86209e-05 1.37361e-07 2.78899e+07 463.442 -1.90546e-05 1.4158e-07 2.69431e+07 458.136 -1.94984e-05 1.45947e-07 2.60284e+07 452.891 -1.99526e-05 1.50468e-07 2.51448e+07 447.706 -2.04174e-05 1.55147e-07 2.42912e+07 442.581 -2.0893e-05 1.59991e-07 2.34666e+07 437.513 -2.13796e-05 1.65005e-07 2.267e+07 432.505 -2.18776e-05 1.70195e-07 2.19004e+07 427.553 -2.23872e-05 1.75568e-07 2.11569e+07 422.658 -2.29087e-05 1.81129e-07 2.04387e+07 417.819 -2.34423e-05 1.86886e-07 1.97449e+07 413.035 -2.39883e-05 1.92846e-07 1.90746e+07 408.307 -2.45471e-05 1.99014e-07 1.84271e+07 403.632 -2.51189e-05 2.054e-07 1.78015e+07 399.011 -2.57039e-05 2.12009e-07 1.71972e+07 394.442 -2.63027e-05 2.18852e-07 1.66134e+07 389.926 -2.69153e-05 2.25934e-07 1.60495e+07 385.462 -2.75423e-05 2.33266e-07 1.55046e+07 381.049 -2.81838e-05 2.40855e-07 1.49783e+07 376.686 -2.88403e-05 2.48711e-07 1.44698e+07 372.373 -2.95121e-05 2.56843e-07 1.39786e+07 368.11 -3.01995e-05 2.65261e-07 1.35041e+07 363.895 -3.09029e-05 2.73975e-07 1.30457e+07 359.729 -3.16228e-05 2.82995e-07 1.26028e+07 355.61 -3.23594e-05 2.92332e-07 1.2175e+07 351.538 -3.31131e-05 3.01997e-07 1.17617e+07 347.513 -3.38844e-05 3.12002e-07 1.13624e+07 343.534 -3.46737e-05 3.22358e-07 1.09767e+07 339.601 -3.54814e-05 3.33079e-07 1.06041e+07 335.713 -3.63078e-05 3.44176e-07 1.02441e+07 331.869 -3.71535e-05 3.55663e-07 9.89638e+06 328.069 -3.80189e-05 3.67554e-07 9.56045e+06 324.313 -3.89045e-05 3.79863e-07 9.2359e+06 320.599 -3.98107e-05 3.92604e-07 8.92239e+06 316.928 -4.0738e-05 4.05793e-07 8.6195e+06 313.299 -4.16869e-05 4.19446e-07 8.32691e+06 309.712 -4.2658e-05 4.33579e-07 8.04424e+06 306.166 -4.36516e-05 4.48208e-07 7.77118e+06 302.66 -4.46684e-05 4.63352e-07 7.50738e+06 299.194 -4.57088e-05 4.79027e-07 7.25254e+06 295.768 -4.67735e-05 4.95254e-07 7.00635e+06 292.382 -4.7863e-05 5.1205e-07 6.76852e+06 289.034 -4.89779e-05 5.29438e-07 6.53876e+06 285.724 -5.01187e-05 5.47436e-07 6.3168e+06 282.452 -5.12861e-05 5.66067e-07 6.10237e+06 279.218 -5.24807e-05 5.85353e-07 5.89523e+06 276.02 -5.37032e-05 6.05317e-07 5.69511e+06 272.86 -5.49541e-05 6.25981e-07 5.5018e+06 269.735 -5.62341e-05 6.47373e-07 5.31504e+06 266.646 -5.7544e-05 6.69517e-07 5.13462e+06 263.593 -5.88844e-05 6.92438e-07 4.96033e+06 260.574 -6.0256e-05 7.16166e-07 4.79195e+06 257.59 -6.16595e-05 7.40727e-07 4.62929e+06 254.64 -6.30958e-05 7.66152e-07 4.47215e+06 251.724 -6.45654e-05 7.9247e-07 4.32035e+06 248.842 -6.60694e-05 8.19713e-07 4.1737e+06 245.992 -6.76083e-05 8.47913e-07 4.03203e+06 243.175 -6.91831e-05 8.77106e-07 3.89516e+06 240.39 -7.07946e-05 9.07323e-07 3.76295e+06 237.637 -7.24436e-05 9.38604e-07 3.63522e+06 234.915 -7.4131e-05 9.70982e-07 3.51183e+06 232.225 -7.58578e-05 1.0045e-06 3.39262e+06 229.565 -7.76247e-05 1.0392e-06 3.27747e+06 226.936 -7.94328e-05 1.07511e-06 3.16622e+06 224.337 -8.1283e-05 1.11229e-06 3.05875e+06 221.768 -8.31764e-05 1.15077e-06 2.95492e+06 219.228 -8.51138e-05 1.19061e-06 2.85463e+06 216.717 -8.70964e-05 1.23185e-06 2.75773e+06 214.235 -8.91251e-05 1.27453e-06 2.66413e+06 211.781 -9.12011e-05 1.31872e-06 2.5737e+06 209.355 -9.33254e-05 1.36446e-06 2.48634e+06 206.957 -9.54993e-05 1.41181e-06 2.40195e+06 204.587 -9.77237e-05 1.46082e-06 2.32043e+06 202.243 -0.0001 1.51155e-06 2.24166e+06 199.927 -0.000102329 1.56407e-06 2.16558e+06 197.637 -0.000104713 1.61844e-06 2.09207e+06 195.373 -0.000107152 1.67471e-06 2.02107e+06 193.135 -0.000109648 1.73297e-06 1.95247e+06 190.923 -0.000112202 1.79327e-06 1.8862e+06 188.736 -0.000114815 1.85569e-06 1.82218e+06 186.574 -0.00011749 1.9203e-06 1.76033e+06 184.436 -0.000120226 1.98719e-06 1.70059e+06 182.324 -0.000123027 2.05642e-06 1.64287e+06 180.235 -0.000125893 2.12809e-06 1.58711e+06 178.17 -0.000128825 2.20228e-06 1.53324e+06 176.129 -0.000131826 2.27908e-06 1.48121e+06 174.111 -0.000134896 2.35858e-06 1.43093e+06 172.116 -0.000138038 2.44087e-06 1.38237e+06 170.145 -0.000141254 2.52606e-06 1.33545e+06 168.195 -0.000144544 2.61424e-06 1.29013e+06 166.268 -0.000147911 2.70552e-06 1.24634e+06 164.363 -0.000151356 2.80001e-06 1.20404e+06 162.48 -0.000154882 2.89782e-06 1.16318e+06 160.618 -0.000158489 2.99906e-06 1.1237e+06 158.778 -0.000162181 3.10387e-06 1.08557e+06 156.959 -0.000165959 3.21236e-06 1.04872e+06 155.16 -0.000169824 3.32467e-06 1.01313e+06 153.382 -0.00017378 3.44092e-06 978751 151.625 -0.000177828 3.56126e-06 945535 149.887 -0.00018197 3.68583e-06 913446 148.169 -0.000186209 3.81478e-06 882447 146.471 -0.000190546 3.94827e-06 852500 144.793 -0.000194984 4.08644e-06 823570 143.134 -0.000199526 4.22948e-06 795621 141.493 -0.000204174 4.37754e-06 768621 139.872 -0.00020893 4.53081e-06 742538 138.268 -0.000213796 4.68947e-06 717340 136.684 -0.000218776 4.8537e-06 692997 135.117 -0.000223872 5.02371e-06 669481 133.568 -0.000229087 5.1997e-06 646763 132.037 -0.000234423 5.38187e-06 624815 130.524 -0.000239883 5.57045e-06 603613 129.028 -0.000245471 5.76566e-06 583130 127.549 -0.000251189 5.96773e-06 563343 126.086 -0.00025704 6.17691e-06 544228 124.641 -0.000263027 6.39344e-06 525761 123.212 -0.000269154 6.61759e-06 507921 121.8 -0.000275423 6.84961e-06 490686 120.403 -0.000281838 7.0898e-06 474037 119.023 -0.000288403 7.33843e-06 457952 117.658 -0.000295121 7.5958e-06 442414 116.309 -0.000301995 7.86222e-06 427403 114.975 -0.00030903 8.13801e-06 412901 113.657 -0.000316228 8.4235e-06 398892 112.354 -0.000323594 8.71902e-06 385358 111.065 -0.000331131 9.02494e-06 372283 109.791 -0.000338844 9.34162e-06 359652 108.532 -0.000346737 9.66943e-06 347450 107.288 -0.000354813 1.00088e-05 335662 106.057 -0.000363078 1.036e-05 324274 104.841 -0.000371535 1.07237e-05 313273 103.638 -0.000380189 1.11001e-05 302645 102.449 -0.000389045 1.14897e-05 292378 101.274 -0.000398107 1.18931e-05 282459 100.112 -0.00040738 1.23106e-05 272877 98.9637 -0.000416869 1.27428e-05 263620 97.8283 -0.00042658 1.31903e-05 254676 96.7058 -0.000436516 1.36534e-05 246037 95.5963 -0.000446684 1.41329e-05 237691 94.4993 -0.000457088 1.46292e-05 229628 93.415 -0.000467735 1.5143e-05 221839 92.343 -0.00047863 1.56748e-05 214314 91.2833 -0.000489779 1.62254e-05 207044 90.2356 -0.000501187 1.67953e-05 200022 89.2 -0.000512861 1.73852e-05 193237 88.1762 -0.000524807 1.79959e-05 186683 87.164 -0.000537032 1.86281e-05 180351 86.1635 -0.000549541 1.92826e-05 174234 85.1743 -0.000562341 1.996e-05 168325 84.1965 -0.00057544 2.06613e-05 162616 83.2298 -0.000588844 2.13872e-05 157101 82.2742 -0.00060256 2.21387e-05 151773 81.3295 -0.000616595 2.29167e-05 146626 80.3956 -0.000630957 2.3722e-05 141653 79.4723 -0.000645654 2.45556e-05 136849 78.5596 -0.000660693 2.54185e-05 132208 77.6573 -0.000676083 2.63119e-05 127725 76.7654 -0.000691831 2.72366e-05 123394 75.8836 -0.000707946 2.81939e-05 119210 75.0118 -0.000724436 2.91848e-05 115168 74.1501 -0.00074131 3.02107e-05 111263 73.2981 -0.000758578 3.12726e-05 107490 72.4559 -0.000776247 3.23719e-05 103845 71.6233 -0.000794328 3.35099e-05 100325 70.8002 -0.000812831 3.46879e-05 96923.1 69.9865 -0.000831764 3.59074e-05 93637.1 69.1821 -0.000851138 3.71698e-05 90462.6 68.3869 -0.000870964 3.84766e-05 87395.7 67.6007 -0.000891251 3.98294e-05 84433 66.8235 -0.000912011 4.12299e-05 81570.7 66.0552 -0.000933254 4.26796e-05 78805.6 65.2957 -0.000954993 4.41803e-05 76134.3 64.5448 -0.000977237 4.57339e-05 73553.6 63.8025 -0.001 4.73421e-05 71060.5 63.0687 -0.00102329 4.9007e-05 68651.9 62.3432 -0.00104713 5.07305e-05 66325.1 61.626 -0.00107152 5.25146e-05 64077.2 60.917 -0.00109648 5.43615e-05 61905.5 60.2161 -0.00112202 5.62735e-05 59807.5 59.5232 -0.00114815 5.82528e-05 57780.7 58.8381 -0.0011749 6.03017e-05 55822.7 58.1609 -0.00120226 6.24228e-05 53931 57.4914 -0.00123027 6.46186e-05 52103.6 56.8296 -0.00125893 6.68917e-05 50338.1 56.1753 -0.00128825 6.92449e-05 48632.5 55.5284 -0.00131826 7.16809e-05 46984.8 54.8889 -0.00134896 7.42027e-05 45392.9 54.2567 -0.00138038 7.68133e-05 43855 53.6318 -0.00141254 7.95159e-05 42369.3 53.0139 -0.00144544 8.23136e-05 40934 52.4031 -0.00147911 8.52099e-05 39547.4 51.7992 -0.00151356 8.82081e-05 38207.8 51.2023 -0.00154882 9.1312e-05 36913.7 50.6121 -0.00158489 9.45253e-05 35663.4 50.0286 -0.00162181 9.78517e-05 34455.5 49.4518 -0.00165959 0.000101295 33288.6 48.8816 -0.00169824 0.00010486 32161.3 48.3179 -0.0017378 0.000108551 31072.2 47.7606 -0.00177828 0.000112371 30020.1 47.2096 -0.0018197 0.000116327 29003.6 46.6649 -0.00186209 0.000120421 28021.5 46.1264 -0.00190546 0.00012466 27072.8 45.5941 -0.00194984 0.000129048 26156.3 45.0678 -0.00199526 0.000133591 25270.8 44.5474 -0.00204174 0.000138294 24415.4 44.0331 -0.0020893 0.000143163 23588.9 43.5245 -0.00213796 0.000148203 22790.5 43.0218 -0.00218776 0.000153422 22019.1 42.5247 -0.00223872 0.000158824 21273.9 42.0334 -0.00229087 0.000164416 20554 41.5476 -0.00234423 0.000170206 19858.5 41.0673 -0.00239883 0.0001762 19186.5 40.5925 -0.00245471 0.000182405 18537.3 40.1231 -0.00251189 0.000188829 17910.1 39.6591 -0.0025704 0.00019548 17304.2 39.2003 -0.00263027 0.000202365 16718.9 38.7467 -0.00269154 0.000209493 16153.3 38.2983 -0.00275423 0.000216873 15607 37.855 -0.00281838 0.000224512 15079.1 37.4167 -0.00288403 0.000232422 14569.2 36.9834 -0.00295121 0.00024061 14076.5 36.5551 -0.00301995 0.000249088 13600.5 36.1316 -0.0030903 0.000257864 13140.7 35.7129 -0.00316228 0.00026695 12696.4 35.2989 -0.00323594 0.000276357 12267.2 34.8897 -0.00331131 0.000286096 11852.5 34.4851 -0.00338844 0.000296179 11451.9 34.0851 -0.00346737 0.000306618 11064.9 33.6896 -0.00354813 0.000317425 10690.9 33.2986 -0.00363078 0.000328614 10329.7 32.9121 -0.00371535 0.000340198 9980.64 32.5299 -0.00380189 0.000352191 9643.44 32.1521 -0.00389045 0.000364608 9317.67 31.7786 -0.00398107 0.000377463 9002.92 31.4093 -0.0040738 0.000390773 8698.84 31.0442 -0.00416869 0.000404552 8405.06 30.6832 -0.0042658 0.000418819 8121.22 30.3263 -0.00436516 0.00043359 7847.01 29.9735 -0.00446684 0.000448883 7582.07 29.6247 -0.00457088 0.000464716 7326.11 29.2798 -0.00467735 0.000481109 7078.82 28.9388 -0.0047863 0.000498081 6839.9 28.6017 -0.00489779 0.000515653 6609.07 28.2684 -0.00501187 0.000533847 6386.05 27.9389 -0.00512861 0.000552684 6170.59 27.6131 -0.00524807 0.000572188 5962.42 27.291 -0.00537032 0.000592381 5761.29 26.9725 -0.00549541 0.000613289 5566.98 26.6577 -0.00562341 0.000634936 5379.24 26.3464 -0.0057544 0.00065735 5197.86 26.0386 -0.00588844 0.000680557 5022.61 25.7343 -0.0060256 0.000704585 4853.29 25.4334 -0.00616595 0.000729464 4689.71 25.136 -0.00630957 0.000755223 4531.66 24.8419 -0.00645654 0.000781895 4378.95 24.5511 -0.00660693 0.000809511 4231.42 24.2636 -0.00676083 0.000838106 4088.87 23.9794 -0.00691831 0.000867713 3951.15 23.6984 -0.00707946 0.000898369 3818.09 23.4205 -0.00724436 0.000930112 3689.52 23.1458 -0.0074131 0.00096298 3565.31 22.8742 -0.00758578 0.000997014 3445.29 22.6056 -0.00776247 0.00103225 3329.33 22.3401 -0.00794328 0.00106874 3217.3 22.0775 -0.00812831 0.00110652 3109.05 21.818 -0.00831764 0.00114565 3004.46 21.5613 -0.00851138 0.00118616 2903.41 21.3075 -0.00870964 0.00122811 2805.77 21.0567 -0.00891251 0.00127155 2711.44 20.8086 -0.00912011 0.00131652 2620.29 20.5633 -0.00933254 0.0013631 2532.22 20.3208 -0.00954993 0.00141133 2447.13 20.081 -0.00977237 0.00146127 2364.91 19.8439 -0.01 0.00151299 2285.47 19.6095 -0.0102329 0.00156654 2208.72 19.3777 -0.0104713 0.001622 2134.55 19.1485 -0.0107152 0.00167943 2062.9 18.9219 -0.0109648 0.0017389 1993.66 18.6978 -0.0112202 0.00180049 1926.76 18.4763 -0.0114815 0.00186426 1862.12 18.2572 -0.011749 0.00193031 1799.66 18.0406 -0.0120226 0.0019987 1739.31 17.8265 -0.0123027 0.00206953 1680.99 17.6147 -0.0125893 0.00214288 1624.65 17.4053 -0.0128825 0.00221884 1570.21 17.1983 -0.0131826 0.00229751 1517.6 16.9936 -0.0134896 0.00237898 1466.77 16.7911 -0.0138038 0.00246335 1417.65 16.591 -0.0141254 0.00255073 1370.19 16.3931 -0.0144544 0.00264123 1324.33 16.1974 -0.0147911 0.00273495 1280.02 16.0039 -0.0151356 0.00283202 1237.2 15.8125 -0.0154882 0.00293255 1195.83 15.6233 -0.0158489 0.00303666 1155.85 15.4362 -0.0162181 0.0031445 1117.21 15.2513 -0.0165959 0.00325618 1079.88 15.0683 -0.0169824 0.00337186 1043.81 14.8875 -0.017378 0.00349166 1008.95 14.7086 -0.0177828 0.00361575 975.269 14.5317 -0.018197 0.00374428 942.72 14.3568 -0.0186209 0.0038774 911.268 14.1839 -0.0190546 0.00401529 880.874 14.0129 -0.0194984 0.00415811 851.504 13.8438 -0.0199526 0.00430604 823.122 13.6766 -0.0204174 0.00445927 795.696 13.5112 -0.020893 0.00461799 769.192 13.3477 -0.0213796 0.00478239 743.58 13.186 -0.0218776 0.00495269 718.83 13.0261 -0.0223872 0.0051291 694.912 12.868 -0.0229087 0.00531184 671.799 12.7117 -0.0234423 0.00550113 649.462 12.5571 -0.0239883 0.00569722 627.877 12.4042 -0.0245471 0.00590035 607.017 12.253 -0.0251189 0.00611078 586.858 12.1034 -0.025704 0.00632877 567.376 11.9556 -0.0263027 0.0065546 548.549 11.8094 -0.0269153 0.00678855 530.353 11.6648 -0.0275423 0.00703092 512.769 11.5218 -0.0281838 0.00728202 495.775 11.3804 -0.0288403 0.00754216 479.352 11.2405 -0.0295121 0.00781167 463.479 11.1023 -0.0301995 0.0080909 448.139 10.9655 -0.030903 0.0083802 433.314 10.8303 -0.0316228 0.00867993 418.986 10.6965 -0.0323594 0.00899049 405.138 10.5642 -0.0331131 0.00931226 391.754 10.4334 -0.0338844 0.00964566 378.819 10.3041 -0.0346737 0.0099911 366.318 10.1762 -0.0354813 0.010349 354.235 10.0496 -0.0363078 0.0107199 342.556 9.92453 -0.0371535 0.0111043 331.269 9.80079 -0.0380189 0.0115025 320.359 9.67842 -0.0389045 0.0119152 309.815 9.5574 -0.0398107 0.0123428 299.623 9.43771 -0.040738 0.012786 289.772 9.31934 -0.0416869 0.0132452 280.25 9.20228 -0.042658 0.0137211 271.047 9.0865 -0.0436516 0.0142143 262.151 8.972 -0.0446684 0.0147255 253.553 8.85876 -0.0457088 0.0152552 245.242 8.74677 -0.0467735 0.0158043 237.208 8.63601 -0.047863 0.0163733 229.442 8.52646 -0.0489779 0.0169631 221.935 8.41812 -0.0501187 0.0175744 214.679 8.31097 -0.0512861 0.018208 207.665 8.20499 -0.0524807 0.0188647 200.885 8.10018 -0.0537032 0.0195455 194.33 7.99652 -0.0549541 0.0202511 187.994 7.894 -0.0562341 0.0209826 181.869 7.7926 -0.057544 0.0217408 175.948 7.69231 -0.0588844 0.0225268 170.224 7.59312 -0.060256 0.0233417 164.691 7.49502 -0.0616595 0.0241865 159.341 7.39799 -0.0630957 0.0250622 154.169 7.30203 -0.0645654 0.0259702 149.17 7.20712 -0.0660693 0.0269116 144.336 7.11324 -0.0676083 0.0278877 139.663 7.02039 -0.0691831 0.0288997 135.145 6.92856 -0.0707946 0.029949 130.777 6.83773 -0.0724436 0.031037 126.553 6.74789 -0.074131 0.0321652 122.47 6.65904 -0.0758578 0.0333351 118.522 6.57115 -0.0776247 0.0345483 114.705 6.48423 -0.0794328 0.0358064 111.015 6.39826 -0.081283 0.0371112 107.446 6.31322 -0.0831764 0.0384643 103.996 6.22912 -0.0851138 0.0398677 100.659 6.14593 -0.0870964 0.0413232 97.4334 6.06365 -0.0891251 0.0428329 94.314 5.98226 -0.0912011 0.0443987 91.2976 5.90177 -0.0933254 0.046023 88.3808 5.82215 -0.0954993 0.0477078 85.5603 5.74341 -0.0977237 0.0494555 82.8328 5.66552 -0.1 0.0512686 80.1952 5.58848 -0.102329 0.0531495 77.6445 5.51228 -0.104713 0.0551009 75.1778 5.43691 -0.107152 0.0571254 72.7923 5.36237 -0.109648 0.0592259 70.4853 5.28863 -0.112202 0.0614054 68.2541 5.2157 -0.114815 0.0636669 66.0963 5.14357 -0.11749 0.0660136 64.0093 5.07223 -0.120226 0.0684487 61.9909 5.00166 -0.123027 0.0709758 60.0387 4.93186 -0.125893 0.0735984 58.1504 4.86282 -0.128825 0.0763202 56.3241 4.79454 -0.131826 0.0791452 54.5576 4.72701 -0.134896 0.0820773 52.849 4.66021 -0.138038 0.0851208 51.1962 4.59414 -0.141254 0.0882801 49.5975 4.52879 -0.144544 0.0915597 48.051 4.46416 -0.147911 0.0949644 46.555 4.40023 -0.151356 0.098499 45.1078 4.33701 -0.154882 0.102169 43.7078 4.27447 -0.158489 0.105979 42.3535 4.21262 -0.162181 0.109936 41.0433 4.15145 -0.165959 0.114044 39.7757 4.09095 -0.169824 0.11831 38.5493 4.03111 -0.17378 0.122741 37.3628 3.97193 -0.177828 0.127343 36.2149 3.9134 -0.18197 0.132122 35.1042 3.85551 -0.186209 0.137086 34.0296 3.79826 -0.190546 0.142243 32.9898 3.74164 -0.194984 0.147599 31.9836 3.68565 -0.199526 0.153164 31.0101 3.63027 -0.204174 0.158946 30.068 3.5755 -0.20893 0.164953 29.1564 3.52134 -0.213796 0.171194 28.2742 3.46778 -0.218776 0.17768 27.4205 3.41481 -0.223872 0.18442 26.5943 3.36243 -0.229087 0.191425 25.7947 3.31063 -0.234423 0.198706 25.0209 3.25941 -0.239883 0.206273 24.2719 3.20876 -0.245471 0.214139 23.5471 3.15867 -0.251189 0.222316 22.8455 3.10914 -0.25704 0.230818 22.1664 3.06017 -0.263027 0.239657 21.5091 3.01175 -0.269153 0.248847 20.8729 2.96387 -0.275423 0.258404 20.257 2.91653 -0.281838 0.268343 19.6609 2.86972 -0.288403 0.27868 19.0837 2.82344 -0.295121 0.289431 18.525 2.77768 -0.301995 0.300614 17.9842 2.73245 -0.30903 0.312247 17.4605 2.68772 -0.316228 0.324351 16.9535 2.64351 -0.323594 0.336943 16.4627 2.59981 -0.331131 0.350047 15.9874 2.5566 -0.338844 0.363682 15.5273 2.51389 -0.346737 0.377874 15.0817 2.47167 -0.354813 0.392644 14.6503 2.42994 -0.363078 0.408019 14.2325 2.38869 -0.371535 0.424024 13.8279 2.34793 -0.380189 0.440687 13.4361 2.30763 -0.389045 0.458037 13.0567 2.26781 -0.398107 0.476103 12.6893 2.22846 -0.40738 0.494918 12.3334 2.18957 -0.416869 0.514515 11.9887 2.15114 -0.42658 0.534926 11.6549 2.11316 -0.436516 0.55619 11.3315 2.07564 -0.446684 0.578344 11.0183 2.03857 -0.457088 0.601427 10.715 2.00194 -0.467735 0.625481 10.4211 1.96575 -0.47863 0.650551 10.1364 1.93001 -0.489779 0.67668 9.86059 1.89469 -0.501187 0.703919 9.59341 1.85981 -0.512861 0.732316 9.33455 1.82536 -0.524807 0.761925 9.08377 1.79133 -0.537032 0.792801 8.84078 1.75773 -0.549541 0.825002 8.60535 1.72454 -0.562341 0.85859 8.37723 1.69177 -0.57544 0.893628 8.15619 1.65941 -0.588844 0.930183 7.942 1.62747 -0.60256 0.968327 7.73444 1.59593 -0.616595 1.00813 7.53329 1.56479 -0.630957 1.04968 7.33836 1.53406 -0.645654 1.09305 7.14945 1.50373 -0.660693 1.13833 6.96636 1.47379 -0.676083 1.1856 6.7889 1.44425 -0.691831 1.23497 6.61691 1.4151 -0.707946 1.28654 6.4502 1.38634 -0.724436 1.34041 6.28861 1.35796 -0.74131 1.39668 6.13197 1.32997 -0.758578 1.45548 5.98013 1.30236 -0.776247 1.51694 5.83293 1.27513 -0.794328 1.58117 5.69024 1.24828 -0.812831 1.64832 5.5519 1.2218 -0.831764 1.71852 5.41778 1.19569 -0.851138 1.79194 5.28775 1.16996 -0.870964 1.86873 5.16168 1.14459 -0.891251 1.94905 5.03944 1.11959 -0.912011 2.03308 4.92092 1.09495 -0.933254 2.12102 4.80599 1.07068 -0.954993 2.21305 4.69455 1.04676 -0.977237 2.30938 4.58648 1.0232 -1 2.41024 4.48169 1 -1.02329 2.51585 4.38006 0.977151 -1.04713 2.62646 4.2815 0.954653 -1.07152 2.74232 4.18592 0.932504 -1.09648 2.86371 4.09322 0.910703 -1.12202 2.99091 4.00331 0.889247 -1.14815 3.12423 3.91611 0.868135 -1.1749 3.26398 3.83153 0.847364 -1.20226 3.41051 3.74949 0.826932 -1.23027 3.56417 3.66991 0.806838 -1.25893 3.72534 3.59272 0.787079 -1.28825 3.89443 3.51784 0.767653 -1.31826 4.07185 3.44521 0.748559 -1.34896 4.25805 3.37475 0.729794 -1.38038 4.4535 3.30639 0.711356 -1.41254 4.65872 3.24008 0.693244 -1.44544 4.87422 3.17575 0.675454 -1.47911 5.10057 3.11333 0.657986 -1.51356 5.33837 3.05278 0.640836 -1.54882 5.58825 2.99403 0.624002 -1.58489 5.85087 2.93702 0.607483 -1.62181 6.12695 2.88172 0.591276 -1.65959 6.41724 2.82805 0.575379 -1.69824 6.72252 2.77598 0.55979 -1.7378 7.04365 2.72546 0.544505 -1.77828 7.38153 2.67643 0.529524 -1.8197 7.73709 2.62886 0.514842 -1.86209 8.11135 2.5827 0.500459 -1.90546 8.50538 2.5379 0.486371 -1.94984 8.92032 2.49443 0.472575 -1.99526 9.35736 2.45224 0.45907 -2.04174 9.81779 2.4113 0.445852 -2.0893 10.303 2.37157 0.432919 -2.13796 10.8143 2.33301 0.420267 -2.18776 11.3534 2.29558 0.407895 -2.23872 11.9219 2.25926 0.395799 -2.29087 12.5214 2.22401 0.383977 -2.34423 13.1539 2.1898 0.372424 -2.39883 13.8212 2.15659 0.361139 -2.45471 14.5255 2.12436 0.350118 -2.51189 15.269 2.09307 0.339358 -2.5704 16.054 2.0627 0.328856 -2.63027 16.883 2.03322 0.318608 -2.69153 17.7587 2.00461 0.308611 -2.75423 18.684 1.97683 0.298863 -2.81838 19.6618 1.94987 0.289359 -2.88403 20.6954 1.92369 0.280095 -2.95121 21.7881 1.89828 0.27107 -3.01995 22.9436 1.87361 0.262278 -3.0903 24.1659 1.84966 0.253717 -3.16228 25.4589 1.8264 0.245383 -3.23594 26.8271 1.80383 0.237272 -3.31131 28.2752 1.7819 0.229381 -3.38844 29.8082 1.76062 0.221706 -3.46737 31.4314 1.73995 0.214242 -3.54813 33.1505 1.71988 0.206987 -3.63078 34.9714 1.70039 0.199937 -3.71535 36.9006 1.68146 0.193087 -3.80189 38.9451 1.66309 0.186435 -3.89045 41.112 1.64524 0.179976 -3.98107 43.4092 1.6279 0.173706 -4.0738 45.845 1.61106 0.167621 -4.16869 48.4282 1.59471 0.161719 -4.2658 51.1684 1.57882 0.155994 -4.36516 54.0756 1.56339 0.150444 -4.46684 57.1606 1.5484 0.145064 -4.57088 60.4348 1.53384 0.13985 -4.67735 63.9105 1.51969 0.134799 -4.7863 67.6007 1.50595 0.129907 -4.89779 71.5195 1.49259 0.12517 -5.01187 75.6817 1.47961 0.120585 -5.12861 80.1031 1.46701 0.116148 -5.24808 84.8009 1.45475 0.111854 -5.37032 89.7929 1.44285 0.107702 -5.49541 95.0987 1.43127 0.103686 -5.62341 100.739 1.42003 0.0998034 -5.7544 106.735 1.40909 0.096051 -5.88844 113.112 1.39847 0.092425 -6.0256 119.893 1.38814 0.0889219 -6.16595 127.106 1.3781 0.0855386 -6.30957 134.78 1.36833 0.0822715 -6.45654 142.945 1.35884 0.0791175 -6.60693 151.634 1.34961 0.0760733 -6.76083 160.882 1.34064 0.0731357 -6.91831 170.726 1.33191 0.0703017 -7.07946 181.207 1.32343 0.0675681 -7.24436 192.366 1.31517 0.064932 -7.4131 204.249 1.30715 0.0623904 -7.58578 216.905 1.29934 0.0599404 -7.76247 230.386 1.29174 0.0575792 -7.94328 244.748 1.28435 0.0553041 -8.12831 260.051 1.27716 0.0531123 -8.31764 276.356 1.27017 0.0510011 -8.51138 293.734 1.26337 0.048968 -8.70964 312.255 1.25675 0.0470105 -8.91251 331.998 1.2503 0.045126 -9.12011 353.046 1.24403 0.0433122 -9.33254 375.488 1.23793 0.0415668 -9.54993 399.417 1.23199 0.0398873 -9.77237 424.937 1.22621 0.0382717 -10 452.154 1.22058 0.0367177 -10.2329 481.186 1.2151 0.0352233 -10.4713 512.157 1.20977 0.0337863 -10.7152 545.198 1.20457 0.0324047 -10.9648 580.452 1.19952 0.0310767 -11.2202 618.07 1.19459 0.0298004 -11.4815 658.216 1.1898 0.0285738 -11.749 701.062 1.18513 0.0273953 -12.0226 746.794 1.18058 0.0262631 -12.3027 795.611 1.17615 0.0251755 -12.5893 847.725 1.17184 0.0241309 -12.8825 903.364 1.16764 0.0231278 -13.1826 962.771 1.16354 0.0221646 -13.4896 1026.21 1.15955 0.0212398 -13.8038 1093.95 1.15567 0.0203521 -14.1254 1166.29 1.15188 0.0194999 -14.4544 1243.56 1.14819 0.0186821 -14.7911 1326.09 1.1446 0.0178973 -15.1356 1414.26 1.1411 0.0171443 -15.4882 1508.44 1.13768 0.0164218 -15.8489 1609.06 1.13436 0.0157286 -16.2181 1716.56 1.13112 0.0150638 -16.5959 1831.43 1.12796 0.0144261 -16.9824 1954.17 1.12488 0.0138146 -17.378 2085.33 1.12187 0.0132281 -17.7828 2225.51 1.11895 0.0126658 -18.197 2375.32 1.11609 0.0121267 -18.6209 2535.44 1.11331 0.0116098 -19.0546 2706.59 1.1106 0.0111144 -19.4984 2889.54 1.10795 0.0106395 -19.9526 3085.12 1.10538 0.0101843 -20.4174 3294.2 1.10286 0.00974812 -20.893 3517.73 1.10041 0.00933012 -21.3796 3756.72 1.09802 0.0089296 -21.8776 4012.25 1.09569 0.00854585 -22.3872 4285.48 1.09341 0.0081782 -22.9087 4577.65 1.09119 0.007826 -23.4423 4890.09 1.08903 0.00748862 -23.9883 5224.22 1.08692 0.00716546 -24.5471 5581.56 1.08486 0.00685595 -25.1189 5963.74 1.08285 0.00655952 -25.704 6372.5 1.08089 0.00627565 -26.3027 6809.72 1.07898 0.00600381 -26.9153 7277.38 1.07712 0.00574352 -27.5423 7777.65 1.0753 0.00549429 -28.1838 8312.79 1.07353 0.00525568 -28.8403 8885.27 1.07179 0.00502725 -29.5121 9497.71 1.0701 0.00480856 -30.1995 10152.9 1.06846 0.00459923 -30.903 10853.9 1.06685 0.00439885 -31.6228 11604 1.06528 0.00420706 -32.3594 12406.5 1.06375 0.0040235 -33.1131 13265.2 1.06225 0.00384782 -33.8844 14184 1.06079 0.00367969 -34.6737 15167.2 1.05937 0.0035188 -35.4813 16219.3 1.05798 0.00336484 -36.3078 17345.2 1.05662 0.00321752 -37.1535 18550.1 1.0553 0.00307656 -38.0189 19839.5 1.05401 0.00294169 -38.9045 21219.5 1.05275 0.00281266 -39.8107 22696.5 1.05152 0.00268921 -40.738 24277.2 1.05032 0.00257112 -41.6869 25969.1 1.04914 0.00245814 -42.658 27780.1 1.048 0.00235007 -43.6516 29718.4 1.04688 0.0022467 -44.6684 31793.2 1.04579 0.00214782 -45.7088 34014 1.04472 0.00205324 -46.7735 36391.3 1.04368 0.00196279 -47.863 38936.1 1.04267 0.00187627 -48.9779 41660.4 1.04168 0.00179353 -50.1187 44576.7 1.04071 0.0017144 -51.2861 47698.7 1.03977 0.00163873 -52.4808 51040.9 1.03884 0.00156636 -53.7032 54619.1 1.03794 0.00149716 -54.9541 58449.9 1.03706 0.00143099 -56.2341 62551.2 1.0362 0.00136772 -57.544 66942.3 1.03537 0.00130722 -58.8844 71643.6 1.03455 0.00124937 -60.256 76677.2 1.03375 0.00119406 -61.6595 82066.8 1.03297 0.00114118 -63.0957 87837.4 1.0322 0.00109063 -64.5654 94016.3 1.03146 0.00104229 -66.0694 100632 1.03073 0.000996083 -67.6083 107717 1.03002 0.000951907 -69.1831 115302 1.02933 0.000909677 -70.7946 123425 1.02865 0.000869306 -72.4436 132124 1.02799 0.000830715 -74.131 141438 1.02735 0.000793825 -75.8578 151412 1.02672 0.000758563 -77.6247 162093 1.0261 0.000724857 -79.4328 173532 1.0255 0.000692639 -81.283 185781 1.02491 0.000661845 -83.1764 198899 1.02434 0.000632411 -85.1138 212947 1.02378 0.000604278 -87.0964 227992 1.02323 0.00057739 -89.1251 244105 1.02269 0.000551691 -91.2011 261361 1.02217 0.00052713 -93.3254 279841 1.02166 0.000503657 -95.4993 299634 1.02116 0.000481223 -97.7237 320832 1.02068 0.000459784 -100 343535 1.0202 0.000439295 -102.329 367850 1.01974 0.000419714 -104.713 393894 1.01928 0.000401003 -107.152 421787 1.01884 0.000383121 -109.648 451662 1.01841 0.000366033 -112.202 483661 1.01798 0.000349704 -114.815 517934 1.01757 0.0003341 -117.49 554644 1.01717 0.00031919 -120.226 593963 1.01677 0.000304942 -123.027 636079 1.01639 0.000291328 -125.893 681189 1.01601 0.000278319 -128.825 729509 1.01565 0.000265889 -131.826 781265 1.01529 0.000254011 -134.896 836705 1.01494 0.000242662 -138.038 896089 1.01459 0.000231819 -141.254 959698 1.01426 0.000221459 -144.544 1.02784e+06 1.01393 0.000211559 -147.911 1.10082e+06 1.01361 0.000202101 -151.356 1.179e+06 1.0133 0.000193065 -154.882 1.26275e+06 1.013 0.000184431 -158.489 1.35246e+06 1.0127 0.000176182 -162.181 1.44856e+06 1.01241 0.000168301 -165.959 1.5515e+06 1.01212 0.000160771 -169.824 1.66178e+06 1.01185 0.000153577 -173.78 1.77991e+06 1.01157 0.000146705 -177.828 1.90645e+06 1.01131 0.000140139 -181.97 2.04201e+06 1.01105 0.000133866 -186.209 2.18723e+06 1.0108 0.000127873 -190.546 2.3428e+06 1.01055 0.000122147 -194.984 2.50945e+06 1.01031 0.000116678 -199.526 2.68797e+06 1.01007 0.000111452 -204.174 2.87922e+06 1.00984 0.00010646 -208.93 3.08411e+06 1.00962 0.000101692 -213.796 3.3036e+06 1.0094 9.71358e-05 -218.776 3.53873e+06 1.00918 9.27837e-05 -223.872 3.79063e+06 1.00897 8.86262e-05 -229.087 4.06049e+06 1.00877 8.46546e-05 -234.423 4.34959e+06 1.00857 8.08606e-05 -239.883 4.65931e+06 1.00837 7.72362e-05 -245.471 4.99111e+06 1.00818 7.3774e-05 -251.189 5.34658e+06 1.00799 7.04667e-05 -257.04 5.7274e+06 1.00781 6.73073e-05 -263.027 6.13539e+06 1.00763 6.42894e-05 -269.154 6.57247e+06 1.00746 6.14065e-05 -275.423 7.04074e+06 1.00729 5.86527e-05 -281.838 7.54241e+06 1.00712 5.60222e-05 -288.403 8.07987e+06 1.00696 5.35094e-05 -295.121 8.65568e+06 1.0068 5.11092e-05 -301.995 9.27258e+06 1.00664 4.88164e-05 -309.03 9.93349e+06 1.00649 4.66263e-05 -316.228 1.06416e+07 1.00634 4.45344e-05 -323.594 1.14002e+07 1.0062 4.25361e-05 -331.131 1.22129e+07 1.00606 4.06274e-05 -338.844 1.30837e+07 1.00592 3.88042e-05 -346.737 1.40166e+07 1.00578 3.70627e-05 -354.813 1.50161e+07 1.00565 3.53993e-05 -363.078 1.60869e+07 1.00552 3.38104e-05 -371.535 1.72342e+07 1.0054 3.22927e-05 -380.189 1.84633e+07 1.00527 3.08431e-05 -389.045 1.97803e+07 1.00515 2.94584e-05 -398.107 2.11912e+07 1.00504 2.81359e-05 -407.38 2.27029e+07 1.00492 2.68726e-05 -416.869 2.43225e+07 1.00481 2.5666e-05 -426.58 2.60578e+07 1.0047 2.45135e-05 -436.516 2.79169e+07 1.00459 2.34127e-05 -446.684 2.99088e+07 1.00449 2.23613e-05 -457.088 3.2043e+07 1.00438 2.13571e-05 -467.735 3.43295e+07 1.00429 2.03979e-05 -478.63 3.67793e+07 1.00419 1.94817e-05 -489.779 3.94041e+07 1.00409 1.86067e-05 -501.187 4.22164e+07 1.004 1.77709e-05 -512.861 4.52295e+07 1.00391 1.69726e-05 -524.807 4.84578e+07 1.00382 1.62101e-05 -537.032 5.19167e+07 1.00373 1.54819e-05 -549.541 5.56227e+07 1.00365 1.47864e-05 -562.341 5.95934e+07 1.00356 1.4122e-05 -575.44 6.38477e+07 1.00348 1.34875e-05 -588.844 6.84059e+07 1.0034 1.28815e-05 -602.56 7.32897e+07 1.00332 1.23027e-05 -616.595 7.85224e+07 1.00325 1.17499e-05 -630.957 8.4129e+07 1.00317 1.12219e-05 -645.654 9.0136e+07 1.0031 1.07176e-05 -660.693 9.65723e+07 1.00303 1.02359e-05 -676.083 1.03468e+08 1.00296 9.77592e-06 -691.831 1.10857e+08 1.00289 9.33656e-06 -707.946 1.18774e+08 1.00283 8.91693e-06 -724.436 1.27256e+08 1.00276 8.51615e-06 -741.31 1.36344e+08 1.0027 8.13337e-06 -758.578 1.46082e+08 1.00264 7.76778e-06 -776.247 1.56516e+08 1.00258 7.41861e-06 -794.328 1.67695e+08 1.00252 7.08514e-06 -812.831 1.79673e+08 1.00246 6.76664e-06 -831.764 1.92507e+08 1.00241 6.46245e-06 -851.138 2.06258e+08 1.00235 6.17193e-06 -870.964 2.20992e+08 1.0023 5.89446e-06 -891.251 2.36778e+08 1.00225 5.62946e-06 -912.011 2.53693e+08 1.0022 5.37637e-06 -933.254 2.71816e+08 1.00215 5.13465e-06 -954.993 2.91235e+08 1.0021 4.9038e-06 -977.237 3.12042e+08 1.00205 4.68331e-06 -1000 3.34336e+08 1.002 4.47274e-06 -1023.29 3.58222e+08 1.00196 4.27163e-06 -1047.13 3.83817e+08 1.00191 4.07955e-06 -1071.52 4.1124e+08 1.00187 3.89611e-06 -1096.48 4.40624e+08 1.00182 3.72092e-06 -1122.02 4.72107e+08 1.00178 3.5536e-06 -1148.15 5.05841e+08 1.00174 3.3938e-06 -1174.9 5.41986e+08 1.0017 3.24118e-06 -1202.26 5.80715e+08 1.00166 3.09542e-06 -1230.27 6.22212e+08 1.00163 2.95622e-06 -1258.93 6.66675e+08 1.00159 2.82327e-06 -1288.25 7.14316e+08 1.00155 2.6963e-06 -1318.26 7.65363e+08 1.00152 2.57504e-06 -1348.96 8.2006e+08 1.00148 2.45922e-06 -1380.38 8.78665e+08 1.00145 2.34862e-06 -1412.54 9.4146e+08 1.00142 2.24299e-06 -1445.44 1.00874e+09 1.00138 2.14211e-06 -1479.11 1.08084e+09 1.00135 2.04576e-06 -1513.56 1.15808e+09 1.00132 1.95374e-06 -1548.82 1.24085e+09 1.00129 1.86587e-06 -1584.89 1.32954e+09 1.00126 1.78194e-06 -1621.81 1.42457e+09 1.00123 1.70179e-06 -1659.59 1.52639e+09 1.00121 1.62524e-06 -1698.24 1.63548e+09 1.00118 1.55214e-06 -1737.8 1.75238e+09 1.00115 1.48232e-06 -1778.28 1.87764e+09 1.00112 1.41564e-06 -1819.7 2.01185e+09 1.0011 1.35196e-06 -1862.09 2.15565e+09 1.00107 1.29115e-06 -1905.46 2.30974e+09 1.00105 1.23307e-06 -1949.84 2.47484e+09 1.00103 1.1776e-06 -1995.26 2.65175e+09 1.001 1.12462e-06 -2041.74 2.8413e+09 1.00098 1.07403e-06 -2089.3 3.04441e+09 1.00096 1.02571e-06 -2137.96 3.26203e+09 1.00094 9.79571e-07 diff --git a/MakeCloud.py b/MakeCloud.py index 15c9653..f872c41 100755 --- a/MakeCloud.py +++ b/MakeCloud.py @@ -1,10 +1,10 @@ #!/usr/bin/env python -""" +""" MakeCloud: "Believe me, we've got some very turbulent clouds, the best clouds. You're gonna love it." Usage: MakeCloud.py [options] -Options: +Options: -h --help Show this screen. --R= Outer radius of the cloud in pc [default: 10.0] --M= Mass of the cloud in msun [default: 2e4] @@ -27,7 +27,7 @@ --x_star= Position of the star, defaults to center of the box --star_stage= Evolutionary stage of the star/black hole [default: 7] --derefinement Apply radial derefinement to ambient cells outside of 3* cloud radius - --no_diffuse_gas Remove diffuse ISM envelope fills the rest of the box with uniform density. + --no_diffuse_gas Remove diffuse ISM envelope fills the rest of the box with uniform density. --phimode= Relative amplitude of m=2 density perturbation (e.g. for Boss-Bodenheimer test) [default: 0.0] --localdir Changes directory defaults assuming all files are used from local directory. --B_unit= Unit of magnetic field in gauss [default: 1e4] @@ -49,7 +49,6 @@ --Z= Metallicity of the cloud in Solar units (just for params file) [default: 1.0] --ISRF= Interstellar radiation background of the cloud in Solar neighborhood units (just for params file) [default: 1.0] """ -# Example: python MakeCloud.py --M=1000 --N=1e7 --R=1.0 --localdir --param_only import os import numpy as np @@ -59,8 +58,11 @@ from docopt import docopt -def get_glass_coords(N_gas, glass_path): - x = np.load(glass_path) +def get_glass_coords(N_gas, glass_path, center_on_cell=False): + x = h5py.File(glass_path)["Coordinates"][:] + if not center_on_cell: # if we don't want a cell at the exact box center + np.random.seed(42) + x = (x + np.random.rand(x.shape[1])) % 1.0 Nx = len(x) while len(x) * np.pi * 4 / 3 / 8 < N_gas: @@ -70,18 +72,15 @@ def get_glass_coords(N_gas, glass_path): ) x = np.concatenate( [ - x / 2 - + i * np.array([0.5, 0, 0]) - + j * np.array([0, 0.5, 0]) - + k * np.array([0, 0, 0.5]) + x / 2 + i * np.array([0.5, 0, 0]) + j * np.array([0, 0.5, 0]) + k * np.array([0, 0, 0.5]) for i in range(2) for j in range(2) for k in range(2) ] ) - Nx = len(x) print("Glass loaded!") - return x + order = x.max(axis=1).argsort() + return x[order] def TurbField(res=256, minmode=2, maxmode=64, slope=2.0, sol_weight=1.0, seed=42): @@ -115,25 +114,13 @@ def TurbField(res=256, minmode=2, maxmode=64, slope=2.0, sol_weight=1.0, seed=42 for j in range(3): if i == j: vk_new[i] += sol_weight * VK[j] - vk_new[i] += ( - (1 - 2 * sol_weight) - * freq3d[i] - * freq3d[j] - / (kSqr + 1e-300) - * VK[j] - ) + vk_new[i] += (1 - 2 * sol_weight) * freq3d[i] * freq3d[j] / (kSqr + 1e-300) * VK[j] vk_new[:, kSqr == 0] = 0.0 VK = vk_new - vel = np.array( - [fftpack.ifftn(vk).real for vk in VK] - ) # transform back to real space - vel -= np.average(vel, axis=(1, 2, 3))[ - :, np.newaxis, np.newaxis, np.newaxis - ] - vel = vel / np.sqrt( - np.sum(vel**2, axis=0).mean() - ) # normalize so that RMS is 1 + vel = np.array([fftpack.ifftn(vk).real for vk in VK]) # transform back to real space + vel -= np.average(vel, axis=(1, 2, 3))[:, np.newaxis, np.newaxis, np.newaxis] + vel = vel / np.sqrt(np.sum(vel**2, axis=0).mean()) # normalize so that RMS is 1 return np.array(vel) @@ -156,7 +143,6 @@ def TurbField(res=256, minmode=2, maxmode=64, slope=2.0, sol_weight=1.0, seed=42 minmode = int(arguments["--minmode"]) filename = arguments["--filename"] diffuse_gas = not arguments["--no_diffuse_gas"] -localdir = arguments["--localdir"] param_only = arguments["--param_only"] B_unit = float(arguments["--B_unit"]) length_unit = float(arguments["--length_unit"]) @@ -182,21 +168,19 @@ def TurbField(res=256, minmode=2, maxmode=64, slope=2.0, sol_weight=1.0, seed=42 if arguments["--glass_path"]: glass_path = arguments["--glass_path"] else: - glass_path = os.path.expanduser("~") + "/glass_orig.npy" + prefix = os.path.expanduser("~") + "/.makecloud_glass" + glass_path = prefix + "/glass_256.hdf5" if not os.path.exists(glass_path): + if not os.path.isdir(prefix): + os.mkdir(prefix) import urllib.request print("Downloading glass file...") urllib.request.urlretrieve( - "http://www.tapir.caltech.edu/~mgrudich/glass_orig.npy", + "https://users.flatironinstitute.org/~mgrudic/glass/glass_256.hdf5", glass_path, - # "https://data.obs.carnegiescience.edu/starforge/glass_orig.npy", glass_path ) -if localdir: - turb_path = "turb" - glass_path = "glass_256.npy" - if arguments["--boxsize"] is not None: boxsize = float(arguments["--boxsize"]) else: @@ -230,24 +214,17 @@ def TurbField(res=256, minmode=2, maxmode=64, slope=2.0, sol_weight=1.0, seed=42 turb_sol, ) + ("_%d" % seed) - + ( - "_collision_%g_%g_%g_%s" - % (impact_dist, impact_param, v_impact, impact_axis) - if impact_dist > 0 - else "" - ) + + ("_collision_%g_%g_%g_%s" % (impact_dist, impact_param, v_impact, impact_axis) if impact_dist > 0 else "") + ".hdf5" ) filename = filename.replace("+", "").replace("e0", "e") filename = "".join(filename.split()) -delta_m = M_gas / N_gas -delta_m_solar = delta_m / mass_unit +dm = M_gas / N_gas +dm_solar = dm / mass_unit rho_avg = 3 * M_gas / R**3 / (4 * np.pi) -if delta_m_solar < 0.1: # if we're doing something marginally IMF-resolving - softening = ( - 3.11e-5 # ~6.5 AU, minimum sink radius is 2.8 times that (~18 AU) - ) +if dm_solar < 0.1: # if we're doing something marginally IMF-resolving + softening = 3.11e-5 # ~6.5 AU, minimum sink radius is 2.8 times that (~18 AU) ncrit = 1e13 # ~100x the opacity limit else: # something more FIRE-like, where we rely on a sub-grid prescription turning gas into star particles softening = 0.1 @@ -269,13 +246,9 @@ def TurbField(res=256, minmode=2, maxmode=64, slope=2.0, sol_weight=1.0, seed=42 0.019111097819633344 * vrms**3 / L ) # ST_Energy sets the dissipation rate of SPECIFIC energy ~ v^2 / (L/v) ~ v^3/L -paramsfile = str( - open( - os.path.realpath(__file__).replace("MakeCloud.py", "params.txt"), "r" - ).read() -) +paramsfile = str(open(os.path.realpath(__file__).replace("MakeCloud.py", "params.txt"), "r").read()) -jet_particle_mass = min(delta_m, max(1e-4, delta_m / 10.0)) +jet_particle_mass = min(dm, max(1e-4, dm / 10.0)) MS_wind_particle_mass = ( jet_particle_mass / 10 ) # MS winds have lower mdot than jets, so we should be able to better resolve them this way @@ -292,7 +265,7 @@ def TurbField(res=256, minmode=2, maxmode=64, slope=2.0, sol_weight=1.0, seed=42 "OUTFOLDER": "output", "JET_PART_MASS": jet_particle_mass, "MS_WIND_PART_MASS": MS_wind_particle_mass, - "BH_SEED_MASS": delta_m / 2.0, + "BH_SEED_MASS": dm / 2.0, "TURBDECAY": tcross / 2, "TURBENERGY": turbenergy, "TURBFREQ": tcross / 20, @@ -311,9 +284,7 @@ def TurbField(res=256, minmode=2, maxmode=64, slope=2.0, sol_weight=1.0, seed=42 } for k, r in replacements.items(): - paramsfile = paramsfile.replace( - k, (r if isinstance(r, str) else "{:.2e}".format(r)) - ) + paramsfile = paramsfile.replace(k, (r if isinstance(r, str) else "{:.2e}".format(r))) open("params_" + filename.replace(".hdf5", "") + ".txt", "w").write(paramsfile) if makebox: @@ -330,27 +301,18 @@ def TurbField(res=256, minmode=2, maxmode=64, slope=2.0, sol_weight=1.0, seed=42 ) for k in replacements_box.keys(): paramsfile = paramsfile.replace(k, str(replacements_box[k])) - open("params_" + filename.replace(".hdf5", "") + "_BOX.txt", "w").write( - paramsfile - ) + open("params_" + filename.replace(".hdf5", "") + "_BOX.txt", "w").write(paramsfile) if makecylinder: # Get cylinder params - R_cyl = R * np.sqrt( - np.pi / (4 * cyl_aspect_ratio) - ) # surface density equivalent cylinder + R_cyl = R * np.sqrt(np.pi / (4 * cyl_aspect_ratio)) # surface density equivalent cylinder L_cyl = R_cyl * 2 * cyl_aspect_ratio vrms_cyl = ( - 2 * G * M_gas / L_cyl - ) ** 0.5 * turbulence**0.5 # the potential is different for a cylinder than for a sphere, so we need to rescale vrms to get the right alpha, using E_grav_cyl = -GM**2/L + (2 * G * M_gas / L_cyl) ** 0.5 * turbulence** 0.5 + ) # the potential is different for a cylinder than for a sphere, so we need to rescale vrms to get the right alpha, using E_grav_cyl = -GM**2/L vrms_cyl *= 0.71 # additional scaling found numerically to make the stirring run reproduce the right alpha and filament length (similarly determined numerical factor added to GIZMO) tcross_cyl = 2 * R_cyl / vrms_cyl - boxsize_cyl = ( - L_cyl * 1.5 + R_cyl * 5 - ) # the box should fit the cylinder and be many times bigger than its width - print( - "Cylinder params: L=%g R=%g boxsize=%g vrms=%g" - % (L_cyl, R_cyl, boxsize_cyl, vrms_cyl) - ) + boxsize_cyl = L_cyl * 1.5 + R_cyl * 5 # the box should fit the cylinder and be many times bigger than its width + print("Cylinder params: L=%g R=%g boxsize=%g vrms=%g" % (L_cyl, R_cyl, boxsize_cyl, vrms_cyl)) replacements_cyl = replacements.copy() replacements_cyl["NAME"] = filename.replace(".hdf5", "_CYL") replacements_cyl["BOXSIZE"] = boxsize_cyl @@ -373,15 +335,12 @@ def TurbField(res=256, minmode=2, maxmode=64, slope=2.0, sol_weight=1.0, seed=42 ) for k in replacements_cyl.keys(): paramsfile = paramsfile.replace(k, str(replacements_cyl[k])) - open("params_" + filename.replace(".hdf5", "") + "_CYL.txt", "w").write( - paramsfile - ) + open("params_" + filename.replace(".hdf5", "") + "_CYL.txt", "w").write(paramsfile) if param_only: print("Parameters only run, exiting...") exit() -dm = M_gas / N_gas mgas = np.repeat(dm, N_gas) x = get_glass_coords(N_gas, glass_path) @@ -395,9 +354,15 @@ def TurbField(res=256, minmode=2, maxmode=64, slope=2.0, sol_weight=1.0, seed=42 x *= (float(Nx) / N_gas * 4 * np.pi / 3 / 8) ** (1.0 / 3) * R print("Done! Recomupting radii...") r = cdist(x, [np.zeros(3)])[:, 0] +if np.any(r == 0): + raise ValueError( + "found point with r=0 in the glass file, we don't handle this case throughout our calculations yet. Stopping." + ) x, r = x / r.max(), r / r.max() print("Doing density profile...") rnew = r ** (3.0 / (3 + density_exponent)) * R + + x = x * (rnew / r)[:, None] r = np.sum(x**2, axis=1) ** 0.5 r_order = r.argsort() @@ -426,7 +391,7 @@ def TurbField(res=256, minmode=2, maxmode=64, slope=2.0, sol_weight=1.0, seed=42 Mr = mgas.cumsum() ugrav = G * np.sum(Mr / r * mgas) v -= np.average(v, axis=0) -Eturb = 0.5 * M_gas / N_gas * np.sum(v**2) +Eturb = 0.5 * dm * np.sum(v**2) v *= np.sqrt(turbulence * ugrav / Eturb) E_rot_target = spin * ugrav Rcyl = np.sqrt(x[:, 0] ** 2 + x[:, 1] ** 2) @@ -438,22 +403,13 @@ def TurbField(res=256, minmode=2, maxmode=64, slope=2.0, sol_weight=1.0, seed=42 B = np.c_[np.zeros(N_gas), np.zeros(N_gas), np.ones(N_gas)] vA_unit = ( - 3.429e8 - * B_unit - * (M_gas) ** -0.5 - * R**1.5 - * np.sqrt(4 * np.pi / 3) - / v_unit + 3.429e8 * B_unit * (M_gas) ** -0.5 * R**1.5 * np.sqrt(4 * np.pi / 3) / v_unit ) # alfven speed for unit magnetic field -uB = ( - 0.5 * M_gas * vA_unit**2 -) # magnetic energy we would have for unit magnetic field +uB = 0.5 * M_gas * vA_unit**2 # magnetic energy we would have for unit magnetic field if bfixed > 0: B = B * bfixed else: - B = B * np.sqrt( - magnetic_field * ugrav / uB - ) # renormalize to desired magnetic energy + B = B * np.sqrt(magnetic_field * ugrav / uB) # renormalize to desired magnetic energy v = v - np.average(v, axis=0) x = x - np.average(x, axis=0) @@ -461,25 +417,16 @@ def TurbField(res=256, minmode=2, maxmode=64, slope=2.0, sol_weight=1.0, seed=42 r, phi = np.sum(x**2, axis=1) ** 0.5, np.arctan2(x[:, 1], x[:, 0]) theta = np.arccos(x[:, 2] / r) phi += phimode * np.sin(2 * phi) / 2 -x = ( - r[:, np.newaxis] - * np.c_[ - np.cos(phi) * np.sin(theta), np.sin(phi) * np.sin(theta), np.cos(theta) - ] -) +x = r[:, np.newaxis] * np.c_[np.cos(phi) * np.sin(theta), np.sin(phi) * np.sin(theta), np.cos(theta)] if makecylinder: def ind_in_cylinder(x, L_cyl, R_cyl): - return (np.abs(x[:, 0]) < L_cyl / 2) & ( - np.sum(x[:, 1:] ** 2, axis=1) < R_cyl**2 - ) + return (np.abs(x[:, 0]) < L_cyl / 2) & (np.sum(x[:, 1:] ** 2, axis=1) < R_cyl**2) # Just get a roughly homogeneous cylinder along the x axis, we will stir it anyway N_cyl = 0 - while ( - N_cyl <= N_gas - ): # should be very unlikely that we need to repeat, but let's check to be sure + while N_cyl <= N_gas: # should be very unlikely that we need to repeat, but let's check to be sure x_cyl = np.random.rand(2 * N_gas, 3) * 2 - 1 x_cyl[:, 0] *= L_cyl / 2 x_cyl[:, 1] *= R_cyl @@ -526,23 +473,17 @@ def ind_in_cylinder(x, L_cyl, R_cyl): if diffuse_gas: # assuming 10K vs 10^4K gas: factor of ~10^3 density contrast rho_warm = M_gas * 3 / (4 * np.pi * R**3) / 1000 - M_warm = ( - boxsize**3 - (4 * np.pi * R**3 / 3) - ) * rho_warm # mass of diffuse box-filling medium - N_warm = int(M_warm / (M_gas / N_gas)) if derefinement: + M_warm = (boxsize**3 - (4 * np.pi * R**3 / 3)) * rho_warm # mass of diffuse box-filling medium + N_warm = int(M_warm / (dm)) x0 = get_glass_coords(N_gas, glass_path) Nx = len(x0) x0 = 2 * (x0 - 0.5) r0 = (x0 * x0).sum(1) ** 0.5 x0, r0 = x0[r0.argsort()], r0[r0.argsort()] # first lay down the stuff within 3*R - N_warm = int( - 4 * np.pi * rho_warm * (3 * R) ** 3 / 3 / dm - ) # number of cells within 3R - x_warm = ( - x0[:N_warm] * 3 * R / r0[N_warm - 1] - ) # uniform density of cells within 3R + N_warm = int(4 * np.pi * rho_warm * (3 * R) ** 3 / 3 / dm) # number of cells within 3R + x_warm = x0[:N_warm] * 3 * R / r0[N_warm - 1] # uniform density of cells within 3R x0 = x0[ N_warm: ] # now we take the ones outside the initial sphere and map them to a n(R) ~ R^-3 profile so that we get constant number of cells per log radius interval @@ -552,39 +493,31 @@ def ind_in_cylinder(x, L_cyl, R_cyl): x_warm = x_warm[np.max(np.abs(x_warm), axis=1) < boxsize / 2] N_warm = len(x_warm) R_warm = (x_warm * x_warm).sum(1) ** 0.5 - mgas = np.concatenate( - [mgas, np.clip(dm * (R_warm / (3 * R)) ** 3, dm, np.inf)] - ) + mgas = np.concatenate([mgas, np.clip(dm * (R_warm / (3 * R)) ** 3, dm, np.inf)]) else: - x_warm = boxsize * np.random.rand(N_warm, 3) - boxsize / 2 + M_warm = boxsize**3 * rho_warm # mass of diffuse box-filling medium + N_warm = int(M_warm / dm) # get glass with N_warm particles + x_warm = get_glass_coords(N_warm, glass_path)[:N_warm] + x_warm /= x_warm.max() + x_warm = boxsize * x_warm - boxsize / 2 if impact_dist == 0: x_warm = x_warm[np.sum(x_warm**2, axis=1) > R**2] N_warm = len(x_warm) - mgas = np.concatenate( - [mgas, np.repeat(mgas.sum() / len(mgas), N_warm)] - ) + mgas = np.concatenate([mgas, np.repeat(mgas.sum() / len(mgas), N_warm)]) x = np.concatenate([x, x_warm]) v = np.concatenate([v, np.zeros((N_warm, 3))]) Bmag = np.average(np.sum(B**2, axis=1)) ** 0.5 - B = np.concatenate( - [B, np.repeat(Bmag, N_warm)[:, np.newaxis] * np.array([0, 0, 1])] - ) + B = np.concatenate([B, np.repeat(Bmag, N_warm)[:, np.newaxis] * np.array([0, 0, 1])]) u = np.concatenate([u, np.repeat(101.0, N_warm)]) if makecylinder: # The magnetic field is paralell to the cylinder (true at low densities, so probably fine for IC) - B_cyl = np.concatenate( - [B, np.repeat(Bmag, N_warm)[:, np.newaxis] * np.array([1, 0, 0])] - ) + B_cyl = np.concatenate([B, np.repeat(Bmag, N_warm)[:, np.newaxis] * np.array([1, 0, 0])]) # Add diffuse medium M_warm_cyl = (boxsize_cyl**3 - (4 * np.pi * R**3 / 3)) * rho_warm - N_warm_cyl = int(M_warm_cyl / (M_gas / N_gas)) - x_warm = ( - boxsize_cyl * np.random.rand(N_warm_cyl, 3) - boxsize_cyl / 2 - ) # will be recentered later - x_warm = x_warm[ - ~ind_in_cylinder(x_warm, L_cyl, R_cyl) - ] # keep only warm gas outside the cylinder + N_warm_cyl = int(M_warm_cyl / (dm)) + x_warm = boxsize_cyl * np.random.rand(N_warm_cyl, 3) - boxsize_cyl / 2 # will be recentered later + x_warm = x_warm[~ind_in_cylinder(x_warm, L_cyl, R_cyl)] # keep only warm gas outside the cylinder # print("N_warm_cyl: %g N_warm_cyl_kept %g "%(N_warm_cyl,len(x_warm))) N_warm_cyl = len(x_warm) x_cyl = np.concatenate([x_cyl, x_warm]) @@ -635,34 +568,18 @@ def ind_in_cylinder(x, L_cyl, R_cyl): F.create_group("PartType5") # Let's add the sink at the center F["PartType5"].create_dataset("Masses", data=np.array([M_star])) - F["PartType5"].create_dataset( - "Coordinates", data=[x_star] - ) # at the center + F["PartType5"].create_dataset("Coordinates", data=[x_star]) # at the center F["PartType5"].create_dataset("Velocities", data=[v_star]) # at rest - F["PartType5"].create_dataset( - "ParticleIDs", data=np.array([F["PartType0/ParticleIDs"][:].max() + 1]) - ) + F["PartType5"].create_dataset("ParticleIDs", data=np.array([F["PartType0/ParticleIDs"][:].max() + 1])) # Advanced properties for sinks - F["PartType5"].create_dataset( - "BH_Mass", data=M_star - ) # all the mass in the sink/protostar/star - F["PartType5"].create_dataset( - "BH_Mass_AlphaDisk", data=np.array([0.0]) - ) # starts with no disk - F["PartType5"].create_dataset( - "BH_Mdot", data=np.array([0.0]) - ) # starts with no mdot - F["PartType5"].create_dataset( - "BH_Specific_AngMom", data=np.array([0.0]) - ) # starts with no angular momentum - F["PartType5"].create_dataset( - "SinkRadius", data=np.array([softening]) - ) # Sinkradius set to softening + F["PartType5"].create_dataset("BH_Mass", data=M_star) # all the mass in the sink/protostar/star + F["PartType5"].create_dataset("BH_Mass_AlphaDisk", data=np.array([0.0])) # starts with no disk + F["PartType5"].create_dataset("BH_Mdot", data=np.array([0.0])) # starts with no mdot + F["PartType5"].create_dataset("BH_Specific_AngMom", data=np.array([0.0])) # starts with no angular momentum + F["PartType5"].create_dataset("SinkRadius", data=np.array([softening])) # Sinkradius set to softening F["PartType5"].create_dataset("StellarFormationTime", data=np.array([0.0])) F["PartType5"].create_dataset("ProtoStellarAge", data=np.array([0.0])) - F["PartType5"].create_dataset( - "ProtoStellarStage", data=np.array([5], dtype=np.int32), dtype=np.int32 - ) + F["PartType5"].create_dataset("ProtoStellarStage", data=np.array([5], dtype=np.int32), dtype=np.int32) # Stellar properties # if (central_star or central_SN): # if central_star: @@ -676,12 +593,8 @@ def ind_in_cylinder(x, L_cyl, R_cyl): R_ZAMS = (M_star) ** 0.57 else: R_ZAMS = (M_star) ** 0.8 - F["PartType5"].create_dataset( - "ProtoStellarRadius_inSolar", data=np.array([R_ZAMS]) - ) # Sinkradius set to softening - F["PartType5"].create_dataset( - "StarLuminosity_Solar", data=np.array([0.0]) - ) # dummy + F["PartType5"].create_dataset("ProtoStellarRadius_inSolar", data=np.array([R_ZAMS])) # Sinkradius set to softening + F["PartType5"].create_dataset("StarLuminosity_Solar", data=np.array([0.0])) # dummy F["PartType5"].create_dataset("Mass_D", data=np.array([0.0])) # No D left if magnetic_field > 0.0: @@ -715,15 +628,13 @@ def ind_in_cylinder(x, L_cyl, R_cyl): F.create_group("Header") F["Header"].attrs["NumPart_ThisFile"] = [N_gas + N_warm_cyl, 0, 0, 0, 0, 0] F["Header"].attrs["NumPart_Total"] = [N_gas + N_warm_cyl, 0, 0, 0, 0, 0] - F["Header"].attrs["MassTable"] = [M_gas / N_gas, 0, 0, 0, 0, 0] + F["Header"].attrs["MassTable"] = [dm, 0, 0, 0, 0, 0] F["Header"].attrs["BoxSize"] = boxsize_cyl F["Header"].attrs["Time"] = 0.0 F["PartType0"].create_dataset("Masses", data=mgas) F["PartType0"].create_dataset("Coordinates", data=x_cyl) F["PartType0"].create_dataset("Velocities", data=v_cyl) - F["PartType0"].create_dataset( - "ParticleIDs", data=1 + np.arange(N_gas + N_warm_cyl) - ) + F["PartType0"].create_dataset("ParticleIDs", data=1 + np.arange(N_gas + N_warm_cyl)) F["PartType0"].create_dataset("InternalEnergy", data=u) if magnetic_field > 0.0: F["PartType0"].create_dataset("MagneticField", data=B_cyl) diff --git a/MakeCloud_Bondi.py b/MakeCloud_Bondi.py deleted file mode 100644 index aaeb10f..0000000 --- a/MakeCloud_Bondi.py +++ /dev/null @@ -1,310 +0,0 @@ -#!/usr/bin/env python -""" -MakeCloud: "Believe me, we've got some very turbulent clouds, the best clouds. You're gonna love it." -This is an alternate version that creates a homogeneous spherical medium around a sink particle, meant to do a Bondi-accretion test. By default it uses default GIZMO units - -Usage: MakeCloud_Bondi.py [options] - -Options: - -h --help Show this screen. - --R_over_Rsonic= Outer radius of the cloud in Rsonic [default: 16.0] - --Rs_over_Rsonic= Ratio of sink radius to sonic radius [default: 1.0] - --rho_gas= Density of the gas around the sink at inifinity, in msolar/pc^3 [default: 0.01] - --filename= Name of the IC file to be generated - --N= Number of gas particles [default: 583200] - --MBH= Mass of the sink at center of the sphere, should be much larger than the gas mass, in msolar [default: 1.0] - --boxsize= Simulation box size in pc - --length_unit= Unit of length in pc [default: 1000] - --mass_unit= Unit of mass in M_sun [default: 1e10] - --v_unit= Unit of velocity in m/s [default: 1000] - --GMC_units Sets units appropriate for GMCs, so pc, m/s, m_sun, tesla - --localdir Changes directory defaults assuming all files are used from local directory. - --R_cut_out= Distance around the sink to be left empty in pc [default: 0.0001] -""" - -from __future__ import print_function -import numpy as np -from scipy import fftpack, interpolate, ndimage -from scipy.integrate import quad, odeint, solve_bvp -from scipy.spatial import cKDTree -from scipy.spatial.distance import cdist -from scipy.optimize import minimize -import sys -import h5py -import os -from docopt import docopt - -# from pykdgrav import Potential -arguments = docopt(__doc__) -GMC_units = arguments["--GMC_units"] -length_unit = float(arguments["--length_unit"]) -mass_unit = float(arguments["--mass_unit"]) -v_unit = float(arguments["--v_unit"]) -if GMC_units: - length_unit = 1.0 - mass_unit = 1.0 - v_unit = 1.0 - - -R_over_Rsonic = float(arguments["--R_over_Rsonic"]) -Rs_over_Rsonic = float(arguments["--Rs_over_Rsonic"]) -rho_gas = float(arguments["--rho_gas"]) / ( - mass_unit / (length_unit**3.0) -) # msolar/pc^3 to code units -N_gas = int(float(arguments["--N"]) + 0.5) -M_BH = float(arguments["--MBH"]) / mass_unit # mstar to code units -R_cut_out = float(arguments["--R_cut_out"]) / length_unit -filename = arguments["--filename"] -localdir = arguments["--localdir"] - - -# set gravitational constant -G = 4325.69 / length_unit / (v_unit**2) * (mass_unit) - -if localdir: - turb_path = "turb" - glass_path = "glass_256.npy" - -# get rsink radius from sonic radius -csound = 200 / v_unit # 200 m/s -rsonic_Bondi = G * M_BH / 2.0 / (csound**2) / length_unit -msonic = 4.0 * np.pi / 3.0 * rho_gas * (rsonic_Bondi**3.0) -Rsink = Rs_over_Rsonic * rsonic_Bondi -R = R_over_Rsonic * rsonic_Bondi -print("Radius set as %g pc" % (R * length_unit)) -if arguments["--boxsize"] is not None: - # print((arguments["--boxsize"])) - boxsize = float(arguments["--boxsize"]) * length_unit -else: - boxsize = 10 * R -center = np.ones(3) * boxsize / 2.0 -res_effective = int(N_gas ** (1.0 / 3.0) + 0.5) - -# Load Bondi Spherical solution -IC_data = np.loadtxt("BONDI_SOLUTION.dat") # load IC data from Hubber -data_r = IC_data[:, 0] * rsonic_Bondi # coordinate originally in Rsonic -data_within_R_index = data_r <= R -data_r = data_r[data_within_R_index] # restrict to within R radius -data_mass_in_r = ( - IC_data[data_within_R_index, 1] * msonic -) # gas mass in units of 4*pi*rho0*Rsonic^3/3 -data_density = ( - IC_data[data_within_R_index, 2] * rho_gas -) # density originally in rho0 -data_vr = ( - IC_data[data_within_R_index, 3] * csound -) # velocity originally in csound -M_gas = np.max(data_mass_in_r) # normalize to total mass -# Get glass -mgas = np.repeat(M_gas / N_gas, N_gas) # gas mass init -x = 2 * (np.load(glass_path) - 0.5) # load glass to have basic structure -Nx = len(x[:, 0]) -# original glass size -r = np.sum(x**2, axis=1) ** 0.5 # calculate radius -x = x[r.argsort()][ - :N_gas -] # sort the particles by radius, it should be spherically symmetric, now take the right number of particles -x *= (float(Nx) / N_gas * 4 * np.pi / 3 / 8) ** ( - 1.0 / 3 -) # scale up the sphere to have unit radius -r = np.sum(x**2, axis=1) ** 0.5 # recalculate radius -# Strech the particles based on how much mass is in a given radius -int_mass = np.cumsum(mgas) # integrated mass -newr = np.interp( - int_mass, data_mass_in_r, data_r -) # calculate the radius the particles should be at to have the sane mass within the radius a sthe solution -stretch = newr / r # stretch factor -x[:, 0] *= stretch -x[:, 1] *= stretch -x[:, 2] *= stretch -u = ( - np.ones_like(mgas) * 0.101 * ((1000 / v_unit) ** 2) / 2.0 -) # /2 needed because it is molecular -rho = np.interp( - int_mass, data_mass_in_r, data_density -) # interpolate the density -h = (32 * mgas / rho) ** (1.0 / 3) -vr = np.interp(int_mass, data_mass_in_r, data_vr) # interpolate the velocity -v = -x -v[:, 0] *= vr / newr -v[:, 1] *= vr / newr -v[:, 2] *= vr / newr -# print chek data -print( - "mdot avg: %g std %g" - % ( - np.mean(data_density * data_r * data_r * data_vr * np.pi * 4.0), - np.std(data_density * data_r * data_r * data_vr * np.pi * 4.0), - ) -) -print("density") -print(data_density) -print("density_alt") -altdens = ( - np.diff(data_mass_in_r) - / np.diff(data_r) - / (4.0 * np.pi * data_r[1:] * data_r[1:]) -) -print(altdens) -print("relative") -print(data_density[1:] / altdens) -# Calculate NEWSNK parameters -print("Calculating t_radial...") -for z in [0.125, 0.5, 1.0, 2.0, 4.0]: - print("\t At %g Rsonic = %g pc" % (z, (rsonic_Bondi * z * length_unit))) - ind = newr <= (z * rsonic_Bondi) - print("\t mass enclosed %g in msun" % (np.sum(mgas[ind]) * mass_unit)) - # tot_weight=np.sum(mgas[ind]/rho[ind]) - tot_weight = 4.0 / 3.0 * np.pi * np.max(newr[ind]) ** 3 - x_dot_v = np.sum(x * v, axis=1) - t_rad = ( - -np.sum(mgas[ind]) - * tot_weight - / np.sum(4.0 * np.pi * (x_dot_v * newr * mgas)[ind]) - ) - print("\t t_rad=%g in code units" % (t_rad)) - print("\t m/t_rad=%g in code units" % (np.sum(mgas[ind]) / t_rad)) - # from IC - ind = np.arange(len(data_r))[(data_r <= (z * rsonic_Bondi))] - ind2 = np.argmax(data_r[ind]) - dm_ic = np.diff(data_mass_in_r) - t_rad_IC = ( - data_mass_in_r[ind2] * (4.0 / 3.0 * np.pi) * data_r[ind2] ** 3 - ) / ( - 4.0 - * np.pi - * np.sum(data_r[ind] * data_r[ind] * data_vr[ind] * dm_ic[ind]) - ) - print("\t t_rad_IC=%g in code units" % (t_rad_IC)) - print( - "\t m_IC/t_rad_IC=%g in code units" % (data_mass_in_r[ind2] / t_rad_IC) - ) -# Calculate flux -print("Calculating density ...") -dr = np.diff(newr) -rho_alt = (mgas[1:] / dr) / (4.0 * np.pi * (newr[1:] ** 2)) -print(rho) -print(rho_alt) -print(rho[1:] / rho_alt) -# keep the one sthat are not too close -ind = newr > R_cut_out -print( - "Removing %d particles that are inside R_cut_out of %g" - % ((N_gas - np.sum(ind)), R_cut_out) -) -N_gas = np.sum(ind) -M_gas = np.sum(mgas[ind]) - -NGBvals = [1, 2, 5, 10, 20, 50, 100, 200, 400, 800, 1600] -for ngb in NGBvals: - print("Radius at Ngbfactor of %g is %g pc" % (ngb, newr[32 * ngb])) - -# center coordinates -x = x + boxsize / 2.0 -print("Writing snapshot...") - -if filename is None: - filename = ( - "rho%3.2g_" % (rho_gas * (mass_unit / (length_unit**3))) - + ("MBH%g_" % (mass_unit * M_BH) if M_BH > 0 else "") - + "R_over_Rsonic%g_Res%d_Rs_over_Rsonic%g" - % (R_over_Rsonic, res_effective, Rs_over_Rsonic) - + ".hdf5" - ) - filename = filename.replace("+", "").replace("e0", "e") - filename = "".join(filename.split()) - -F = h5py.File(filename, "w") -F.create_group("PartType0") -F.create_group("Header") -F["Header"].attrs["NumPart_ThisFile"] = [ - N_gas, - 0, - 0, - 0, - 0, - (1 if M_BH > 0 else 0), -] -F["Header"].attrs["NumPart_Total"] = [ - N_gas, - 0, - 0, - 0, - 0, - (1 if M_BH > 0 else 0), -] -F["Header"].attrs["MassTable"] = [M_gas / N_gas, 0, 0, 0, 0, M_BH] -F["Header"].attrs["BoxSize"] = boxsize -F["Header"].attrs["Time"] = 0.0 -F["PartType0"].create_dataset("Masses", data=mgas[ind]) -F["PartType0"].create_dataset("Coordinates", data=x[ind, :]) -F["PartType0"].create_dataset("Velocities", data=v[ind, :]) -F["PartType0"].create_dataset("ParticleIDs", data=np.arange(N_gas) + 1) -F["PartType0"].create_dataset("InternalEnergy", data=u[ind]) -F["PartType0"].create_dataset("Density", data=rho[ind]) -F["PartType0"].create_dataset("SmoothingLength", data=h[ind]) -if M_BH > 0: - F.create_group("PartType5") - F["PartType5"].create_dataset("Masses", data=np.array([M_BH])) - F["PartType5"].create_dataset("BH_Mass", data=np.array([M_BH])) - F["PartType5"].create_dataset("Coordinates", data=center) - F["PartType5"].create_dataset("Velocities", data=v[0, :] * 0.0) - F["PartType5"].create_dataset("ParticleIDs", data=np.array([N_gas + 1])) - F["PartType5"].create_dataset("SinkRadius", data=np.array([Rsink])) - -F.close() - -if GMC_units: - delta_m = M_gas / N_gas - rhocrit = 421 / delta_m**2 - rho_avg = 3 * M_gas / (R**3) / (4 * np.pi) - softening = 0.000173148 # 100AU/2.8 #(delta_m/rhocrit)**(1./3) - ncrit = 1.0e11 # 8920 / delta_m**2 - tff = 8.275e-3 * rho_avg**-0.5 - mdot_Bondi = ( - np.exp(1.5) * np.pi * (G**2) * (M_BH**2) * rho_gas / (csound**3) - ) - tend_Bondi = 2.0 * (2.0 * G * M_BH / (csound**3)) - print( - "Bondi accretion parameters: \n \t Mgas:\t\t", - M_gas, - "\n \t Mdot:\t\t", - mdot_Bondi, - "\n \t Rsonic:\t", - rsonic_Bondi, - "\n \t Rsink:\t\t", - Rsink, - "\n \t t_end:\t\t", - tend_Bondi, - "\n \t m_acc:\t\t", - tend_Bondi * mdot_Bondi, - "\n \t f_acc:\t\t", - tend_Bondi * mdot_Bondi, - ) - paramsfile = str( - open( - os.path.realpath(__file__).replace( - "MakeCloud_Bondi.py", "params.txt" - ), - "r", - ).read() - ) - - replacements = { - "NAME": "../ICs/" + filename.replace(".hdf5", ""), - "DTSNAP": tend_Bondi / 200, - "SOFTENING": softening, - "GASSOFT": 2.0e-8, - "TMAX": tend_Bondi, - "RHOMAX": ncrit, - "BOXSIZE": boxsize, - "OUTFOLDER": "output_%g" % (Rs_over_Rsonic), - } - - print(replacements["NAME"]) - # print(paramsfile) - for k in replacements.keys(): - paramsfile = paramsfile.replace(k, str(replacements[k])) - open("params_" + filename.replace(".hdf5", "") + ".txt", "w").write( - paramsfile - )