diff --git a/examples/rmvs_swifter_comparison/8pl_16tp_encounters/init_cond.py b/examples/rmvs_swifter_comparison/8pl_16tp_encounters/init_cond.py index 49d017b86..5ce3cec71 100755 --- a/examples/rmvs_swifter_comparison/8pl_16tp_encounters/init_cond.py +++ b/examples/rmvs_swifter_comparison/8pl_16tp_encounters/init_cond.py @@ -19,12 +19,13 @@ swiftest_cb = "cb.swiftest.in" swiftest_bin = "bin.swiftest.dat" swiftest_enc = "enc.swiftest.dat" +swiftest_disc = "discard.swiftest.out" sim = swiftest.Simulation() sim.param['T0'] = 0.0 sim.param['DT'] = 1.0 -sim.param['TSTOP'] = 365.25e1 +sim.param['TSTOP'] = 365.25e2 sim.param['ISTEP_OUT'] = 10 sim.param['ISTEP_DUMP'] = 10 sim.param['CHK_QMIN_COORD'] = "HELIO" @@ -84,7 +85,7 @@ for i in pl.id: pli = pl.sel(id=i) rstart = 2 * np.double(pli['Radius']) # Start the test particles at a multiple of the planet radius away - vstart = 1.5 * np.sqrt(2 * np.double(pli['Mass']) / rstart) # Start the test particle velocities at a multiple of the escape speed + vstart = 1.5 * np.sqrt(2 * np.double(pli['GMass']) / rstart) # Start the test particle velocities at a multiple of the escape speed xvstart = np.array([rstart / np.sqrt(2.0), rstart / np.sqrt(2.0), 0.0, vstart, 0.0, 0.0]) # The positions and velocities of each pair of test particles will be in reference to a planet plvec = np.array([np.double(pli['px']), @@ -124,6 +125,7 @@ sim.param['CB_IN'] = swiftest_cb sim.param['BIN_OUT'] = swiftest_bin sim.param['ENC_OUT'] = swiftest_enc +sim.param['DISCARD_OUT'] = swiftest_disc sim.save(swiftest_input) sim.param['PL_IN'] = swifter_pl diff --git a/examples/rmvs_swifter_comparison/8pl_16tp_encounters/param.swifter.in b/examples/rmvs_swifter_comparison/8pl_16tp_encounters/param.swifter.in index 36dd2060f..741d5ac99 100644 --- a/examples/rmvs_swifter_comparison/8pl_16tp_encounters/param.swifter.in +++ b/examples/rmvs_swifter_comparison/8pl_16tp_encounters/param.swifter.in @@ -1,6 +1,6 @@ ! VERSION Swifter parameter file converted from Swiftest T0 0.0 -TSTOP 3652.5 +TSTOP 36525.0 DT 1.0 ISTEP_OUT 10 ISTEP_DUMP 10 @@ -11,13 +11,13 @@ IN_TYPE ASCII PL_IN pl.swifter.in TP_IN tp.in BIN_OUT bin.swifter.dat -ENC_OUT enc.swifter.dat CHK_QMIN 0.004650467260962157 CHK_RMIN 0.004650467260962157 CHK_RMAX 1000.0 CHK_EJECT 1000.0 CHK_QMIN_COORD HELIO CHK_QMIN_RANGE 0.004650467260962157 1000.0 +ENC_OUT enc.swifter.dat EXTRA_FORCE NO BIG_DISCARD NO CHK_CLOSE YES diff --git a/examples/rmvs_swifter_comparison/8pl_16tp_encounters/param.swiftest.in b/examples/rmvs_swifter_comparison/8pl_16tp_encounters/param.swiftest.in index b08b66850..633613b13 100644 --- a/examples/rmvs_swifter_comparison/8pl_16tp_encounters/param.swiftest.in +++ b/examples/rmvs_swifter_comparison/8pl_16tp_encounters/param.swiftest.in @@ -1,6 +1,6 @@ ! VERSION Swiftest parameter input T0 0.0 -TSTOP 3652.5 +TSTOP 36525.0 DT 1.0 ISTEP_OUT 10 ISTEP_DUMP 10 @@ -12,7 +12,6 @@ PL_IN pl.swiftest.in TP_IN tp.in CB_IN cb.swiftest.in BIN_OUT bin.swiftest.dat -ENC_OUT enc.swiftest.dat CHK_QMIN 0.004650467260962157 CHK_RMIN 0.004650467260962157 CHK_RMAX 1000.0 @@ -22,7 +21,9 @@ CHK_QMIN_RANGE 0.004650467260962157 1000.0 MU2KG 1.988409870698051e+30 TU2S 86400 DU2M 149597870700.0 +ENC_OUT enc.swiftest.dat EXTRA_FORCE NO +DISCARD_OUT discard.swiftest.out BIG_DISCARD NO CHK_CLOSE YES RHILL_PRESENT YES @@ -31,6 +32,3 @@ ROTATION NO TIDES NO ENERGY NO GR NO -YARKOVSKY NO -YORP NO -MTINY 0.0 diff --git a/examples/rmvs_swifter_comparison/8pl_16tp_encounters/pl.in b/examples/rmvs_swifter_comparison/8pl_16tp_encounters/pl.in index 207dd84f6..e506d4743 100644 --- a/examples/rmvs_swifter_comparison/8pl_16tp_encounters/pl.in +++ b/examples/rmvs_swifter_comparison/8pl_16tp_encounters/pl.in @@ -1,33 +1,33 @@ 8 -1 4.9125474498983623693e-11 0.001475124456355905224 +1 4.9125474498983623693e-11 0.0014751244935834097723 1.6306381826061645943e-05 --0.30949970210807342674 0.1619004125820537876 0.041620272188990829754 --0.018820805516945871005 -0.023743802865467341506 -0.00021385162925667799668 -2 7.243452483873646905e-10 0.006759069616556246028 +-0.37723213134178928962 0.033583006382262102496 0.037347605390151236704 +-0.008342765357522348782 -0.026815373891535949219 -0.0014259978126226569104 +2 7.243452483873646905e-10 0.006759064884572926305 4.0453784346544178454e-05 --0.5567137338251560985 -0.46074173273652380134 0.02580196630219121906 -0.012753121506668980284 -0.015678149412530151263 -0.0009510907726656827677 -3 8.9970113821660187435e-10 0.010044908171483009529 +-0.48772324783973242113 -0.53438159383885197595 0.020810210727290629623 +0.014797158412775200664 -0.013730742274988789325 -0.0010423161661001279551 +3 8.9970113821660187435e-10 0.01004492566381686821 4.25875607065040958e-05 -0.6978790186886838498 -0.73607603319120218366 3.261671020506711323e-05 -0.012205130808798069983 0.0117727888369263504476 -6.0385404652521189453e-07 -4 9.549535102761465607e-11 0.0072466797341124641736 +0.75635861236241797023 -0.674672320414933413 2.9482164814591560215e-05 +0.0111724732704936106226 0.0127745538308708106445 -6.4922715710692383536e-07 +4 9.549535102761465607e-11 0.007246661019460395855 2.265740805092889601e-05 --1.617661473167097963 0.38314370807747849534 0.04771055403546069218 --0.0027036789764029569086 -0.012421968497550240837 -0.00019400613558421780209 -5 2.825345908631354893e-07 0.35527079166215922855 +-1.6298746849997449715 0.3207423054955682029 0.04670239471945906301 +-0.0021807443711779391643 -0.012535287802970449672 -0.00020920945169970240392 +5 2.825345908631354893e-07 0.3552709189640263194 0.00046732617030490929307 -4.1527454588897487753 -2.8347492039446908763 -0.081136554176388195336 -0.0041683967523185880624 0.0065946899141205552256 -0.00012065009272080269359 -6 8.459715183006415395e-08 0.43765832419088212185 +4.173466317007961557 -2.8016935252076420326 -0.081737437087388009616 +0.0041198995131368296857 0.0066275220486694394126 -0.000119701756809504006665 +6 8.459715183006415395e-08 0.43765971211410006393 0.00038925687730393611812 -6.39471595410062843 -7.621162747287802297 -0.121992225877669294154 -0.0039680130835247464163 0.0035798698934692090544 -0.00022010758050265331019 -7 1.2920249163736673626e-08 0.46960112247450473807 +6.414531809550852337 -7.603234575912382276 -0.12309230194363039723 +0.0039583234912346160553 0.003591394820962775982 -0.00021992254044702681153 +7 1.2920249163736673626e-08 0.4696151691794380732 0.00016953449859497231466 -14.793135356927480828 13.074218343364380601 -0.14311846037737518955 --0.0026297294662822792016 0.0027702756265410048361 4.4212949669357180555e-05 -8 1.5243589003230834323e-08 0.78136567314580814177 +14.779979482167510341 13.088063528768900667 -0.14289732557134240953 +-0.0026326202792006140295 0.002767798188238951903 4.424098699891899271e-05 +8 1.5243589003230834323e-08 0.7813830605782720197 0.000164587904124493665 -29.568629894896030663 -4.5543028991960081697 -0.58771107137394917874 -0.00046181040300440859715 0.0031288137434451902125 -7.498349850432879627e-05 +29.57093474400743105 -4.53865809805519671 -0.5880859062837571205 +0.00046012900387533010914 0.0031291067120909890273 -7.4950441013985698475e-05 diff --git a/examples/rmvs_swifter_comparison/8pl_16tp_encounters/pl.swifter.in b/examples/rmvs_swifter_comparison/8pl_16tp_encounters/pl.swifter.in index 3179473c0..1dbc75e21 100644 --- a/examples/rmvs_swifter_comparison/8pl_16tp_encounters/pl.swifter.in +++ b/examples/rmvs_swifter_comparison/8pl_16tp_encounters/pl.swifter.in @@ -2,35 +2,35 @@ 0 0.00029591220819207775568 0.0 0.0 0.0 0.0 0.0 0.0 -1 4.9125474498983623693e-11 0.001475124456355905224 +1 4.9125474498983623693e-11 0.0014751244935834097723 1.6306381826061645943e-05 --0.30949970210807342674 0.1619004125820537876 0.041620272188990829754 --0.018820805516945871005 -0.023743802865467341506 -0.00021385162925667799668 -2 7.243452483873646905e-10 0.006759069616556246028 +-0.37723213134178928962 0.033583006382262102496 0.037347605390151236704 +-0.008342765357522348782 -0.026815373891535949219 -0.0014259978126226569104 +2 7.243452483873646905e-10 0.006759064884572926305 4.0453784346544178454e-05 --0.5567137338251560985 -0.46074173273652380134 0.02580196630219121906 -0.012753121506668980284 -0.015678149412530151263 -0.0009510907726656827677 -3 8.9970113821660187435e-10 0.010044908171483009529 +-0.48772324783973242113 -0.53438159383885197595 0.020810210727290629623 +0.014797158412775200664 -0.013730742274988789325 -0.0010423161661001279551 +3 8.9970113821660187435e-10 0.01004492566381686821 4.25875607065040958e-05 -0.6978790186886838498 -0.73607603319120218366 3.261671020506711323e-05 -0.012205130808798069983 0.0117727888369263504476 -6.0385404652521189453e-07 -4 9.549535102761465607e-11 0.0072466797341124641736 +0.75635861236241797023 -0.674672320414933413 2.9482164814591560215e-05 +0.0111724732704936106226 0.0127745538308708106445 -6.4922715710692383536e-07 +4 9.549535102761465607e-11 0.007246661019460395855 2.265740805092889601e-05 --1.617661473167097963 0.38314370807747849534 0.04771055403546069218 --0.0027036789764029569086 -0.012421968497550240837 -0.00019400613558421780209 -5 2.825345908631354893e-07 0.35527079166215922855 +-1.6298746849997449715 0.3207423054955682029 0.04670239471945906301 +-0.0021807443711779391643 -0.012535287802970449672 -0.00020920945169970240392 +5 2.825345908631354893e-07 0.3552709189640263194 0.00046732617030490929307 -4.1527454588897487753 -2.8347492039446908763 -0.081136554176388195336 -0.0041683967523185880624 0.0065946899141205552256 -0.00012065009272080269359 -6 8.459715183006415395e-08 0.43765832419088212185 +4.173466317007961557 -2.8016935252076420326 -0.081737437087388009616 +0.0041198995131368296857 0.0066275220486694394126 -0.000119701756809504006665 +6 8.459715183006415395e-08 0.43765971211410006393 0.00038925687730393611812 -6.39471595410062843 -7.621162747287802297 -0.121992225877669294154 -0.0039680130835247464163 0.0035798698934692090544 -0.00022010758050265331019 -7 1.2920249163736673626e-08 0.46960112247450473807 +6.414531809550852337 -7.603234575912382276 -0.12309230194363039723 +0.0039583234912346160553 0.003591394820962775982 -0.00021992254044702681153 +7 1.2920249163736673626e-08 0.4696151691794380732 0.00016953449859497231466 -14.793135356927480828 13.074218343364380601 -0.14311846037737518955 --0.0026297294662822792016 0.0027702756265410048361 4.4212949669357180555e-05 -8 1.5243589003230834323e-08 0.78136567314580814177 +14.779979482167510341 13.088063528768900667 -0.14289732557134240953 +-0.0026326202792006140295 0.002767798188238951903 4.424098699891899271e-05 +8 1.5243589003230834323e-08 0.7813830605782720197 0.000164587904124493665 -29.568629894896030663 -4.5543028991960081697 -0.58771107137394917874 -0.00046181040300440859715 0.0031288137434451902125 -7.498349850432879627e-05 +29.57093474400743105 -4.53865809805519671 -0.5880859062837571205 +0.00046012900387533010914 0.0031291067120909890273 -7.4950441013985698475e-05 diff --git a/examples/rmvs_swifter_comparison/8pl_16tp_encounters/pl.swiftest.in b/examples/rmvs_swifter_comparison/8pl_16tp_encounters/pl.swiftest.in index 207dd84f6..e506d4743 100644 --- a/examples/rmvs_swifter_comparison/8pl_16tp_encounters/pl.swiftest.in +++ b/examples/rmvs_swifter_comparison/8pl_16tp_encounters/pl.swiftest.in @@ -1,33 +1,33 @@ 8 -1 4.9125474498983623693e-11 0.001475124456355905224 +1 4.9125474498983623693e-11 0.0014751244935834097723 1.6306381826061645943e-05 --0.30949970210807342674 0.1619004125820537876 0.041620272188990829754 --0.018820805516945871005 -0.023743802865467341506 -0.00021385162925667799668 -2 7.243452483873646905e-10 0.006759069616556246028 +-0.37723213134178928962 0.033583006382262102496 0.037347605390151236704 +-0.008342765357522348782 -0.026815373891535949219 -0.0014259978126226569104 +2 7.243452483873646905e-10 0.006759064884572926305 4.0453784346544178454e-05 --0.5567137338251560985 -0.46074173273652380134 0.02580196630219121906 -0.012753121506668980284 -0.015678149412530151263 -0.0009510907726656827677 -3 8.9970113821660187435e-10 0.010044908171483009529 +-0.48772324783973242113 -0.53438159383885197595 0.020810210727290629623 +0.014797158412775200664 -0.013730742274988789325 -0.0010423161661001279551 +3 8.9970113821660187435e-10 0.01004492566381686821 4.25875607065040958e-05 -0.6978790186886838498 -0.73607603319120218366 3.261671020506711323e-05 -0.012205130808798069983 0.0117727888369263504476 -6.0385404652521189453e-07 -4 9.549535102761465607e-11 0.0072466797341124641736 +0.75635861236241797023 -0.674672320414933413 2.9482164814591560215e-05 +0.0111724732704936106226 0.0127745538308708106445 -6.4922715710692383536e-07 +4 9.549535102761465607e-11 0.007246661019460395855 2.265740805092889601e-05 --1.617661473167097963 0.38314370807747849534 0.04771055403546069218 --0.0027036789764029569086 -0.012421968497550240837 -0.00019400613558421780209 -5 2.825345908631354893e-07 0.35527079166215922855 +-1.6298746849997449715 0.3207423054955682029 0.04670239471945906301 +-0.0021807443711779391643 -0.012535287802970449672 -0.00020920945169970240392 +5 2.825345908631354893e-07 0.3552709189640263194 0.00046732617030490929307 -4.1527454588897487753 -2.8347492039446908763 -0.081136554176388195336 -0.0041683967523185880624 0.0065946899141205552256 -0.00012065009272080269359 -6 8.459715183006415395e-08 0.43765832419088212185 +4.173466317007961557 -2.8016935252076420326 -0.081737437087388009616 +0.0041198995131368296857 0.0066275220486694394126 -0.000119701756809504006665 +6 8.459715183006415395e-08 0.43765971211410006393 0.00038925687730393611812 -6.39471595410062843 -7.621162747287802297 -0.121992225877669294154 -0.0039680130835247464163 0.0035798698934692090544 -0.00022010758050265331019 -7 1.2920249163736673626e-08 0.46960112247450473807 +6.414531809550852337 -7.603234575912382276 -0.12309230194363039723 +0.0039583234912346160553 0.003591394820962775982 -0.00021992254044702681153 +7 1.2920249163736673626e-08 0.4696151691794380732 0.00016953449859497231466 -14.793135356927480828 13.074218343364380601 -0.14311846037737518955 --0.0026297294662822792016 0.0027702756265410048361 4.4212949669357180555e-05 -8 1.5243589003230834323e-08 0.78136567314580814177 +14.779979482167510341 13.088063528768900667 -0.14289732557134240953 +-0.0026326202792006140295 0.002767798188238951903 4.424098699891899271e-05 +8 1.5243589003230834323e-08 0.7813830605782720197 0.000164587904124493665 -29.568629894896030663 -4.5543028991960081697 -0.58771107137394917874 -0.00046181040300440859715 0.0031288137434451902125 -7.498349850432879627e-05 +29.57093474400743105 -4.53865809805519671 -0.5880859062837571205 +0.00046012900387533010914 0.0031291067120909890273 -7.4950441013985698475e-05 diff --git a/examples/rmvs_swifter_comparison/8pl_16tp_encounters/tp.in b/examples/rmvs_swifter_comparison/8pl_16tp_encounters/tp.in index c1e239467..e8eeff31e 100644 --- a/examples/rmvs_swifter_comparison/8pl_16tp_encounters/tp.in +++ b/examples/rmvs_swifter_comparison/8pl_16tp_encounters/tp.in @@ -1,49 +1,49 @@ 16 101 --0.30947664140174180325 0.16192347328838543885 0.041620272188990829754 --0.01621725604493672035 -0.023743802865467341506 -0.00021385162925667799668 +-0.37720907063545766613 0.033606067088593753744 0.037347605390151236704 +-0.0057392158855131963913 -0.026815373891535949219 -0.0014259978126226569104 102 --0.30952276281440505024 0.16187735187572213635 0.041620272188990829754 --0.021424354988955021661 -0.023743802865467341506 -0.00021385162925667799668 +-0.37725519204812091312 0.033559945675930451248 0.037347605390151236704 +-0.010946314829531501173 -0.026815373891535949219 -0.0014259978126226569104 103 --0.55665652353468386693 -0.46068452244605162527 0.02580196630219121906 -0.019100355212014374223 -0.015678149412530151263 -0.0009510907726656827677 +-0.48766603754926024505 -0.53432438354837974437 0.020810210727290629623 +0.021144392118120594604 -0.013730742274988789325 -0.0010423161661001279551 104 --0.5567709441156283301 -0.4607989430269959774 0.02580196630219121906 -0.0064058878013235863447 -0.015678149412530151263 -0.0009510907726656827677 +-0.4877804581302045972 -0.53443880412932420754 0.020810210727290629623 +0.008449924707429806725 -0.013730742274988789325 -0.0010423161661001279551 105 -0.6979392465946233637 -0.7360158052852626698 3.261671020506711323e-05 -0.019099571043071944532 0.0117727888369263504476 -6.0385404652521189453e-07 +0.7564188402683574841 -0.6746120925089938991 2.9482164814591560215e-05 +0.018066913504767485171 0.0127745538308708106445 -6.4922715710692383536e-07 106 -0.6978187907827443359 -0.73613626109714169754 3.261671020506711323e-05 -0.005310690574524194567 0.0117727888369263504476 -6.0385404652521189453e-07 +0.75629838445647845635 -0.67473254832087292687 2.9482164814591560215e-05 +0.0042780330362197352065 0.0127745538308708106445 -6.4922715710692383536e-07 107 --1.6176294307533440886 0.38317575049123231423 0.04771055403546069218 -0.00037580012182093606998 -0.012421968497550240837 -0.00019400613558421780209 +-1.6298426425859910971 0.3207743479093220218 0.04670239471945906301 +0.00089873472704595381427 -0.012535287802970449672 -0.00020920945169970240392 108 --1.6176935155808518374 0.38311166566372467646 0.04771055403546069218 --0.005783158074626849887 -0.012421968497550240837 -0.00019400613558421780209 +-1.6299067274134988459 0.32071026308181438402 0.04670239471945906301 +-0.005260223469401832143 -0.012535287802970449672 -0.00020920945169970240392 109 -4.1534063578978459574 -2.834088304936593694 -0.081136554176388195336 -0.041050613953966016978 0.0065946899141205552256 -0.00012065009272080269359 +4.174127216016058739 -2.8010326261995448505 -0.081737437087388009616 +0.041002116714784257734 0.0066275220486694394126 -0.000119701756809504006665 110 -4.152084559881651593 -2.8354101029527880584 -0.081136554176388195336 --0.032713820449328842588 0.0065946899141205552256 -0.00012065009272080269359 +4.1728054179998643747 -2.8023544242157392148 -0.081737437087388009616 +-0.032762317688510601832 0.0066275220486694394126 -0.000119701756809504006665 111 -6.395266446455758924 -7.620612254932671803 -0.121992225877669294154 -0.026081181967058334609 0.0035798698934692090544 -0.00022010758050265331019 +6.415082301905982831 -7.602684083557251782 -0.12309230194363039723 +0.026071492374768204248 0.003591394820962775982 -0.00021992254044702681153 112 -6.394165461745497936 -7.621713239642932791 -0.121992225877669294154 --0.01814515580000884351 0.0035798698934692090544 -0.00022010758050265331019 +6.413981317195721843 -7.60378506826751277 -0.12309230194363039723 +-0.018154845392298973872 0.003591394820962775982 -0.00021992254044702681153 113 -14.793375114914683266 13.074458101351583039 -0.14311846037737518955 -0.0104650340723796142495 0.0027702756265410048361 4.4212949669357180555e-05 +14.780219240154712779 13.088303286756103105 -0.14289732557134240953 +0.010462143259461278988 0.002767798188238951903 4.424098699891899271e-05 114 -14.79289559894027839 13.073978585377178163 -0.14311846037737518955 --0.015724493004944172653 0.0027702756265410048361 4.4212949669357180555e-05 +14.779739724180307903 13.087823770781698229 -0.14289732557134240953 +-0.015727383817862507914 0.002767798188238951903 4.424098699891899271e-05 115 -29.568862657342247502 -4.5540701367497931074 -0.58771107137394917874 -0.0148974462162825404404 0.0031288137434451902125 -7.498349850432879627e-05 +29.571167506453647889 -4.538425335608981648 -0.5880859062837571205 +0.01489576481715346179 0.0031291067120909890273 -7.4950441013985698475e-05 116 -29.568397132449813824 -4.554535661642223232 -0.58771107137394917874 --0.013973825410273721728 0.0031288137434451902125 -7.498349850432879627e-05 +29.570701981561214211 -4.5388908605014117725 -0.5880859062837571205 +-0.013975506809402800379 0.0031291067120909890273 -7.4950441013985698475e-05 diff --git a/examples/symba_energy_momentum/param.supercatastrophic_headon.in b/examples/symba_energy_momentum/param.supercatastrophic_headon.in index e9b60e7da..a8ff7cbde 100644 --- a/examples/symba_energy_momentum/param.supercatastrophic_headon.in +++ b/examples/symba_energy_momentum/param.supercatastrophic_headon.in @@ -17,6 +17,7 @@ CHK_RMIN 0.005 CHK_RMAX 1e6 CHK_EJECT -1.0 ! ignore this check CHK_QMIN -1.0 ! ignore this check +DISCARD_OUT discard.supercatastrophic_headon.out ENC_OUT enc.supercatastrophic_headon.dat EXTRA_FORCE no ! no extra user-defined forces BIG_DISCARD no ! output all planets if anything discarded diff --git a/examples/symba_mars_disk/param.in b/examples/symba_mars_disk/param.in index 08a0979e6..e2be7eece 100644 --- a/examples/symba_mars_disk/param.in +++ b/examples/symba_mars_disk/param.in @@ -19,7 +19,6 @@ CHK_EJECT 3389500000.0 CHK_QMIN 3389500.0 CHK_QMIN_COORD HELIO CHK_QMIN_RANGE 3389500.0 338950000000.0 -ENC_OUT /dev/null EXTRA_FORCE no BIG_DISCARD no RHILL_PRESENT yes diff --git a/examples/symba_swifter_comparison/1pl_1pl_encounter/swiftest_vs_swifter.ipynb b/examples/symba_swifter_comparison/1pl_1pl_encounter/swiftest_vs_swifter.ipynb index 3a80eebd1..dee66b3b6 100644 --- a/examples/symba_swifter_comparison/1pl_1pl_encounter/swiftest_vs_swifter.ipynb +++ b/examples/symba_swifter_comparison/1pl_1pl_encounter/swiftest_vs_swifter.ipynb @@ -2,7 +2,7 @@ "cells": [ { "cell_type": "code", - "execution_count": 8, + "execution_count": 1, "metadata": {}, "outputs": [], "source": [ @@ -13,7 +13,7 @@ }, { "cell_type": "code", - "execution_count": 9, + "execution_count": 2, "metadata": {}, "outputs": [ { @@ -35,7 +35,7 @@ }, { "cell_type": "code", - "execution_count": 10, + "execution_count": 3, "metadata": {}, "outputs": [ { @@ -57,7 +57,7 @@ }, { "cell_type": "code", - "execution_count": 11, + "execution_count": 4, "metadata": {}, "outputs": [], "source": [ @@ -66,7 +66,7 @@ }, { "cell_type": "code", - "execution_count": 12, + "execution_count": 5, "metadata": {}, "outputs": [], "source": [ @@ -75,23 +75,23 @@ }, { "cell_type": "code", - "execution_count": 13, + "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "text/plain": [ - "[,\n", - " ]" + "[,\n", + " ]" ] }, - "execution_count": 13, + "execution_count": 6, "metadata": {}, "output_type": "execute_result" }, { "data": { - "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAEGCAYAAABy53LJAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAaRklEQVR4nO3df5BV5Z3n8fenG5AYUVHBNDYRJvQ4os4S7EX8USalwxYwWYkxSclmAxozDP4qs052h5qtnUlqaxJmMu64blgtjD/AzYbK5ocSC3UYY9asCUZMHAUZBkaNNLZCMBqJo/y43/3jnG6ulwvcc/uevpdzP6+qW33POc9z7vdcOP3t5zznPI8iAjMzs1p1NDsAMzM7ujhxmJlZJk4cZmaWiROHmZll4sRhZmaZjGh2AMPhlFNOiUmTJjU7DDOzo8rTTz/9q4gYV7m+LRLHpEmTWL9+fbPDMDM7qkj6ZbX1vlRlZmaZOHGYmVkmThxmZpZJW/RxVLN37176+vp45513mh1KU4wePZru7m5GjhzZ7FDM7CjTtomjr6+PMWPGMGnSJCQ1O5xhFRHs2rWLvr4+Jk+e3OxwzOwo07aXqt555x1OPvnktksaAJI4+eST27a1ZWZD07aJA2jLpDGgnY/dzIambS9VmZnl7V/27Oetd/ZSCtgfQakUREApIn2VvS8l72OgbFp+7/5gfynYWyqxf3+wrxTsK5WSdfuD/aXSgTL7k/X7SsG+dNvl07uZfMr7G3pcThzD6IILLuAnP/nJQeuvuuoqPvaxj/HJT36yCVGZWR727i9x4V/9kNd/u6epcUw/fawTx9GsWtIws2Las6/E67/dw+yzPsBHzhhHh5JLxB0SnR3QIaXLyfuO8vcdSdkRHaKzQ4zs7Eh+diQ/R3Qm20Z0dBx4n5YZMbg9Wc6DE8cwOu6449i9ezcRwY033sgPf/hDJk+ejGdhNCueUnpen3v6WObP+GCTo2mstu4cb5bvf//7bN68meeee44777zTLRGzAiqlfw8W8T4UJ44mePzxx5k/fz6dnZ1MmDCBSy65pNkhmVmDDVxJ6Chg5nDiaBLfDmtWbANXoHPqZmgqJ44muPjii1m1ahX79++nv7+fxx57rNkhmVmDDfRxdBQwc+SaOCTNlrRZ0lZJS6psl6Tb0u3PSpqerp8o6TFJmyRtlHRTWZ2TJK2VtCX9OTbPY8jD5ZdfTk9PD+eccw7XXnstH/nIR5odkpk12GAfR3PDyEVud1VJ6gSWAbOAPuApSasj4vmyYnOAnvR1HnB7+nMf8CcR8XNJY4CnJa1N6y4BHo2IpWkyWgL8aV7H0Ui7d+8GkstUX//615scjZnlaaCPo4iXpfNsccwAtkbECxGxB1gFzKsoMw9YGYl1wImSuiKiPyJ+DhARbwGbgNPK6qxI368APp7jMZiZ1WXgJnt3jmdzGrCtbLmPA7/8ay4jaRLwYeDJdNWpEdEPkP4cX+3DJS2StF7S+p07d9Z7DGZmdRns4yhe3sg1cVT7uiqfdDtsGUnHAd8FvhARv8ny4RGxPCJ6I6J33LiD5lo3M8tVafCuquJljjwTRx8wsWy5G3il1jKSRpIkjW9GxPfKyrwmqSst0wXsaHDcZmZDVipw73ieieMpoEfSZEmjgCuB1RVlVgML0rurZgJvRkS/kt6ku4BNEfHfqtRZmL5fCDyQ3yGYmdUnCtziyO2uqojYJ+kG4BGgE7g7IjZKWpxuvwNYA8wFtgJvA1en1S8EPgs8J+mZdN2fRcQaYCnwbUnXAC8Dn8rrGMzM6hUUt48j10EO01/0ayrW3VH2PoDrq9T7fxyigRcRu4BLGxtpc2zbto0FCxbw6quv0tHRwaJFi7jpppveUyYiuOmmm1izZg3HHnss9957L9OnT29SxGZWqyL3cXh03CYaMWIEt9xyC9OnT+ett97i3HPPZdasWUydOnWwzEMPPcSWLVvYsmULTz75JNdeey1PPvnkYfZqZq2gNPgcR5MDyYGHHGmirq6uwdbDmDFjOPPMM9m+fft7yjzwwAMsWLAAScycOZM33niD/v7+ZoRrZhkU+QFAtziAL/9gI8+/kulu3yOaOuF4/uLfnlVz+Zdeeolf/OIXnHfeee9Zv337diZOPHDjWXd3N9u3b6erq6thsZpZ45U8yKHlaffu3VxxxRXceuutHH/88e/ZVm2SpyL+BWNWNL6rquCytAwabe/evVxxxRV85jOf4ROf+MRB27u7u9m27cDD9X19fUyYMGE4QzSzOvjJcctFRHDNNddw5plncvPNN1ctc9lll7Fy5UoignXr1nHCCSf4MpXZUaA0eLWgeJnDLY4meuKJJ7jvvvs455xzmDZtGgBf+cpXePnllwFYvHgxc+fOZc2aNUyZMoVjjz2We+65p4kRm1mtijyRkxNHE1100UVV+zDKSWLZsmXDFJGZNUrJU8eamVkWgy2OAv6WLeAhmZk1X6nAz3E4cZiZ5aDAg+M6cZiZ5SHcx2FmZlkUeZBDJw4zsxyEHwC0PHzuc59j/PjxnH322YPrXn/9dWbNmkVPTw+zZs3i17/+9eC2r371q0yZMoUzzjiDRx55pOo+D1ffzIbPYB+HWxzWSFdddRUPP/zwe9YtXbqUSy+9lC1btnDppZeydOlSAJ5//nlWrVrFxo0befjhh7nuuuvYv3//Qfs8VH0zG17hYdUtDxdffDEnnXTSe9Y98MADLFyYzIy7cOFC7r///sH1V155JccccwyTJ09mypQp/OxnPzton4eqb2bDq8h9HH5yHOChJfDqc43d5wfOgTnZ/9p/7bXXBsei6urqYseOHUAyvPrMmTMHyw0Mr15rfTMbXh7k0JrOw6ubHV0Ghzgs4HnqFgfU1TLIy6mnnkp/fz9dXV309/czfvx4oPbh1Q9V38yGl1scNmwuu+wyVqxYAcCKFSuYN2/e4PpVq1bx7rvv8uKLL7JlyxZmzJhRc30zG15FnjrWiaOJ5s+fz/nnn8/mzZvp7u7mrrvuYsmSJaxdu5aenh7Wrl3LkiVLADjrrLP49Kc/zdSpU5k9ezbLli2js7MTgM9//vOsX78e4JD1zWx4lUrJzyK2OHSkYb2LoLe3NwZ+sQ7YtGkTZ555ZpMiag3+Dszy83cbX2XRfU/z4I0XcfZpJzQ7nLpIejoieivXu8VhZpaDA53jTQ0jF04cZmY58CCHBdUOl+kOpZ2P3Ww4HBhypLlx5KFtE8fo0aPZtWtXW/4CjQh27drF6NGjmx2KWWEVeerYtn2Oo7u7m76+Pnbu3NnsUJpi9OjRdHd3NzsMs8I6MORIc+PIQ9smjpEjRzJ58uRmh2FmBeXnOMzMLJMo8CCHThxmZjkY6OMoXtpw4jAzy0WRh1V34jAzy0HJEzmZmVkmAy2OAt5W5cRhZpYDD6tuZmaZDD45XsDucScOM7McuMVhZmaZ+AHAOkmaLWmzpK2SDppRSInb0u3PSppetu1uSTskbaio8yVJ2yU9k77m5nkMZmb1GBgFzy2ODCR1AsuAOcBUYL6kqRXF5gA96WsRcHvZtnuB2YfY/d9GxLT0taahgZuZNUCpVNxBDvNsccwAtkbECxGxB1gFVE6APQ9YGYl1wImSugAi4nHg9RzjMzPLjYdVr89pwLay5b50XdYy1dyQXtq6W9LYagUkLZK0XtL6dh0B18yap+Q+jrpU+7YqJ7+opUyl24EPAdOAfuCWaoUiYnlE9EZE77hx446wSzOzxooCD6ueZ+LoAyaWLXcDr9RR5j0i4rWI2B8RJeBOkktiZmYtJXAfRz2eAnokTZY0CrgSWF1RZjWwIL27aibwZkT0H26nA30gqcuBDYcqa2bWLEUe5DC3iZwiYp+kG4BHgE7g7ojYKGlxuv0OYA0wF9gKvA1cPVBf0reAjwKnSOoD/iIi7gL+WtI0kktaLwF/nNcxmJnVq8iDHOY6A2B6q+yainV3lL0P4PpD1J1/iPWfbWSMZmZ5CN9VZWZmWUS4j8PMzDIoch+HE4eZWQ48yKGZmWVy4Mnx4mUOJw4zsxxERCE7xsGJw8wsFxHF7N8AJw4zs1yUIgrZvwFOHGZmuShFMfs3wInDzCwXEVHA2cYTThxmZjlILlUVM3U4cZiZ5SDpHG92FPlw4jAzy0HJd1WZmVkWpYjqU9UVgBOHmVkOwn0cZmaWRcl9HGZmlkXgFoeZmWXgBwDNzCwTD3JoZmaZlEru4zAzswz85LiZmWUS+AFAMzPLoOQ+DjMzyyICJw4zM6ud+zjMzCwTD3JoZmaZ+DkOMzPLJNziMDOzLEqeOtbMzLJw57iZmWVS8u24ZmaWhfs4zMwsk4igo6C/YQt6WGZmzZV0jrdpi0PS+CrrzsgnHDOzYmj3qWN/LOnTAwuS/gT4fn4hmZkd/ZJBDouZOUbUUOajwHJJnwJOBTYBM/IMysysCNq2xRER/cDDwPnAJGBlROzOOS4zs6NakZ/jOGKLQ9JaoB84G+gG7pb0eER8Me/gzMyOVqVSez/H8RDwZxHxRkRsAC4A3qxl55JmS9osaaukJVW2S9Jt6fZnJU0v23a3pB2SNlTUOUnSWklb0p9ja4nFzGw4FbmPo5bEMQZ4RNKPJV0PnBwR//VIlSR1AsuAOcBUYL6kqRXF5gA96WsRcHvZtnuB2VV2vQR4NCJ6gEfTZTOzlhLtfFdVRHw5Is4CrgcmAP9X0t/XsO8ZwNaIeCEi9gCrgHkVZeaR9JlERKwDTpTUlX7u48DrVfY7D1iRvl8BfLyGWMzMhlVQ3D6OLA8A7gBeBXYBBz3bUcVpwLay5b50XdYylU5NO+wHOu6rxiJpkaT1ktbv3LmzhnDNzBqnrSdyknStpB+RXBY6BfijiPj9GvZd7RuLOsrUJSKWR0RvRPSOGzeuEbs0M6tZqcATOdXyHMfpwBci4pmM++4DJpYtdwOv1FGm0muSuiKiP72stSNjXGZmuUtGxy1m5qilj2NJHUkD4CmgR9JkSaOAK4HVFWVWAwvSu6tmAm8OXIY6jNXAwvT9QuCBOmIzM8tVRLRv53i9ImIfcAPwCMnT5t+OiI2SFktanBZbA7wAbAXuBK4bqC/pW8BPgTMk9Um6Jt20FJglaQswK102M2spRR5WvZZLVXWLiDUkyaF83R1l74Pkbq1qdecfYv0u4NIGhmlm1nAltzjMzCyLUkD1+3+Ofk4cZmY5cB+HmZllUuRBDp04zMxyEIGnjjUzs9q1+yCHZmaWUURRu8adOMzMcuE+DjMzyyRo42HVzcwsO7c4zMwsk2TqWCcOMzOrURR4WHUnDjOzHJTaeepYMzPLzlPHmplZJm09kZOZmWXnQQ7NzCyTpMXR7Cjy4cRhZpYDP8dhZmaZFHnqWCcOM7MclPwch5mZZZGMjlvMzOHEYWaWg5LvqjIzsyxKEXQUNHM4cZiZ5SB8O66ZmWXhu6rMzCyTUkRBu8adOMzMcuEHAM3MLBMPq25mZjWLCMCj45qZWY3SvOFLVWZmVpvSYIujyYHkxInDzKzBSoMtjubGkRcnDjOzBiu5j8PMzOrhPg4zM6vJQIvDl6rMzKwmA30cBW1wOHGYmTXagRZHMTOHE4eZWYNFKfnpznEzM6tJ4D6OukmaLWmzpK2SllTZLkm3pduflTT9SHUlfUnSdknPpK+5eR6DmVlWJT85Xh9JncAyYA4wFZgvaWpFsTlAT/paBNxeY92/jYhp6WtNXsdgZlYPPzlevxnA1oh4ISL2AKuAeRVl5gErI7EOOFFSV411zcxakh8ArN9pwLay5b50XS1ljlT3hvTS1t2Sxlb7cEmLJK2XtH7nzp31HoOZWWbhIUfqVu0rixrLHK7u7cCHgGlAP3BLtQ+PiOUR0RsRvePGjaspYDOzRij66Lgjctx3HzCxbLkbeKXGMqMOVTciXhtYKelO4MHGhWxmNnR+crx+TwE9kiZLGgVcCayuKLMaWJDeXTUTeDMi+g9XN+0DGXA5sCHHYzAzy2ywj6Ogs47n1uKIiH2SbgAeATqBuyNio6TF6fY7gDXAXGAr8DZw9eHqprv+a0nTSC5dvQT8cV7HYGZWjyj4kCN5XqoivVV2TcW6O8reB3B9rXXT9Z9tcJhmZg3lIUfMzCyTwc7xgv6GLehhmZk1j1scZmaWSanywYOCceIwM2uwcIvDzMyy8CCHZmaWiYdVNzOzTEqeyMnMzLLwsOpmZpZJ0Qc5dOIwM2swD3JoZmaZDDzG4RaHmZnVxH0cZmaWSXjqWDMzy6LkqWPNzCwL31VlZmaZuI/DzMwy8bDqZmaWyeDUsc0NIzdOHGZmDTbY4iho77gTh5lZg4XvqjIzsyxKfo7DzMyy8O24ZmaWyWCLo8lx5MWJw8yswTx1rJmZZRJ+ANDMzLJwi8PMzDJxi8PMzDJxi8PMzDLx1LFmZpbJwNSxfgDQzMxqEm5xmJlZFh5yxMzMMimVkp9ucZiZWU08kZOZmWVyoHO8qWHkxonDzKzBwi0OMzPLYuABwILmDScOM7NGcx/HEEiaLWmzpK2SllTZLkm3pduflTT9SHUlnSRpraQt6c+xeR6DmVlWbnHUSVInsAyYA0wF5kuaWlFsDtCTvhYBt9dQdwnwaET0AI+my2ZmraPgLY4ROe57BrA1Il4AkLQKmAc8X1ZmHrAykp6kdZJOlNQFTDpM3XnAR9P6K4AfAX+axwE8s3wxeu25PHZtZgU2tRT8+YiJiD9odii5yDNxnAZsK1vuA86rocxpR6h7akT0A0REv6Tx1T5c0iKSVgwf/OAH6zqA943qRKM666prZu3trHHHc9L7RzU7jFzkmTiqtdGixjK11D2siFgOLAfo7e3NVHfAGVctq6eamVmh5dk53gdMLFvuBl6psczh6r6WXs4i/bmjgTGbmdkR5Jk4ngJ6JE2WNAq4ElhdUWY1sCC9u2om8GZ6GepwdVcDC9P3C4EHcjwGMzOrkNulqojYJ+kG4BGgE7g7IjZKWpxuvwNYA8wFtgJvA1cfrm6666XAtyVdA7wMfCqvYzAzs4Np4NH4Iuvt7Y3169c3Owwzs6OKpKcjordyvZ8cNzOzTJw4zMwsEycOMzPLxInDzMwyaYvOcUk7gV/WWf0U4FcNDCcPrR6j4xu6Vo+x1eOD1o+xFeM7PSLGVa5si8QxFJLWV7uroJW0eoyOb+haPcZWjw9aP8ZWj6+cL1WZmVkmThxmZpaJE8eRLW92ADVo9Rgd39C1eoytHh+0foytHt8g93GYmVkmbnGYmVkmThxmZpZJWycOSbMlbZa0VdJBc5enw73flm5/VtL0Wus2Mz5JEyU9JmmTpI2Sbmql+Mq2d0r6haQH84hvqDGmUxl/R9I/pt/l+S0W339I/303SPqWpNGNjq/GGH9P0k8lvSvpi1nqNjO+4TpPhhJj2fbcz5VMIqItXyTDtf8z8DvAKOAfgKkVZeYCD5HMSDgTeLLWuk2OrwuYnr4fA/xTK8VXtv1m4H8DD7bav3G6bQXw+fT9KODEVomPZHrlF4H3pcvfBq5q0nc4HvjXwF8CX8xSt8nx5X6eDDXG4TpXsr7aucUxA9gaES9ExB5gFTCvosw8YGUk1gEnKpl1sJa6TYsvIvoj4ucAEfEWsInkF01LxAcgqRv4Q+AbDY6rITFKOh64GLgLICL2RMQbrRJfum0E8D5JI4BjOXiGzWGJMSJ2RMRTwN6sdZsZ3zCdJ0OKEYbtXMmknRPHacC2suU+Dv5Pc6gytdRtZnyDJE0CPgw82WLx3Qr8J6DU4Lhq/fwjlfkdYCdwT3qJ4BuS3t8q8UXEduBvSCYz6yeZPfPvGhxfrTHmUbdWDfmMHM8TGHqMt5L/uZJJOycOVVlXeW/yocrUUneohhJfslE6Dvgu8IWI+E0DYzviZx+ujKSPATsi4ukGx1RpKN/hCGA6cHtEfBj4LdDoa/RD+Q7HkvzVOhmYALxf0r9vcHyH/PxhqFurIX9GzucJDCHGYTxXMmnnxNEHTCxb7ubgpv6hytRSt5nxIWkkycnwzYj4XoNjG2p8FwKXSXqJpNl+iaT/1WIx9gF9ETHwF+h3SBJJq8T3B8CLEbEzIvYC3wMuaHB8tcaYR91aDekzhuE8gaHFOFznSjbN7mRp1ovkL8oXSP5iG+iwOquizB/y3o7Jn9Vat8nxCVgJ3NqK319FmY+SX+f4kGIEfgyckb7/EvC1VokPOA/YSNK3IZKO/Bub8R2Wlf0S7+18bonz5DDx5X6eDDXGim25nSuZj6nZATT14JM7Vv6J5I6H/5yuWwwsTt8LWJZufw7oPVzdVokPuIikKfws8Ez6mtsq8VXsI9eTYYj/xtOA9en3eD8wtsXi+zLwj8AG4D7gmCZ9hx8g+av6N8Ab6fvjW+g8qRrfcJ0nQ/0Oh+tcyfLykCNmZpZJO/dxmJlZHZw4zMwsEycOMzPLxInDzMwyceIwM7NMnDjMMkpHzb2ubHmCpO/k9Fkfl/TnRyjzN5IuyePzzarx7bhmGaXjGj0YEWcPw2f9BLgsIn51mDKnA3dGxL/JOx4zcIvDrB5LgQ9JekbS1yRNkrQBQNJVku6X9ANJL0q6QdLN6UCJ6ySdlJb7kKSHJT0t6ceSfq/yQyT9LvBuRPxK0ph0fyPTbcdLeknSyIj4JXCypA8M43dgbcyJwyy7JcA/R8S0iPiPVbafDfw7kuG0/xJ4O5KBEn8KLEjLLCcZIuRc4IvA/6yynwuB8mG/f0QyBAnAlcB3IxmnirTchUM8LrOajGh2AGYF9Fj6i/4tSW8CP0jXPwf8fjoa6wXA/5EGB049psp+ukiGdh/wDZLhte8Hrgb+qGzbDpJRcs1y58Rh1njvlr0vlS2XSM65DuCNiJh2hP38C3DCwEJEPJFeFvsI0BkRG8rKjk7Lm+XOl6rMsnuLZKrRukQy58OLkj4Fg/OK/6sqRTcBUyrWrQS+BdxTsf53SQY7NMudE4dZRhGxC3hC0gZJX6tzN58BrpH0DyTDo1ebUvVx4MMqu54FfBMYS5I8gME5JaaQjORrljvfjmvWwiT9d+AHEfH36fIngXkR8dmyMpcD0yPivzQpTGsz7uMwa21fIZm0CUn/A5hDMrdDuRHALcMcl7UxtzjMzCwT93GYmVkmThxmZpaJE4eZmWXixGFmZpk4cZiZWSb/H/QYAJv507pwAAAAAElFTkSuQmCC\n", + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAERCAYAAABl3+CQAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAgUklEQVR4nO3de3Rd5X3m8e9jybZsSeZiXXyRjWVJNjaXOEQ1AbIIBMiA2+IGJhk8mQk0Sb3IrdPJpFN3WNO0zUpCFsmapI3bjpM0A2kSr4TWJQnGYEJS0qSUGMLFYIyNbbBsgWWDMTa+yfrNH+fICHPOsSydc/Y+0vNZS8tn7/3ufX4SHD3a+333uxURmJmZ5TMm6QLMzCzdHBRmZlaQg8LMzApyUJiZWUEOCjMzK8hBYWZmBY3YoJD095J2SVpfpON9SdL67Nd/KsYxzcwqwYgNCuD/AVcX40CSfhu4AFgAXAj8saRJxTi2mVnajdigiIgHgZcHrpPUJmmNpEck/ULS2YM83HzgXyKiNyIOAI9TpBAyM0u7ERsUeawAPhUR7wA+A/zNIPd7HLhG0kRJDcDlwIwS1WhmlirVSRdQLpLqgIuBH0rqXz0+u+064C9z7LYjIv5DRNwn6beAXwE9wL8BvaWv2swseRrJcz1JmgX8JCLOzfYpbIyIqUU47veAf4iI1cM9lplZ2o2aS08RsQ/YKun9AMp422D2lVQlaXL29fnA+cB9JSvWzCxFRuwZhaTvA5cBDcBLwGeBB4C/BaYCY4GVEZHrktOJx6oBHs0u7gNujojHil+1mVn6jNigMDOz4hg1l57MzGxoRuSop4aGhpg1a1bSZZiZVYxHHnlkd0Q05to2IoNi1qxZrFu3LukyzMwqhqTn823zpSczMyvIQWFmZgU5KMzMrCAHhZmZFeSgMDOzghwUZmZWkIPCzMwKSjQoJF0taaOkzZKW5dguSX+V3f6EpAuSqNPM7GR27z/MD9dtT7qMkkgsKCRVAcuBa8g8QW6JpPknNLsG6Mh+LSUzoZ+ZWer8yZ1P8Md3PsG23QeSLqXokrwzeyGwOSK2AEhaCSwGnh7QZjFwR2RmLnxI0umSpkZEdykKemz7Xvo8SaKZnaIN3fv46TO7AHjmxX3MaqhNuKLiSjIopgMDz9O6gAsH0WY6UJKgeOIbNzOHbaU4tJmNYG3AV+rb+cz+JTzz4mtcfe6wn4+WKkkGhXKsO/HP+cG0yTSUlpK5PMXMmTOHVNAV85qpfWXPkPY1s9Ht7TOn8dcbJrLxxdeSLqXokgyKLmDGgOUWYOcQ2gAQESuAFQCdnZ1Dun40fcnXhrKbmRkAZ7/yyIgMiiRHPf0a6JDUKmkccAPwoxPa/Aj4UHb00zuBV0vVP2FmNlxzp9Szbc8BDh09lnQpRZVYUEREL/BJ4F5gA/CDiHhK0s2Sbs42Ww1sATYD3wA+nkixZmaDcPaUevoCNr20P+lSiirR51FExGoyYTBw3d8NeB3AJ8pdl5nZUMydUg9kRj6d13JawtUUj+/MNjMrkrMm11IzdsyI66dwUJiZFUnVGNHRVM8zDgozM8tn7hQHhZmZFXD2lHp27z/Mnv2Hky6laBwUZmZFdPaUSQAjqp/CQWFmVkRvjHxyUJiZWQ6N9eOZXDvOZxRmZpbf3Cn1PPOSg8LMzPKYO6WeZ198jb6+kfHYAgeFmVmRnT2lnoNHj/HCy68nXUpROCjMzIpsbnbk00jp0HZQmJkV2ZzmOqSRM0TWQWFmVmQTx1Vz1pkT2fjSvqRLKQoHhZlZCYykqTwcFGZmJTB3yiS27R4ZDzFyUJiZlUD/Q4w276r8hxg5KMzMSmAkTeXhoDAzK4FZk2sZXz2GZ7orv0PbQWFmVgJVY0RHcx0bR8BUHg4KM7MSOXvKJF96GipJZ0paK2lT9t8z8rTbJulJSY9JWlfuOs3MhuPsKfX0vHaYlw8cSbqUYUnqjGIZ8NOI6AB+ml3O5/KIWBARneUpzcysON7o0K7sfoqkgmIxcHv29e3A7yVUh5lZyfQHRaVP5ZFUUDRHRDdA9t+mPO0CuE/SI5KWFjqgpKWS1kla19PTU+RyzcxOXWPdeM4cAQ8xqi7VgSXdD0zJsemWUzjMJRGxU1ITsFbSMxHxYK6GEbECWAHQ2dk5MiaBN7OKJok5I2DkU8mCIiKuzLdN0kuSpkZEt6SpwK48x9iZ/XeXpFXAQiBnUJiZpdGc5npWPbqDiEBS0uUMSVKXnn4E3Jh9fSNw14kNJNVKqu9/DbwXWF+2Cs3MiqCjuZ7XDvfS/eqhpEsZsqSC4lbgKkmbgKuyy0iaJml1tk0z8K+SHgceBu6OiDWJVGtmNkRzmzMd2s9W8OWnkl16KiQi9gBX5Fi/E1iUfb0FeFuZSzMzK6o5zXVAJigum5tv3E66+c5sM7MSOn3iOBrrx/PsS5U7i6yDwsysxOY017Gpgi89OSjMzEpsTnM9z760n76+yhy576AwMyuxOc31HDx6jB17DyZdypA4KMzMSmxgh3YlclCYmZVYx/EhspXZoe2gMDMrsUk1Y5l6Wo3PKMzMLL+O5noHhZmZ5TenqY7Nu/ZzrAJHPjkozMzKYM6Ueg739vHCy68nXcopc1CYmZXBnAqe88lBYWZWBh1NmSGylXiHtoPCzKwMasdXM/30CRU5RNZBYWZWJnOnVObIJweFmVmZdDTXsaXnAL3H+pIu5ZQ4KMzMymROUz1HjvWxbU9ljXxyUJiZlUmljnxyUJiZlUl7Ux2Sg8LMzPKYMK6KmWdOZFOFjXxyUJiZlVFHU+WNfEokKCS9X9JTkvokdRZod7WkjZI2S1pWzhrNzEphTnMdW3cf4Ehv5Yx8SuqMYj1wHfBgvgaSqoDlwDXAfGCJpPnlKc/MrDTmTqmnty/YuvtA0qUMWiJBEREbImLjSZotBDZHxJaIOAKsBBaXvjozs9LpaKq8kU9p7qOYDmwfsNyVXZeTpKWS1kla19PTU/LizMyGYnZjLWMqbORTdakOLOl+YEqOTbdExF2DOUSOdXknco+IFcAKgM7Ozsqb8N3MRoWasVXMaqh1UABExJXDPEQXMGPAcguwc5jHNDNL3JwKG/mU5ktPvwY6JLVKGgfcAPwo4ZrMzIZtTnMd2/Yc4NDRY0mXMihJDY99n6Qu4CLgbkn3ZtdPk7QaICJ6gU8C9wIbgB9ExFNJ1GtmVkwdzfX0BTzXUxk33pXs0lMhEbEKWJVj/U5g0YDl1cDqMpZmZlZyHc2Zhxht3rWfc6adlnA1J5fmS09mZiNSa0Nm5NNzuyrjjMJBYWZWZuOrqzhrci2bHBRmZpZPe1Odg8LMzPLraKpj2+4DHK2Ap905KMzMEtDeVEdvX/D8nvTP+eSgMDNLQP+cT5XwbAoHhZlZAtqaaoHMENm0c1CYmSVg4rhqpp8+oSI6tB0UZmYJ6WiujJFPDgozs4R0NNWxpWc/x/rSPeG1g8LMLCHtTXUc7u2j65XXky6lIAeFmVlC2itk5JODwswsIe1N2ckBUz6LrIPCzCwhp00YS1P9eJ9RmJlZfh3NdWzele6n3TkozMwS1NFUz+Zd+4lI78gnB4WZWYLam+o4cOQY3a8eSrqUvBwUZmYJ6u/QTvONdw4KM7MEdTS98VjUtHJQmJklaHLdeM6sHZfqDu1EgkLS+yU9JalPUmeBdtskPSnpMUnrylmjmVm5tDfWpXqIbFJnFOuB64AHB9H28ohYEBF5A8XMrJK1ZycHTOvIp0SCIiI2RMTGJN7bzCxtOprqePXgUXbvP5J0KTmlvY8igPskPSJpaaGGkpZKWidpXU9PT5nKMzMbvjdGPqWzn6JkQSHpfknrc3wtPoXDXBIRFwDXAJ+QdGm+hhGxIiI6I6KzsbFx2PWbmZVL/2NRn0vpyKfqUh04Iq4swjF2Zv/dJWkVsJDB9WuYmVWM5knjqR9fndp7KVJ76UlSraT6/tfAe8l0gpuZjSiSaGuqS+29FEkNj32fpC7gIuBuSfdm10+TtDrbrBn4V0mPAw8Dd0fEmiTqNTMrtY6m9D4WtWSXngqJiFXAqhzrdwKLsq+3AG8rc2lmZonoaK7jh4908errRzlt4tiky3mT1F56MjMbTd54iFH6Rj6dNCgkNeVYN7c05ZiZjU4dKX4s6mDOKH4h6QP9C5L+BzkuG5mZ2dBNP30CNWPHpLKfYjB9FJcBKyS9n0wH8wYyw1TNzKxIxowRbY3p7NA+6RlFRHQDa8iMUJoF3BER6ftOzMwqXFtjHVt60vfrdTB9FGuBC4FzyYxI+j+SvlzqwszMRpv2pjp27D3IwSPHki7lTQbTR3EP8L8iYm9ErAcuBl4tbVlmZqNPW2MdEbBld7rOKgYTFPXAvZJ+IekTwOSI+FyJ6zIzG3XammoBeK7nQMKVvNlg+ij+IiLOAT4BTAP+RdL9Ja/MzGyUmTW5ljFK3+SAp3LD3S7gRWAP8JZ7K8zMbHhqxlYx48yJPJeyDu3BdGZ/TNLPgZ8CDcAfRMT5pS7MzGw0amtM3+SAg7mP4izgjyLisRLXYmY26rU11vLLzbs51hdUjVHS5QCD66NY5pAwMyuPtsY6Dvf2sXPvwaRLOc6TApqZpUjb8ckB03P5yUFhZpYi7Y2ZoEjTyCcHhZlZipxRO44za8elauSTg8LMLGXaGmt5bld6brpzUJiZpUxbY53PKMzMLL+2xjr2HDjCKweOJF0KkFBQSLpN0jOSnpC0StLpedpdLWmjpM2SlpW5TDOzRPQ/FjUtZxVJnVGsBc7N3uH9LPCnJzaQVAUsB64B5gNLJM0va5VmZgloa3RQEBH3RURvdvEhoCVHs4XA5ojYEhFHgJXA4nLVaGaWlOlnTGBc9ZjUzCKbhj6KD5N55sWJpgPbByx3ZdflJGmppHWS1vX09BS5RDOz8qkaI2Y31KbmXoqSBYWk+yWtz/G1eECbW4Be4Lu5DpFjXeR7v4hYERGdEdHZ2Ng4/G/AzCxBbU11qbk7ezCTAg5JRFxZaLukG4HfAa6IiFwB0AXMGLDcAuwsXoVmZunV1ljHPU92c+joMWrGViVaS1Kjnq4G/gS4NiJez9Ps10CHpFZJ44AbgB+Vq0YzsyS1NdbSF/D8nny/IssnqT6Kr5N5xOpaSY9J+jsASdMkrQbIdnZ/ErgX2AD8ICKeSqheM7OyStPIp5JdeiokItrzrN8JLBqwvBpYXa66zMzSYnZj5vnZaXiIURpGPZmZ2Qkmjqtm+ukTUnFG4aAwM0uptqZ0zPnkoDAzS6n+WWT7+vLeGVAWDgozs5Rqa6zj4NFjdO87lGgdDgozs5RqS8nT7hwUZmYplZZZZB0UZmYp1VA3jkk11Q4KMzPLTVJmzidfejIzs3wyj0VNdrpxB4WZWYq1N9XR89phXj14NLEaHBRmZimWhjmfHBRmZinWP+fT1gQvPzkozMxSbMYZE6kaI7bs9hmFmZnlMK56DDPPnMgWn1GYmVk+sxtqHRRmZpbf7MZatu45wLGEJgd0UJiZpdzsxjqO9Paxc+/BRN7fQWFmlnKzGzIjn7bsTubyk4PCzCzlWrNDZLckdC+Fg8LMLOUa68ZTP746sQ7t6iTeVNJtwO8CR4DngN+PiL052m0DXgOOAb0R0VnGMs3MUkESsxtrE7uXIqkzirXAuRFxPvAs8KcF2l4eEQscEmY2ms1urEvsjCKRoIiI+yKiN7v4ENCSRB1mZpVidkMt3a8e4vUjvSdvXGRp6KP4MHBPnm0B3CfpEUlLCx1E0lJJ6ySt6+npKXqRZmZJmp2dHHBrAiOfShYUku6XtD7H1+IBbW4BeoHv5jnMJRFxAXAN8AlJl+Z7v4hYERGdEdHZ2NhY1O/FzCxprf1DZBO4/FSyzuyIuLLQdkk3Ar8DXBEROW83jIid2X93SVoFLAQeLHatZmZpl2RQJHLpSdLVwJ8A10bE63na1Eqq738NvBdYX74qzczSY8K4KqafPiGRkU9J9VF8HagH1kp6TNLfAUiaJml1tk0z8K+SHgceBu6OiDXJlGtmlrzZjclMDpjIfRQR0Z5n/U5gUfb1FuBt5azLzCzNZjfU8o+P7iAikFS2903DqCczMxuE2Y117D/cS89rh8v6vg4KM7MK0d+h/VyZLz85KMzMKkT/87PL3aGdSB9FEo4ePUpXVxeHDh1KupSyq6mpoaWlhbFjxyZdipkNw7TTJlAzdkzZO7RHTVB0dXVRX1/PrFmzytoJlLSIYM+ePXR1ddHa2pp0OWY2DGPGiFmTa8s+3fioufR06NAhJk+ePKpCAjKzTk6ePHlUnkmZjURtjXVln8Zj1AQFMOpCot9o/b7NRqLZjbVsf+UgR3r7yvaeoyoozMwqXWtDLcf6ghdeLt9ZhYOiRC6++OKc62+66SbuvPPOMldjZiNF/yyy5Rwi66AokV/96ldJl2BmI9DxIbJlDIpRM+qp3Orq6ti/fz8Rwac+9SkeeOABWltbyTNRrpnZoEyqGUtD3Xi2lvFeCp9RlNiqVavYuHEjTz75JN/4xjd8pmFmw1buyQEdFCX24IMPsmTJEqqqqpg2bRrvec97ki7JzCpcW2MtW8o4RNZBUQYenmpmxdTaUMvLB46w9/UjZXk/B0WJXXrppaxcuZJjx47R3d3Nz372s6RLMrMKN7uhvCOf3JldYu973/t44IEHOO+885gzZw7vfve7ky7JzCrcGyOf9vOOs84o+fs5KEpk//7MiARJfP3rX0+4GjMbSWacOZHqMSrbVB6+9GRmVmHGVo1h5uSJZRv55KAwM6tAsxvqyvZcCgeFmVkFmt1Yy7Y9r3Osr/Q38SYSFJI+J+kJSY9Juk/StDztrpa0UdJmScvKXaeZWVrNbqjlSG8fO145WPL3SuqM4raIOD8iFgA/Af7sxAaSqoDlwDXAfGCJpPllrdLMLKWOTw5YhstPiQRFROwbsFgL5Dp3WghsjogtEXEEWAksLkd9ZmZp1z9EdmsZOrQT66OQ9HlJ24EPkuOMApgObB+w3JVdl+94SyWtk7Sup6enuMUWwfbt27n88suZN28e55xzDl/72tfe0iYi+MM//EPa29s5//zzefTRRxOo1MwqweTacUyqqS5Lh3bJgkLS/ZLW5/haDBARt0TEDOC7wCdzHSLHury9NhGxIiI6I6KzsbGxON9EEVVXV/OVr3yFDRs28NBDD7F8+XKefvrpN7W555572LRpE5s2bWLFihV87GMfS6haM0s7ScxurCvLENmS3XAXEVcOsun3gLuBz56wvguYMWC5BdhZhNL4ix8/xdM795284SmYP20Sn/3dc/Junzp1KlOnTgWgvr6eefPmsWPHDubPf6Pb5a677uJDH/oQknjnO9/J3r176e7uPr6fmdlAsxtq+dVze0r+PkmNeuoYsHgt8EyOZr8GOiS1ShoH3AD8qBz1ldq2bdv4zW9+w4UXXvim9Tt27GDGjDeysaWlhR07dpS7PDOrELMba3lx3yEOHO4t6fskNYXHrZLmAn3A88DNANlhst+MiEUR0Svpk8C9QBXw9xHxVDHevNBf/qW2f/9+rr/+er761a8yadKkN23L9VAjzzxrZvn0j3zauvsA504/rWTvk0hQRMT1edbvBBYNWF4NrC5XXaV29OhRrr/+ej74wQ9y3XXXvWV7S0sL27e/0X/f1dXFtGk5bzExM6O1ITvyqcRB4TuzyyQi+MhHPsK8efP49Kc/nbPNtddeyx133EFE8NBDD3Haaae5f8LM8po1ORMU20o8OaBnjy2TX/7yl3znO9/hvPPOY8GCBQB84Qtf4IUXXgDg5ptvZtGiRaxevZr29nYmTpzIt7/97QQrNrO0mzCuimmn1ZR8FlkHRZm8613vytkHMZAkli9fXqaKzGwkaC3DY1F96cnMrILNmlzLlp79J/1DdDgcFGZmFay1oZZ9h3p55fWjJXsPB4WZWQU7PudTCafycFCYmVWw1obMvRSlnMrDQWFmVsFazphA9RixbY+DwszMchhbNYaZZ04s6RBZB0UZffjDH6apqYlzzz33+LqXX36Zq666io6ODq666ipeeeWV49u++MUv0t7ezty5c7n33ntzHrPQ/mY2OsxqqPWlp5HipptuYs2aNW9ad+utt3LFFVewadMmrrjiCm699VYAnn76aVauXMlTTz3FmjVr+PjHP86xY8fecsx8+5vZ6NHaUMu2PQfoK9Hzs0fnDXf3LIMXnyzuMaecB9cU/iV96aWXsm3btjetu+uuu/j5z38OwI033shll13Gl770Je666y5uuOEGxo8fT2trK+3t7Tz88MNcdNFFg9rfzEaP1oZaDh3t48V9h5h2+oSiH99nFAl76aWXjs/nNHXqVHbt2gUMfsrxfPub2egxe8DkgKUwOs8oTvKXfxp4ynEzG6zWxjeC4pL2hqIf32cUCWtubqa7uxuA7u5umpqagMFPOZ5vfzMbPZrra5gwtqpkZxQOioRde+213H777QDcfvvtLF68+Pj6lStXcvjwYbZu3cqmTZtYuHDhoPc3s9FjzBgxq6HWQTESLFmyhIsuuoiNGzfS0tLCt771LZYtW8batWvp6Ohg7dq1LFu2DIBzzjmHD3zgA8yfP5+rr76a5cuXU1VVBcBHP/pR1q1bB5B3fzMbXVobSncvhUo542BSOjs7o/8Xab8NGzYwb968hCpK3mj//s1Guu8//AKPb9/LF687b0j9mZIeiYjOXNtGZ2e2mdkIs2ThTJYsnFmSYycSFJI+BywG+oBdwE3Z52Wf2G4b8BpwDOjNl3ZmZlY6SfVR3BYR50fEAuAnwJ8VaHt5RCwoRkiMxMtsgzFav28zK45EgiIi9g1YrAVK/puspqaGPXv2jLpfmhHBnj17qKmpSboUM6tQifVRSPo88CHgVeDyPM0CuE9SAP83IlYM9f1aWlro6uqip6dnqIeoWDU1NbS0tCRdhplVqJKNepJ0PzAlx6ZbIuKuAe3+FKiJiM/mOMa0iNgpqQlYC3wqIh7M835LgaUAM2fOfMfzzz9fjG/DzGxUKDTqKfHhsZLOAu6OiHNP0u7Pgf0R8eWTHTPX8FgzM8uvUFAk0kchqWPA4rXAMzna1Eqq738NvBdYX54KzcysX1J9FLdKmktmeOzzwM2QudQEfDMiFgHNwKrsjSPVwPciYk2e45mZWYkkfumpFCT1kAmgoWgAdhexnGJLe32Q/hpd3/Clvca01wfpq/GsiGjMtWFEBsVwSFqX5hv70l4fpL9G1zd8aa8x7fVBZdTYz5MCmplZQQ4KMzMryEHxVkO+qa9M0l4fpL9G1zd8aa8x7fVBZdQIuI/CzMxOwmcUZmZWkIPCzMwKGjVBIelqSRslbZb0lueFKuOvstufkHTBYPdNukZJMyT9TNIGSU9J+m9pqm/A9ipJv5H0k1LUN9waJZ0u6U5Jz2R/lhelrL7/nv3vu17S9yUVfUrgQdR3tqR/k3RY0mdOZd+ka0zR5yTvzzC7veSfk1MWESP+C6gCngNmA+OAx4H5J7RZBNwDCHgn8O+D3TcFNU4FLsi+rgeeLXaNw6lvwPZPA98DfpK2/87ZbbcDH82+Hgecnpb6gOnAVmBCdvkHZB74Ve76moDfAj4PfOZU9k1BjWn5nOSsr1yfk6F8jZYzioXA5ojYEhFHgJVknrA30GLgjsh4CDhd0tRB7ptojRHRHRGPAkTEa8AGMr9YUlEfgKQW4LeBbxa5rqLUKGkScCnwLYCIOBIRe9NSX3ZbNTBBUjUwEXjLUyFLXV9E7IqIXwNHT3XfpGtMy+ekwM+wXJ+TUzZagmI6sH3Achdv/R8kX5vB7Jt0jcdJmgW8Hfj3lNX3VeB/kpnfq1SGU+NsoAf4dva0/5vKTEaZivoiYgfwZeAFoBt4NSLuS6C+Uux7KoryPgl/Tgr5KqX/nJyy0RIUyrHuxHHB+doMZt9iGE6NmY1SHfCPwB/Fm58iWAxDrk/S7wC7IuKRItd0ouH8DKuBC4C/jYi3AweAYl9nH87P8Awyf5m2AtOAWkn/JYH6SrHvqRj2+6Tgc5J7x/J9Tk7ZaAmKLmDGgOUW3nranq/NYPZNukYkjSXzP/93I+KfUlbfJcC1kraRORV/j6R/SFmNXUBXRPT/hXknmeBIS31XAlsjoicijgL/BFycQH2l2PdUDOt9UvI5yadcn5NTl3QnSTm+yPy1uIXMX2P9HUznnNDmt3lzJ+LDg903BTUKuAP4ahp/hie0uYzSdWYPq0bgF8Dc7Os/B25LS33AhcBTZPomRKbj/VPlrm9A2z/nzR3FqfmcFKgxFZ+TfPWdsK1kn5MhfV9JF1C2bzQzmuRZMiMSbsmuuxm4ecD/RMuz258EOgvtm6YagXeROb19Angs+7UoLfWdcIySfgCG+d95AbAu+3P8Z+CMlNX3F2Qe8rUe+A4wPoH6ppD5q3kfsDf7elLKPic5a0zR5yTvz7Bcn5NT/fIUHmZmVtBo6aMwM7MhclCYmVlBDgozMyvIQWFmZgU5KMzMrCAHhVkB2RllPz5geZqkO0v0Xr8n6c9O0ubLkt5Tivc3y8fDY80KyM4J9JOIOLcM7/Ur4NqI2F2gzVnANyLivaWux6yfzyjMCrsVaJP0mKTbJM2StB5A0k2S/lnSjyVtlfRJSZ/OTir4kKQzs+3aJK2R9IikX0g6+8Q3kTQHOBwRuyXVZ483NrttkqRtksZGxPPAZElTyvgzsFHOQWFW2DLguYhYEBF/nGP7ucB/JjO99OeB1yMzqeC/AR/KtllBZrqNdwCfAf4mx3EuAQZOgf1zMtN5ANwA/GNk5ngi2+6SYX5fZoNWnXQBZhXuZ9lf7K9JehX4cXb9k8D52ZlKLwZ+KB2fWHR8juNMJTPNeb9vkplu+p+B3wf+YMC2XWRmkDUrCweF2fAcHvC6b8ByH5nP1xhgb0QsOMlxDgKn9S9ExC+zl7neDVRFxPoBbWuy7c3KwpeezAp7jcxjM4ckMs872Crp/XD8mdhvy9F0A9B+wro7gO8D3z5h/RwyEwOalYWDwqyAiNgD/FLSekm3DfEwHwQ+IulxMlOF53pE6IPA2zXg+hTwXeAMMmEBHH+eQjuZWW7NysLDY81SQtLXgB9HxP3Z5f8ILI6I/zqgzfuACyLifydUpo1C7qMwS48vkHlAEZL+GriGzLMNBqoGvlLmumyU8xmFmZkV5D4KMzMryEFhZmYFOSjMzKwgB4WZmRXkoDAzs4L+P9l9DczSV9flAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] @@ -108,7 +108,7 @@ }, { "cell_type": "code", - "execution_count": 19, + "execution_count": 14, "metadata": {}, "outputs": [ { @@ -466,140 +466,140 @@ " fill: currentColor;\n", "}\n", "
<xarray.DataArray 'vx' (time (y): 221)>\n",
-       "array([0.        , 0.        , 0.        , 0.        , 0.        ,\n",
-       "       0.        , 0.        , 0.        , 0.        , 0.        ,\n",
-       "       0.        , 0.        , 0.        , 0.        , 0.        ,\n",
-       "       0.        , 0.        , 0.        , 0.        , 0.        ,\n",
-       "       0.        , 0.        , 0.        , 0.        , 0.        ,\n",
-       "       0.        , 0.        , 0.        , 0.        , 0.        ,\n",
-       "       0.        , 0.        , 0.        , 0.        , 0.        ,\n",
-       "       0.        , 0.        , 0.        , 0.        , 0.        ,\n",
-       "       0.        , 0.        , 0.        , 0.        , 0.        ,\n",
-       "       0.        , 0.        , 0.        , 0.        , 0.        ,\n",
-       "       0.        , 0.        , 0.        , 0.        , 0.        ,\n",
-       "       0.        , 0.        , 0.        , 0.        , 0.        ,\n",
-       "       0.        , 0.        , 0.        , 0.        , 0.        ,\n",
-       "       0.        , 0.        , 0.        , 0.        , 0.        ,\n",
-       "       0.        , 0.        , 0.        , 0.        , 0.        ,\n",
-       "       0.        , 0.        , 0.        , 0.        , 0.        ,\n",
-       "       0.        , 0.        , 0.        , 0.        , 0.        ,\n",
-       "       0.        , 0.        , 0.        , 0.        , 0.        ,\n",
-       "       0.        , 0.        , 0.        , 0.        , 0.        ,\n",
-       "       0.        , 0.        , 0.        , 0.        , 0.        ,\n",
+       "array([ 0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,\n",
+       "        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,\n",
+       "        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,\n",
+       "        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,\n",
+       "        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,\n",
+       "        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,\n",
+       "        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,\n",
+       "        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,\n",
+       "        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,\n",
+       "        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,\n",
+       "        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,\n",
+       "        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,\n",
+       "        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,\n",
+       "        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,\n",
+       "        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,\n",
+       "        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,\n",
+       "        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,\n",
+       "        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,\n",
+       "        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,\n",
+       "        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,\n",
        "...\n",
-       "       0.        , 0.        , 0.        , 0.        , 0.        ,\n",
-       "       0.        , 0.        , 0.        , 0.        , 0.        ,\n",
-       "       0.        , 0.        , 0.        , 0.        , 0.        ,\n",
-       "       0.        , 0.        , 0.        , 0.        , 0.        ,\n",
-       "       0.        , 0.        , 0.        , 0.        , 0.        ,\n",
-       "       0.        , 0.        , 0.        , 0.        , 0.        ,\n",
-       "       0.        , 0.        , 0.        , 0.        , 0.        ,\n",
-       "       0.        , 0.        , 0.        , 0.        , 0.        ,\n",
-       "       0.        , 0.        , 0.        , 0.        , 0.        ,\n",
-       "       0.        , 0.        , 0.        , 0.        , 0.        ,\n",
-       "       0.        , 0.        , 0.        , 0.        , 0.        ,\n",
-       "       0.        , 0.        , 0.        , 0.        , 0.        ,\n",
-       "       0.        , 0.        , 0.        , 0.        , 0.        ,\n",
-       "       0.        , 0.        , 0.        , 0.        , 0.        ,\n",
-       "       0.        , 0.        , 0.02101973, 0.02102004, 0.02102057,\n",
-       "       0.0210213 , 0.02102224, 0.02102336, 0.02102465, 0.02102612,\n",
-       "       0.02102774, 0.02102951, 0.02103142, 0.02103346, 0.02103561,\n",
-       "       0.02103787, 0.02104022, 0.02104267, 0.02104519, 0.02104778,\n",
-       "       0.02105043, 0.02105312, 0.02105586, 0.02105862, 0.0210614 ,\n",
-       "       0.02106419])\n",
+       "        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,\n",
+       "        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,\n",
+       "        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,\n",
+       "        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,\n",
+       "        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,\n",
+       "        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,\n",
+       "        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,\n",
+       "        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,\n",
+       "        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,\n",
+       "        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,\n",
+       "        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,\n",
+       "        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,\n",
+       "        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,\n",
+       "        0.00000000e+00,  3.98117095e-11, -9.46069889e-11, -2.30539143e-10,\n",
+       "       -3.67965214e-10, -5.06869213e-10, -6.47232490e-10, -7.89038168e-10,\n",
+       "       -9.32266708e-10, -1.07690123e-09, -1.22292310e-09, -1.37031364e-09,\n",
+       "       -1.51905422e-09, -1.66912617e-09, -1.82051085e-09, -1.97318961e-09,\n",
+       "       -2.12714202e-09, -2.28234942e-09, -2.43879317e-09, -2.59645283e-09,\n",
+       "       -2.75530976e-09, -2.91534441e-09, -3.07653725e-09, -3.23886606e-09,\n",
+       "       -3.40231310e-09])\n",
        "Coordinates:\n",
        "    id        float64 2.0\n",
-       "  * time (y)  (time (y)) float64 0.0 0.0006845 0.001369 ... 0.1492 0.1499 0.1506
    • id
      ()
      float64
      2.0
      array(2.)
    • time (y)
      (time (y))
      float64
      0.0 0.0006845 ... 0.1499 0.1506
      array([0.      , 0.000684, 0.001369, ..., 0.149213, 0.149897, 0.150582])
    • " ], "text/plain": [ "\n", - "array([0. , 0. , 0. , 0. , 0. ,\n", - " 0. , 0. , 0. , 0. , 0. ,\n", - " 0. , 0. , 0. , 0. , 0. ,\n", - " 0. , 0. , 0. , 0. , 0. ,\n", - " 0. , 0. , 0. , 0. , 0. ,\n", - " 0. , 0. , 0. , 0. , 0. ,\n", - " 0. , 0. , 0. , 0. , 0. ,\n", - " 0. , 0. , 0. , 0. , 0. ,\n", - " 0. , 0. , 0. , 0. , 0. ,\n", - " 0. , 0. , 0. , 0. , 0. ,\n", - " 0. , 0. , 0. , 0. , 0. ,\n", - " 0. , 0. , 0. , 0. , 0. ,\n", - " 0. , 0. , 0. , 0. , 0. ,\n", - " 0. , 0. , 0. , 0. , 0. ,\n", - " 0. , 0. , 0. , 0. , 0. ,\n", - " 0. , 0. , 0. , 0. , 0. ,\n", - " 0. , 0. , 0. , 0. , 0. ,\n", - " 0. , 0. , 0. , 0. , 0. ,\n", - " 0. , 0. , 0. , 0. , 0. ,\n", - " 0. , 0. , 0. , 0. , 0. ,\n", + "array([ 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n", + " 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n", + " 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n", + " 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n", + " 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n", + " 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n", + " 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n", + " 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n", + " 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n", + " 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n", + " 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n", + " 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n", + " 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n", + " 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n", + " 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n", + " 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n", + " 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n", + " 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n", + " 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n", + " 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n", "...\n", - " 0. , 0. , 0. , 0. , 0. ,\n", - " 0. , 0. , 0. , 0. , 0. ,\n", - " 0. , 0. , 0. , 0. , 0. ,\n", - " 0. , 0. , 0. , 0. , 0. ,\n", - " 0. , 0. , 0. , 0. , 0. ,\n", - " 0. , 0. , 0. , 0. , 0. ,\n", - " 0. , 0. , 0. , 0. , 0. ,\n", - " 0. , 0. , 0. , 0. , 0. ,\n", - " 0. , 0. , 0. , 0. , 0. ,\n", - " 0. , 0. , 0. , 0. , 0. ,\n", - " 0. , 0. , 0. , 0. , 0. ,\n", - " 0. , 0. , 0. , 0. , 0. ,\n", - " 0. , 0. , 0. , 0. , 0. ,\n", - " 0. , 0. , 0. , 0. , 0. ,\n", - " 0. , 0. , 0.02101973, 0.02102004, 0.02102057,\n", - " 0.0210213 , 0.02102224, 0.02102336, 0.02102465, 0.02102612,\n", - " 0.02102774, 0.02102951, 0.02103142, 0.02103346, 0.02103561,\n", - " 0.02103787, 0.02104022, 0.02104267, 0.02104519, 0.02104778,\n", - " 0.02105043, 0.02105312, 0.02105586, 0.02105862, 0.0210614 ,\n", - " 0.02106419])\n", + " 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n", + " 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n", + " 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n", + " 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n", + " 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n", + " 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n", + " 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n", + " 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n", + " 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n", + " 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n", + " 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n", + " 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n", + " 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n", + " 0.00000000e+00, 3.98117095e-11, -9.46069889e-11, -2.30539143e-10,\n", + " -3.67965214e-10, -5.06869213e-10, -6.47232490e-10, -7.89038168e-10,\n", + " -9.32266708e-10, -1.07690123e-09, -1.22292310e-09, -1.37031364e-09,\n", + " -1.51905422e-09, -1.66912617e-09, -1.82051085e-09, -1.97318961e-09,\n", + " -2.12714202e-09, -2.28234942e-09, -2.43879317e-09, -2.59645283e-09,\n", + " -2.75530976e-09, -2.91534441e-09, -3.07653725e-09, -3.23886606e-09,\n", + " -3.40231310e-09])\n", "Coordinates:\n", " id float64 2.0\n", " * time (y) (time (y)) float64 0.0 0.0006845 0.001369 ... 0.1492 0.1499 0.1506" ] }, - "execution_count": 19, + "execution_count": 14, "metadata": {}, "output_type": "execute_result" } @@ -610,7 +610,7 @@ }, { "cell_type": "code", - "execution_count": 17, + "execution_count": 11, "metadata": {}, "outputs": [ { @@ -967,152 +967,80 @@ " stroke: currentColor;\n", " fill: currentColor;\n", "}\n", - "
    <xarray.DataArray 'vx' (time: 221)>\n",
-       "array([ 0.        , -0.02730963, -0.05461883, -0.08192718, -0.10923426,\n",
-       "       -0.13653965, -0.16384292, -0.19114364, -0.21844141, -0.24573578,\n",
-       "       -0.27302634, -0.30031266, -0.32759433, -0.35487091, -0.38214199,\n",
-       "       -0.40940715, -0.43666596, -0.463918  , -0.49116285, -0.51840009,\n",
-       "       -0.5456293 , -0.57285005, -0.60006193, -0.62726452, -0.6544574 ,\n",
-       "       -0.68164014, -0.70881234, -0.73597358, -0.76312342, -0.79026147,\n",
-       "       -0.8173873 , -0.8445005 , -0.87160064, -0.89868733, -0.92576014,\n",
-       "       -0.95281866, -0.97986247, -1.00689117, -1.03390434, -1.06090158,\n",
-       "       -1.08788246, -1.11484659, -1.14179356, -1.16872296, -1.19563437,\n",
-       "       -1.22252741, -1.24940165, -1.27625671, -1.30309216, -1.32990762,\n",
-       "       -1.35670269, -1.38347696, -1.41023003, -1.43696151, -1.463671  ,\n",
-       "       -1.4903581 , -1.51702243, -1.54366359, -1.57028119, -1.59687484,\n",
-       "       -1.62344416, -1.64998874, -1.67650822, -1.70300221, -1.72947032,\n",
-       "       -1.75591217, -1.78232739, -1.8087156 , -1.83507643, -1.8614095 ,\n",
-       "       -1.88771444, -1.91399088, -1.94023846, -1.96645681, -1.99264557,\n",
-       "       -2.01880437, -2.04493287, -2.0710307 , -2.09709752, -2.12313297,\n",
-       "       -2.1491367 , -2.17510838, -2.20104766, -2.2269542 , -2.25282767,\n",
-       "       -2.27866774, -2.30447407, -2.33024634, -2.35598424, -2.38168744,\n",
-       "       -2.40735563, -2.43298851, -2.45858576, -2.48414708, -2.50967219,\n",
-       "       -2.53516079, -2.56061259, -2.58602731, -2.61140468, -2.63674443,\n",
-       "...\n",
-       "       -3.28166908, -3.30592115, -3.33013189, -3.3543014 , -3.37842979,\n",
-       "       -3.40251722, -3.42656386, -3.45056991, -3.47453562, -3.49846126,\n",
-       "       -3.52234715, -3.54619364, -3.57000113, -3.59377005, -3.61750092,\n",
-       "       -3.64119426, -3.66485068, -3.68847086, -3.71205553, -3.73560548,\n",
-       "       -3.75912161, -3.78260489, -3.80605637, -3.82947723, -3.85286873,\n",
-       "       -3.87623226, -3.89956936, -3.92288168, -3.94617105, -3.96943947,\n",
-       "       -3.99268912, -4.01592242, -4.03914199, -4.06235073, -4.08555183,\n",
-       "       -4.10874879, -4.13194551, -4.15514625, -4.17835577, -4.20157933,\n",
-       "       -4.2248228 , -4.24809272, -4.2713964 , -4.29474206, -4.31813893,\n",
-       "       -4.34159746, -4.36512951, -4.38874856, -4.41247007, -4.43631182,\n",
-       "       -4.46029439, -4.48444172, -4.50878191, -4.53334814, -4.55817998,\n",
-       "       -4.58332498, -4.60884098, -4.63479915, -4.66128825, -4.68842081,\n",
-       "       -4.71634199, -4.7452432 , -4.77538326, -4.80712299, -4.84098468,\n",
-       "       -4.87776098, -4.91873117, -4.96614064, -5.02443531, -5.10428159,\n",
-       "       -5.24263186, -5.9750488 ,         nan,         nan,         nan,\n",
-       "               nan,         nan,         nan,         nan,         nan,\n",
-       "               nan,         nan,         nan,         nan,         nan,\n",
-       "               nan,         nan,         nan,         nan,         nan,\n",
-       "               nan,         nan,         nan,         nan,         nan,\n",
-       "               nan])\n",
+       "
    <xarray.DataArray 'vx' (time (y): 221)>\n",
+       "array([ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,\n",
+       "        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,\n",
+       "        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,\n",
+       "        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,\n",
+       "        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,\n",
+       "        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,\n",
+       "        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,\n",
+       "        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,\n",
+       "        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,\n",
+       "        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,\n",
+       "        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,\n",
+       "        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,\n",
+       "        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,\n",
+       "        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,\n",
+       "        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,\n",
+       "        0.,  0., nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan,\n",
+       "       nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan])\n",
        "Coordinates:\n",
-       "    id       float64 100.0\n",
-       "  * time     (time) float64 0.0 0.0006845 0.001369 ... 0.1492 0.1499 0.1506
    "
+       "    id        float64 100.0\n",
+       "  * time (y)  (time (y)) float64 0.0 0.0006845 0.001369 ... 0.1492 0.1499 0.1506
    " ], "text/plain": [ - "\n", - "array([ 0. , -0.02730963, -0.05461883, -0.08192718, -0.10923426,\n", - " -0.13653965, -0.16384292, -0.19114364, -0.21844141, -0.24573578,\n", - " -0.27302634, -0.30031266, -0.32759433, -0.35487091, -0.38214199,\n", - " -0.40940715, -0.43666596, -0.463918 , -0.49116285, -0.51840009,\n", - " -0.5456293 , -0.57285005, -0.60006193, -0.62726452, -0.6544574 ,\n", - " -0.68164014, -0.70881234, -0.73597358, -0.76312342, -0.79026147,\n", - " -0.8173873 , -0.8445005 , -0.87160064, -0.89868733, -0.92576014,\n", - " -0.95281866, -0.97986247, -1.00689117, -1.03390434, -1.06090158,\n", - " -1.08788246, -1.11484659, -1.14179356, -1.16872296, -1.19563437,\n", - " -1.22252741, -1.24940165, -1.27625671, -1.30309216, -1.32990762,\n", - " -1.35670269, -1.38347696, -1.41023003, -1.43696151, -1.463671 ,\n", - " -1.4903581 , -1.51702243, -1.54366359, -1.57028119, -1.59687484,\n", - " -1.62344416, -1.64998874, -1.67650822, -1.70300221, -1.72947032,\n", - " -1.75591217, -1.78232739, -1.8087156 , -1.83507643, -1.8614095 ,\n", - " -1.88771444, -1.91399088, -1.94023846, -1.96645681, -1.99264557,\n", - " -2.01880437, -2.04493287, -2.0710307 , -2.09709752, -2.12313297,\n", - " -2.1491367 , -2.17510838, -2.20104766, -2.2269542 , -2.25282767,\n", - " -2.27866774, -2.30447407, -2.33024634, -2.35598424, -2.38168744,\n", - " -2.40735563, -2.43298851, -2.45858576, -2.48414708, -2.50967219,\n", - " -2.53516079, -2.56061259, -2.58602731, -2.61140468, -2.63674443,\n", - "...\n", - " -3.28166908, -3.30592115, -3.33013189, -3.3543014 , -3.37842979,\n", - " -3.40251722, -3.42656386, -3.45056991, -3.47453562, -3.49846126,\n", - " -3.52234715, -3.54619364, -3.57000113, -3.59377005, -3.61750092,\n", - " -3.64119426, -3.66485068, -3.68847086, -3.71205553, -3.73560548,\n", - " -3.75912161, -3.78260489, -3.80605637, -3.82947723, -3.85286873,\n", - " -3.87623226, -3.89956936, -3.92288168, -3.94617105, -3.96943947,\n", - " -3.99268912, -4.01592242, -4.03914199, -4.06235073, -4.08555183,\n", - " -4.10874879, -4.13194551, -4.15514625, -4.17835577, -4.20157933,\n", - " -4.2248228 , -4.24809272, -4.2713964 , -4.29474206, -4.31813893,\n", - " -4.34159746, -4.36512951, -4.38874856, -4.41247007, -4.43631182,\n", - " -4.46029439, -4.48444172, -4.50878191, -4.53334814, -4.55817998,\n", - " -4.58332498, -4.60884098, -4.63479915, -4.66128825, -4.68842081,\n", - " -4.71634199, -4.7452432 , -4.77538326, -4.80712299, -4.84098468,\n", - " -4.87776098, -4.91873117, -4.96614064, -5.02443531, -5.10428159,\n", - " -5.24263186, -5.9750488 , nan, nan, nan,\n", - " nan, nan, nan, nan, nan,\n", - " nan, nan, nan, nan, nan,\n", - " nan, nan, nan, nan, nan,\n", - " nan, nan, nan, nan, nan,\n", - " nan])\n", + "\n", + "array([ 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,\n", + " 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,\n", + " 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,\n", + " 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,\n", + " 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,\n", + " 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,\n", + " 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,\n", + " 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,\n", + " 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,\n", + " 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,\n", + " 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,\n", + " 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,\n", + " 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,\n", + " 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,\n", + " 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,\n", + " 0., 0., nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan,\n", + " nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan])\n", "Coordinates:\n", - " id float64 100.0\n", - " * time (time) float64 0.0 0.0006845 0.001369 ... 0.1492 0.1499 0.1506" + " id float64 100.0\n", + " * time (y) (time (y)) float64 0.0 0.0006845 0.001369 ... 0.1492 0.1499 0.1506" ] }, - "execution_count": 17, + "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ - "swiftestsim.ds.sel(id=100)['vx']" + "swiftdiff.sel(id=100)['vx']" ] }, { "cell_type": "code", - "execution_count": 18, + "execution_count": 9, "metadata": {}, "outputs": [ { @@ -1513,7 +1441,7 @@ " nan])\n", "Coordinates:\n", " id int64 100\n", - " * time (time) float64 0.0 0.0006845 0.001369 ... 0.1492 0.1499 0.1506
    • id
      ()
      int64
      100
      array(100)
    • time
      (time)
      float64
      0.0 0.0006845 ... 0.1499 0.1506
      array([0.      , 0.000684, 0.001369, ..., 0.149213, 0.149897, 0.150582])
    • " ], "text/plain": [ "\n", @@ -1603,7 +1531,7 @@ " * time (time) float64 0.0 0.0006845 0.001369 ... 0.1492 0.1499 0.1506" ] }, - "execution_count": 18, + "execution_count": 9, "metadata": {}, "output_type": "execute_result" } diff --git a/python/swiftest/swiftest/io.py b/python/swiftest/swiftest/io.py index 0492a05f8..4e46ad808 100644 --- a/python/swiftest/swiftest/io.py +++ b/python/swiftest/swiftest/io.py @@ -285,8 +285,6 @@ def write_labeled_param(param, param_file_name): 'TP_IN', 'CB_IN', 'BIN_OUT', - 'ENC_OUT', - 'DISCARD_OUT', 'CHK_QMIN', 'CHK_RMIN', 'CHK_RMAX', @@ -303,7 +301,7 @@ def write_labeled_param(param, param_file_name): if val is not None: print(f"{key:<16} {val}", file=outfile) # Print the remaining key/value pairs in whatever order for key, val in ptmp.items(): - print(f"{key:<16} {val}", file=outfile) + if val != "": print(f"{key:<16} {val}", file=outfile) outfile.close() return diff --git a/src/io/io.f90 b/src/io/io.f90 index 2fa42dcc1..ea2b695b0 100644 --- a/src/io/io.f90 +++ b/src/io/io.f90 @@ -195,6 +195,7 @@ module subroutine io_dump_system(self, param) start = clock_count / (count_rate * 1.0_DP) finish = start lfirst = .false. + if (param%lenergy) call self%conservation_report(param, lterminal=.false.) else allocate(dump_param, source=param) param_file_name = trim(adjustl(DUMP_PARAM_FILE(idx))) @@ -224,7 +225,6 @@ module subroutine io_dump_system(self, param) if (param%lenergy) call self%conservation_report(param, lterminal=.true.) end if - return end subroutine io_dump_system @@ -574,7 +574,6 @@ module subroutine io_param_reader(self, unit, iotype, v_list, iostat, iomsg) iostat = -1 return end if - end if if (param%qmin > 0.0_DP) then if ((param%qmin_coord /= "HELIO") .and. (param%qmin_coord /= "BARY")) then @@ -622,6 +621,9 @@ module subroutine io_param_reader(self, unit, iotype, v_list, iostat, iomsg) write(*,*) "MU2KG = ",param%MU2KG write(*,*) "TU2S = ",param%TU2S write(*,*) "DU2M = ",param%DU2M + if (trim(adjustl(param%outfile)) == "") then + param%outfile = BIN_OUTFILE + end if if (trim(adjustl(param%enc_out)) /= "") then write(*,*) "ENC_OUT = ",trim(adjustl(param%enc_out)) else @@ -632,6 +634,7 @@ module subroutine io_param_reader(self, unit, iotype, v_list, iostat, iomsg) write(*,*) "BIG_DISCARD = ",param%lbig_discard else write(*,*) "! DISCARD_OUT not set: Discards will not be recorded to file" + param%lbig_discard = .false. write(*,*) "! BIG_DISCARD = ",param%lbig_discard end if @@ -643,10 +646,11 @@ module subroutine io_param_reader(self, unit, iotype, v_list, iostat, iomsg) ! Calculate the G for the system units param%GU = GC / (param%DU2M**3 / (param%MU2KG * param%TU2S**2)) - - ! Calculate the inverse speed of light in the system units - param%inv_c2 = einsteinC * param%TU2S / param%DU2M - param%inv_c2 = (param%inv_c2)**(-2) + if (param%lgr) then + ! Calculate the inverse speed of light in the system units + param%inv_c2 = einsteinC * param%TU2S / param%DU2M + param%inv_c2 = (param%inv_c2)**(-2) + end if associate(integrator => v_list(1)) if (integrator == RMVS) then @@ -911,7 +915,7 @@ module subroutine io_read_cb_in(self, param) end subroutine io_read_cb_in - function io_read_encounter(t, name1, name2, mass1, mass2, radius1, radius2, & + function io_read_encounter(t, id1, id2, Gmass1, Gmass2, radius1, radius2, & xh1, xh2, vh1, vh2, enc_out, out_type) result(ierr) !! author: David A. Minton !! @@ -921,8 +925,8 @@ function io_read_encounter(t, name1, name2, mass1, mass2, radius1, radius2, & !! Adapted from David E. Kaufmann's Swifter routine: io_read_encounter.f90 implicit none ! Arguments - integer(I4B), intent(out) :: name1, name2 - real(DP), intent(out) :: t, mass1, mass2, radius1, radius2 + integer(I4B), intent(out) :: id1, id2 + real(DP), intent(out) :: t, Gmass1, Gmass2, radius1, radius2 real(DP), dimension(:), intent(out) :: xh1, xh2, vh1, vh2 character(*), intent(in) :: enc_out, out_type ! Result @@ -947,12 +951,12 @@ function io_read_encounter(t, name1, name2, mass1, mass2, radius1, radius2, & return end if - read(iu, iostat = ierr) name1, xh1(1), xh1(2), xh1(3), vh1(1), vh1(2), vh1(3), mass1, radius1 + read(iu, iostat = ierr) id1, xh1(1), xh1(2), xh1(3), vh1(1), vh1(2), vh1(3), Gmass1, radius1 if (ierr /= 0) then close(unit = iu, iostat = ierr) return end if - read(iu, iostat = ierr) name2, xh2(2), xh2(2), xh2(3), vh2(2), vh2(2), vh2(3), mass2, radius2 + read(iu, iostat = ierr) id2, xh2(2), xh2(2), xh2(3), vh2(2), vh2(2), vh2(3), Gmass2, radius2 if (ierr /= 0) then close(unit = iu, iostat = ierr) return @@ -1238,6 +1242,8 @@ module subroutine io_write_discard(self, param) character(*), parameter :: PLNAMEFMT = '(I8, 2(1X, E23.16))' class(swiftest_body), allocatable :: pltemp + if (param%discard_out == "") return + associate(tp_discards => self%tp_discards, nsp => self%tp_discards%nbody, pl => self%pl, npl => self%pl%nbody) if (nsp == 0) return select case(param%out_stat) @@ -1283,30 +1289,23 @@ module subroutine io_write_discard(self, param) end subroutine io_write_discard - module subroutine io_write_encounter(t, name1, name2, mass1, mass2, radius1, radius2, & - xh1, xh2, vh1, vh2, enc_out, out_type) - !! author: David A. Minton - !! - !! Write close encounter data to output binary files - !! There is no direct file output from this subroutine - !! - !! Adapted from David E. Kaufmann's Swifter routine: io_write_encounter.f90 - !! Adapted from Hal Levison's Swift routine io_write_encounter.f + module subroutine io_write_encounter(self, pl, encbody, param) implicit none ! Arguments - integer(I4B), intent(in) :: name1, name2 - real(DP), intent(in) :: t, mass1, mass2, radius1, radius2 - real(DP), dimension(:), intent(in) :: xh1, xh2, vh1, vh2 - character(*), intent(in) :: enc_out, out_type + class(swiftest_encounter), intent(in) :: self !! Swiftest encounter list object + class(swiftest_pl), intent(in) :: pl !! Swiftest massive body object + class(swiftest_body), intent(in) :: encbody !! Encountering body - Swiftest generic body object (pl or tp) + class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters ! Internals - logical , save :: lfirst = .true. - integer(I4B), parameter :: lun = 30 - integer(I4B) :: ierr - integer(I4B), save :: iu = lun + logical , save :: lfirst = .true. + integer(I4B), parameter :: LUN = 30 + integer(I4B) :: k, ierr - open(unit = iu, file = enc_out, status = 'OLD', position = 'APPEND', form = 'UNFORMATTED', iostat = ierr) + if (param%enc_out == "" .or. self%nenc == 0) return + + open(unit = LUN, file = param%enc_out, status = 'OLD', position = 'APPEND', form = 'UNFORMATTED', iostat = ierr) if ((ierr /= 0) .and. lfirst) then - open(unit = iu, file = enc_out, status = 'NEW', form = 'UNFORMATTED', iostat = ierr) + open(unit = LUN, file = param%enc_out, status = 'NEW', form = 'UNFORMATTED', iostat = ierr) end if if (ierr /= 0) then write(*, *) "Swiftest Error:" @@ -1314,21 +1313,36 @@ module subroutine io_write_encounter(t, name1, name2, mass1, mass2, radius1, rad call util_exit(FAILURE) end if lfirst = .false. - write(iu, iostat = ierr) t - if (ierr < 0) then - write(*, *) "Swiftest Error:" - write(*, *) " Unable to write binary file record" - call util_exit(FAILURE) - end if - write(iu) name1, xh1(1), xh1(2), xh1(3), vh1(1), vh1(2), mass1, radius1 - write(iu) name2, xh2(1), xh2(2), xh2(3), vh2(1), vh2(2), mass2, radius2 - close(unit = iu, iostat = ierr) + + associate(ind1 => self%index1, ind2 => self%index2) + select type(encbody) + class is (swiftest_pl) + do k = 1, self%nenc + call io_write_frame_encounter(LUN, self%t(k), & + pl%id(ind1(k)), encbody%id(ind2(k)), & + pl%Gmass(ind1(k)), encbody%Gmass(ind2(k)), & + pl%radius(ind1(k)), encbody%radius(ind2(k)), & + self%x1(:,k), self%x2(:,k), & + self%v1(:,k), self%v2(:,k)) + end do + class is (swiftest_tp) + do k = 1, self%nenc + call io_write_frame_encounter(LUN, self%t(k), & + pl%id(ind1(k)), encbody%id(ind2(k)), & + pl%Gmass(ind1(k)), 0.0_DP, & + pl%radius(ind1(k)), 0.0_DP, & + self%x1(:,k), self%x2(:,k), & + self%v1(:,k), self%v2(:,k)) + end do + end select + end associate + + close(unit = LUN, iostat = ierr) if (ierr /= 0) then write(*, *) "Swiftest Error:" write(*, *) " Unable to close binary encounter file" call util_exit(FAILURE) end if - return end subroutine io_write_encounter @@ -1429,6 +1443,39 @@ module subroutine io_write_frame_cb(self, iu, param) end subroutine io_write_frame_cb + module subroutine io_write_frame_encounter(iu, t, id1, id2, Gmass1, Gmass2, radius1, radius2, xh1, xh2, vh1, vh2) + !! author: David A. Minton + !! + !! Write a single frame of close encounter data to output binary files + !! + !! Adapted from David E. Kaufmann's Swifter routine: io_write_encounter.f90 + !! Adapted from Hal Levison's Swift routine io_write_encounter.f + implicit none + ! Arguments + integer(I4B), intent(in) :: iu !! Open file unit number + real(DP), intent(in) :: t !! Time of encounter + integer(I4B), intent(in) :: id1, id2 !! ids of the two encountering bodies + real(DP), intent(in) :: Gmass1, Gmass2 !! G*mass of the two encountering bodies + real(DP), intent(in) :: radius1, radius2 !! Radii of the two encountering bodies + real(DP), dimension(:), intent(in) :: xh1, xh2 !! Heliocentric position vectors of the two encountering bodies + real(DP), dimension(:), intent(in) :: vh1, vh2 !! Heliocentric velocity vectors of the two encountering bodies + ! Internals + integer(I4B) :: ierr + + write(iu, iostat = ierr) t + write(iu, iostat = ierr) id1, xh1(1), xh1(2), xh1(3), vh1(1), vh1(2), Gmass1, radius1 + write(iu, iostat = ierr) id2, xh2(1), xh2(2), xh2(3), vh2(1), vh2(2), Gmass2, radius2 + + if (ierr /= 0) then + write(*, *) "Swiftest Error:" + write(*, *) " Unable to write binary file record for encounter" + call util_exit(FAILURE) + end if + + return + end subroutine + + module subroutine io_write_frame_system(self, iu, param) !! author: The Purdue Swiftest Team - David A. Minton, Carlisle A. Wishard, Jennifer L.L. Pouplin, and Jacob R. Elliott !! diff --git a/src/kick/kick.f90 b/src/kick/kick.f90 index 3945a91d0..ba3226eb7 100644 --- a/src/kick/kick.f90 +++ b/src/kick/kick.f90 @@ -15,19 +15,25 @@ module pure subroutine kick_getacch_int_pl(self) ! Internals integer(I4B) :: k real(DP) :: rji2, irij3, faci, facj - real(DP), dimension(NDIM) :: dx + real(DP) :: dx, dy, dz associate(pl => self, npl => self%nbody, nplpl => self%nplpl) do k = 1, nplpl associate(i => pl%k_plpl(1, k), j => pl%k_plpl(2, k)) if (pl%lmask(i) .and. pl%lmask(j)) then - dx(:) = pl%xh(:, j) - pl%xh(:, i) - rji2 = dot_product(dx(:), dx(:)) + dx = pl%xh(1, j) - pl%xh(1, i) + dy = pl%xh(2, j) - pl%xh(2, i) + dz = pl%xh(3, j) - pl%xh(3, i) + rji2 = dx**2 + dy**2 + dz**2 irij3 = 1.0_DP / (rji2 * sqrt(rji2)) faci = pl%Gmass(i) * irij3 facj = pl%Gmass(j) * irij3 - pl%ah(:, i) = pl%ah(:, i) + facj * dx(:) - pl%ah(:, j) = pl%ah(:, j) - faci * dx(:) + pl%ah(1, i) = pl%ah(1, i) + facj * dx + pl%ah(2, i) = pl%ah(2, i) + facj * dy + pl%ah(3, i) = pl%ah(3, i) + facj * dz + pl%ah(1, j) = pl%ah(1, j) - faci * dx + pl%ah(2, j) = pl%ah(2, j) - faci * dy + pl%ah(3, j) = pl%ah(3, j) - faci * dz end if end associate end do diff --git a/src/modules/rmvs_classes.f90 b/src/modules/rmvs_classes.f90 index 6ffb7ba1b..0b837ec0a 100644 --- a/src/modules/rmvs_classes.f90 +++ b/src/modules/rmvs_classes.f90 @@ -125,6 +125,14 @@ module function rmvs_encounter_check_tp(self, system, dt) result(lencounter) logical :: lencounter !! Returns true if there is at least one close encounter end function rmvs_encounter_check_tp + module subroutine rmvs_io_write_encounter(t, id1, id2, Gmass1, Gmass2, radius1, radius2, xh1, xh2, vh1, vh2, enc_out) + implicit none + integer(I4B), intent(in) :: id1, id2 + real(DP), intent(in) :: t, Gmass1, Gmass2, radius1, radius2 + real(DP), dimension(:), intent(in) :: xh1, xh2, vh1, vh2 + character(*), intent(in) :: enc_out + end subroutine rmvs_io_write_encounter + module subroutine rmvs_kick_getacch_tp(self, system, param, t, lbeg) use swiftest_classes, only : swiftest_nbody_system, swiftest_parameters implicit none diff --git a/src/modules/swiftest_classes.f90 b/src/modules/swiftest_classes.f90 index 2e4bff8a2..96c35675b 100644 --- a/src/modules/swiftest_classes.f90 +++ b/src/modules/swiftest_classes.f90 @@ -321,11 +321,13 @@ module swiftest_classes real(DP), dimension(:,:), allocatable :: x2 !! the position of body 2 in the encounter real(DP), dimension(:,:), allocatable :: v1 !! the velocity of body 1 in the encounter real(DP), dimension(:,:), allocatable :: v2 !! the velocity of body 2 in the encounter + real(DP), dimension(:), allocatable :: t !! Time of encounter contains procedure :: setup => setup_encounter !! A constructor that sets the number of encounters and allocates and initializes all arrays procedure :: copy => util_copy_encounter !! Copies elements from the source encounter list into self. procedure :: spill => util_spill_encounter !! "Spills" bodies from one object to another depending on the results of a mask (uses the PACK intrinsic) procedure :: resize => util_resize_encounter !! Checks the current size of the encounter list against the required size and extends it by a factor of 2 more than requested if it is too small. + procedure :: write => io_write_encounter !! Write close encounter data to output binary file end type swiftest_encounter abstract interface @@ -656,13 +658,12 @@ module subroutine io_toupper(string) character(*), intent(inout) :: string !! String to make upper case end subroutine io_toupper - module subroutine io_write_encounter(t, name1, name2, mass1, mass2, radius1, radius2, & - xh1, xh2, vh1, vh2, enc_out, out_type) + module subroutine io_write_encounter(self, pl, encbody, param) implicit none - integer(I4B), intent(in) :: name1, name2 - real(DP), intent(in) :: t, mass1, mass2, radius1, radius2 - real(DP), dimension(:), intent(in) :: xh1, xh2, vh1, vh2 - character(*), intent(in) :: enc_out, out_type + class(swiftest_encounter), intent(in) :: self !! Swiftest encounter list object + class(swiftest_pl), intent(in) :: pl !! Swiftest massive body object + class(swiftest_body), intent(in) :: encbody !! Encountering body - Swiftest generic body object (pl or tp) + class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters end subroutine io_write_encounter module subroutine io_write_frame_body(self, iu, param) @@ -679,6 +680,17 @@ module subroutine io_write_frame_cb(self, iu, param) class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters end subroutine io_write_frame_cb + module subroutine io_write_frame_encounter(iu, t, id1, id2, Gmass1, Gmass2, radius1, radius2, xh1, xh2, vh1, vh2) + implicit none + integer(I4B), intent(in) :: iu !! Open file unit number + real(DP), intent(in) :: t !! Time of encounter + integer(I4B), intent(in) :: id1, id2 !! ids of the two encountering bodies + real(DP), intent(in) :: Gmass1, Gmass2 !! G*mass of the two encountering bodies + real(DP), intent(in) :: radius1, radius2 !! Radii of the two encountering bodies + real(DP), dimension(:), intent(in) :: xh1, xh2 !! Heliocentric position vectors of the two encountering bodies + real(DP), dimension(:), intent(in) :: vh1, vh2 !! Heliocentric velocity vectors of the two encountering bodies + end subroutine io_write_frame_encounter + module subroutine io_write_frame_system(self, iu, param) implicit none class(swiftest_nbody_system), intent(in) :: self !! Swiftest system object diff --git a/src/modules/swiftest_globals.f90 b/src/modules/swiftest_globals.f90 index 5ec55f6c6..6cacac789 100644 --- a/src/modules/swiftest_globals.f90 +++ b/src/modules/swiftest_globals.f90 @@ -107,9 +107,6 @@ module swiftest_globals character(*), dimension(2), parameter :: DUMP_PARAM_FILE = ['dump_param1.dat', 'dump_param2.dat'] !> Default file names that can be changed by the user in the parameters file - character(*), parameter :: ENC_OUTFILE = 'encounter.out' - character(*), parameter :: DISCARD_FILE = 'discard.out' - character(*), parameter :: ENERGY_FILE = 'energy.out' character(*), parameter :: CB_INFILE = 'cb.in' character(*), parameter :: PL_INFILE = 'pl.in' character(*), parameter :: TP_INFILE = 'tp.in' diff --git a/src/modules/symba_classes.f90 b/src/modules/symba_classes.f90 index 4628202f8..6a69adcc7 100644 --- a/src/modules/symba_classes.f90 +++ b/src/modules/symba_classes.f90 @@ -149,7 +149,7 @@ module symba_classes !> SyMBA class for tracking pl-pl close encounters in a step type, extends(symba_pltpenc) :: symba_plplenc contains - procedure :: scrub_non_collision => symba_collision_encounter_scrub !! Processes the pl-pl encounter list remove only those encounters that led to a collision + procedure :: extract_collisions => symba_collision_encounter_extract_collisions !! Processes the pl-pl encounter list remove only those encounters that led to a collision procedure :: resolve_fragmentations => symba_collision_resolve_fragmentations !! Process list of collisions, determine the collisional regime, and then create fragments procedure :: resolve_mergers => symba_collision_resolve_mergers !! Process list of collisions and merge colliding bodies together end type symba_plplenc @@ -161,6 +161,7 @@ module symba_classes class(symba_merger), allocatable :: pl_adds !! List of added bodies in mergers or collisions class(symba_pltpenc), allocatable :: pltpenc_list !! List of massive body-test particle encounters in a single step class(symba_plplenc), allocatable :: plplenc_list !! List of massive body-massive body encounters in a single step + class(symba_plplenc), allocatable :: plplcollision_list !! List of massive body-massive body collisions in a single step integer(I4B) :: irec !! System recursion level contains procedure :: write_discard => symba_io_write_discard !! Write out information about discarded and merged planets and test particles in SyMBA @@ -184,11 +185,11 @@ module subroutine symba_collision_check_pltpenc(self, system, param, t, dt, irec integer(I4B), intent(in) :: irec !! Current recursion level end subroutine symba_collision_check_pltpenc - module subroutine symba_collision_encounter_scrub(self, system, param) + module subroutine symba_collision_encounter_extract_collisions(self, system, param) implicit none class(symba_plplenc), intent(inout) :: self !! SyMBA pl-pl encounter list class(symba_nbody_system), intent(inout) :: system !! SyMBA nbody system object - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameterss + class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters end subroutine module subroutine symba_collision_make_family_pl(self,idx) diff --git a/src/rmvs/rmvs_io.f90 b/src/rmvs/rmvs_io.f90 new file mode 100644 index 000000000..b85ce2202 --- /dev/null +++ b/src/rmvs/rmvs_io.f90 @@ -0,0 +1,48 @@ +submodule (rmvs_classes) s_rmvs_io + use swiftest +contains + + module subroutine rmvs_io_write_encounter(t, id1, id2, Gmass1, Gmass2, radius1, radius2, & + xh1, xh2, vh1, vh2, enc_out) + !! author: David A. Minton + !! + !! Write close encounter data from RMVS to output binary files + !! There is no direct file output from this subroutine + !! + !! Adapted from David E. Kaufmann's Swifter routine: io_write_encounter.f90 + !! Adapted from Hal Levison's Swift routine io_write_encounter.f + implicit none + ! Arguments + integer(I4B), intent(in) :: id1, id2 + real(DP), intent(in) :: t, Gmass1, Gmass2, radius1, radius2 + real(DP), dimension(:), intent(in) :: xh1, xh2, vh1, vh2 + character(*), intent(in) :: enc_out + ! Internals + logical , save :: lfirst = .true. + integer(I4B), parameter :: LUN = 30 + integer(I4B) :: ierr + + if (enc_out == "") return + + open(unit = LUN, file = enc_out, status = 'OLD', position = 'APPEND', form = 'UNFORMATTED', iostat = ierr) + if ((ierr /= 0) .and. lfirst) then + open(unit = LUN, file = enc_out, status = 'NEW', form = 'UNFORMATTED', iostat = ierr) + end if + if (ierr /= 0) then + write(*, *) "Swiftest Error:" + write(*, *) " Unable to open binary encounter file" + call util_exit(FAILURE) + end if + lfirst = .false. + call io_write_frame_encounter(LUN, t, id1, id2, Gmass1, Gmass2, radius1, radius2, xh1, xh2, vh1, vh2) + close(unit = LUN, iostat = ierr) + if (ierr /= 0) then + write(*, *) "Swiftest Error:" + write(*, *) " Unable to close binary encounter file" + call util_exit(FAILURE) + end if + + return + end subroutine rmvs_io_write_encounter + +end submodule s_rmvs_io \ No newline at end of file diff --git a/src/rmvs/rmvs_step.f90 b/src/rmvs/rmvs_step.f90 index 0385aeecc..c801ee16f 100644 --- a/src/rmvs/rmvs_step.f90 +++ b/src/rmvs/rmvs_step.f90 @@ -547,8 +547,7 @@ subroutine rmvs_peri_tp(tp, pl, t, dt, lfirst, inner_index, ipleP, param) id2 = tp%id(i) xh2(:) = xpc(:, i) + xh1(:) vh2(:) = xpc(:, i) + vh1(:) - call io_write_encounter(t, id1, id2, mu, 0.0_DP, rpl, 0.0_DP, xh1(:), xh2(:), vh1(:), vh2(:), & - param%enc_out, param%out_type) + call rmvs_io_write_encounter(t, id1, id2, mu, 0.0_DP, rpl, 0.0_DP, xh1(:), xh2(:), vh1(:), vh2(:), param%enc_out) end if if (tp%lperi(i)) then if (peri < tp%peri(i)) then diff --git a/src/setup/setup.f90 b/src/setup/setup.f90 index 8f96c48a1..726da07fc 100644 --- a/src/setup/setup.f90 +++ b/src/setup/setup.f90 @@ -56,8 +56,9 @@ module subroutine setup_construct_system(system, param) allocate(symba_tp :: system%tp_discards) allocate(symba_merger :: system%pl_adds) allocate(symba_merger :: system%pl_discards) - allocate(symba_plplenc :: system%plplenc_list) allocate(symba_pltpenc :: system%pltpenc_list) + allocate(symba_plplenc :: system%plplenc_list) + allocate(symba_plplenc :: system%plplcollision_list) end select case (RINGMOONS) write(*,*) 'RINGMOONS-SyMBA integrator not yet enabled' @@ -91,6 +92,7 @@ module subroutine setup_encounter(self, n) if (allocated(self%x2)) deallocate(self%x2) if (allocated(self%v1)) deallocate(self%v1) if (allocated(self%v2)) deallocate(self%v2) + if (allocated(self%t)) deallocate(self%t) allocate(self%lvdotr(n)) allocate(self%status(n)) @@ -100,6 +102,7 @@ module subroutine setup_encounter(self, n) allocate(self%x2(NDIM,n)) allocate(self%v1(NDIM,n)) allocate(self%v2(NDIM,n)) + allocate(self%t(n)) self%lvdotr(:) = .false. self%status(:) = INACTIVE @@ -109,6 +112,7 @@ module subroutine setup_encounter(self, n) self%x2(:,:) = 0.0_DP self%v1(:,:) = 0.0_DP self%v2(:,:) = 0.0_DP + self%t(:) = 0.0_DP return end subroutine setup_encounter diff --git a/src/symba/symba_collision.f90 b/src/symba/symba_collision.f90 index 1910411b9..b0e588300 100644 --- a/src/symba/symba_collision.f90 +++ b/src/symba/symba_collision.f90 @@ -66,32 +66,33 @@ module subroutine symba_collision_check_pltpenc(self, system, param, t, dt, irec end do end if - if (any(lcollision(:))) then - do k = 1, nenc - if (.not.lcollision(k)) cycle - self%status(k) = COLLISION - self%x1(:,k) = pl%xh(:,ind1(k)) - self%v1(:,k) = pl%vb(:,ind1(k)) - if (isplpl) then - self%x2(:,k) = pl%xh(:,ind2(k)) - self%v2(:,k) = pl%vb(:,ind2(k)) - + do k = 1, nenc + if (lcollision(k)) self%status(k) = COLLISION + self%t(k) = t + self%x1(:,k) = pl%xh(:,ind1(k)) + self%v1(:,k) = pl%vb(:,ind1(k)) - system%cb%vb(:) + if (isplpl) then + self%x2(:,k) = pl%xh(:,ind2(k)) + self%v2(:,k) = pl%vb(:,ind2(k)) - system%cb%vb(:) + if (lcollision(k)) then ! Check to see if either of these bodies has been involved with a collision before, and if so, make this a collisional family if (pl%lcollision(ind1(k)) .or. pl%lcollision(ind2(k))) call pl%make_family([ind1(k),ind2(k)]) - + ! Set the collision flag for these to bodies to true in case they become involved in another collision later in the step pl%lcollision([ind1(k), ind2(k)]) = .true. pl%ldiscard([ind1(k), ind2(k)]) = .true. pl%status([ind1(k), ind2(k)]) = COLLISION - else - self%x2(:,k) = tp%xh(:,ind2(k)) - self%v2(:,k) = tp%vb(:,ind2(k)) + end if + else + self%x2(:,k) = tp%xh(:,ind2(k)) + self%v2(:,k) = tp%vb(:,ind2(k)) - system%cb%vb(:) + if (lcollision(k)) then tp%status(ind2(k)) = DISCARDED_PLR tp%ldiscard(ind2(k)) = .true. write(*,*) 'Test particle ',tp%id(ind2(k)), ' collided with massive body ',pl%id(ind1(k)), ' at time ',t end if - end do - end if + end if + end do end associate end select end select @@ -267,7 +268,7 @@ function symba_collision_consolidate_familes(pl, param, idx_parent, family, x, v end function symba_collision_consolidate_familes - module subroutine symba_collision_encounter_scrub(self, system, param) + module subroutine symba_collision_encounter_extract_collisions(self, system, param) !! author: David A. Minton !! !! Processes the pl-pl encounter list remove only those encounters that led to a collision @@ -283,56 +284,55 @@ module subroutine symba_collision_encounter_scrub(self, system, param) integer(I4B), dimension(:), pointer :: plparent integer(I4B), dimension(:), allocatable :: collision_idx, unique_parent_idx integer(I4B) :: i, index_coll, ncollisions, nunique_parent - type(symba_plplenc) :: plplenc_noncollision select type (pl => system%pl) class is (symba_pl) associate(plplenc_list => self, nplplenc => self%nenc, idx1 => self%index1, idx2 => self%index2, plparent => pl%kin%parent) lplpl_collision(:) = plplenc_list%status(1:nplplenc) == COLLISION - if (any(lplpl_collision)) then ! Collisions have been detected in this step. So we need to determine which of them are between unique bodies. - - ! Get the subset of pl-pl encounters that lead to a collision - ncollisions = count(lplpl_collision(:)) - allocate(collision_idx(ncollisions)) - collision_idx = pack([(i, i=1, nplplenc)], lplpl_collision) - - ! Get the subset of collisions that involve a unique pair of parents - allocate(lplpl_unique_parent(ncollisions)) - - lplpl_unique_parent(:) = plparent(idx1(collision_idx(:))) /= plparent(idx2(collision_idx(:))) - nunique_parent = count(lplpl_unique_parent(:)) - allocate(unique_parent_idx(nunique_parent)) - unique_parent_idx = pack(collision_idx(:), lplpl_unique_parent(:)) - - ! Scrub all pl-pl collisions involving unique pairs of parents, which will remove all duplicates and leave behind - ! all pairs that have themselves as parents but are not part of the unique parent list. This can hapepn in rare cases - ! due to restructuring of parent/child relationships when there are large numbers of multi-body collisions in a single - ! step - lplpl_unique_parent(:) = .true. - do index_coll = 1, ncollisions - associate(ip1 => plparent(idx1(collision_idx(index_coll))), ip2 => plparent(idx2(collision_idx(index_coll)))) - lplpl_unique_parent(:) = .not. ( any(plparent(idx1(unique_parent_idx(:))) == ip1) .or. & - any(plparent(idx2(unique_parent_idx(:))) == ip1) .or. & - any(plparent(idx1(unique_parent_idx(:))) == ip2) .or. & - any(plparent(idx2(unique_parent_idx(:))) == ip2) ) - end associate - end do + if (.not.any(lplpl_collision)) return + ! Collisions have been detected in this step. So we need to determine which of them are between unique bodies. + + ! Get the subset of pl-pl encounters that lead to a collision + ncollisions = count(lplpl_collision(:)) + allocate(collision_idx(ncollisions)) + collision_idx = pack([(i, i=1, nplplenc)], lplpl_collision) + + ! Get the subset of collisions that involve a unique pair of parents + allocate(lplpl_unique_parent(ncollisions)) + + lplpl_unique_parent(:) = plparent(idx1(collision_idx(:))) /= plparent(idx2(collision_idx(:))) + nunique_parent = count(lplpl_unique_parent(:)) + allocate(unique_parent_idx(nunique_parent)) + unique_parent_idx = pack(collision_idx(:), lplpl_unique_parent(:)) + + ! Scrub all pl-pl collisions involving unique pairs of parents, which will remove all duplicates and leave behind + ! all pairs that have themselves as parents but are not part of the unique parent list. This can hapepn in rare cases + ! due to restructuring of parent/child relationships when there are large numbers of multi-body collisions in a single + ! step + lplpl_unique_parent(:) = .true. + do index_coll = 1, ncollisions + associate(ip1 => plparent(idx1(collision_idx(index_coll))), ip2 => plparent(idx2(collision_idx(index_coll)))) + lplpl_unique_parent(:) = .not. ( any(plparent(idx1(unique_parent_idx(:))) == ip1) .or. & + any(plparent(idx2(unique_parent_idx(:))) == ip1) .or. & + any(plparent(idx1(unique_parent_idx(:))) == ip2) .or. & + any(plparent(idx2(unique_parent_idx(:))) == ip2) ) + end associate + end do - ! Reassemble collision index list to include only those containing the unique pairs of parents, plus all the non-unique pairs that don't - ! contain a parent body on the unique parent list. - ncollisions = nunique_parent + count(lplpl_unique_parent) - collision_idx = [unique_parent_idx(:), pack(collision_idx(:), lplpl_unique_parent(:))] + ! Reassemble collision index list to include only those containing the unique pairs of parents, plus all the non-unique pairs that don't + ! contain a parent body on the unique parent list. + ncollisions = nunique_parent + count(lplpl_unique_parent) + collision_idx = [unique_parent_idx(:), pack(collision_idx(:), lplpl_unique_parent(:))] - ! Create a mask that contains only the pl-pl encounters that did not result in a collision, and then discard them - lplpl_collision(:) = .false. - lplpl_collision(collision_idx(:)) = .true. - end if - call plplenc_list%spill(plplenc_noncollision, .not.lplpl_collision, ldestructive=.true.) ! Remove any encounters that are not collisions from the list. + ! Create a mask that contains only the pl-pl encounters that did not result in a collision, and then discard them + lplpl_collision(:) = .false. + lplpl_collision(collision_idx(:)) = .true. + call plplenc_list%spill(system%plplcollision_list, lplpl_collision, ldestructive=.true.) ! Extract any encounters that are not collisions from the list. end associate end select return - end subroutine symba_collision_encounter_scrub + end subroutine symba_collision_encounter_extract_collisions module subroutine symba_collision_make_family_pl(self, idx) diff --git a/src/symba/symba_discard.f90 b/src/symba/symba_discard.f90 index c154b4290..499c527cc 100644 --- a/src/symba/symba_discard.f90 +++ b/src/symba/symba_discard.f90 @@ -291,21 +291,22 @@ module subroutine symba_discard_pl(self, system, param) class is (symba_nbody_system) select type(param) class is (symba_parameters) - associate(pl => self, plplenc_list => system%plplenc_list) + associate(pl => self, plplenc_list => system%plplenc_list, plplcollision_list => system%plplcollision_list) call pl%h2b(system%cb) ! First deal with the non pl-pl collisions call symba_discard_nonplpl(self, system, param) - ! Scrub the pl-pl encounter list of any encounters that did not lead to a collision - call plplenc_list%scrub_non_collision(system, param) + ! Extract the pl-pl encounter list and return the plplcollision_list + call plplenc_list%extract_collisions(system, param) + call plplenc_list%write(pl, pl, param) - if ((plplenc_list%nenc > 0) .and. any(pl%lcollision(:))) then + if ((plplcollision_list%nenc > 0) .and. any(pl%lcollision(:))) then write(*, *) "Collision between massive bodies detected at time t = ",param%t if (param%lfragmentation) then - call plplenc_list%resolve_fragmentations(system, param) + call plplcollision_list%resolve_fragmentations(system, param) else - call plplenc_list%resolve_mergers(system, param) + call plplcollision_list%resolve_mergers(system, param) end if end if diff --git a/src/symba/symba_fragmentation.f90 b/src/symba/symba_fragmentation.f90 index b36c54e9a..335ab9e66 100644 --- a/src/symba/symba_fragmentation.f90 +++ b/src/symba/symba_fragmentation.f90 @@ -30,7 +30,9 @@ module function symba_fragmentation_casedisruption(system, param, family, x, v, real(DP), dimension(:), allocatable :: m_frag, rad_frag integer(I4B), dimension(:), allocatable :: id_frag logical :: lfailure - + + write(*, '("Disruption between bodies ",I8,99(:,",",I8))') system%pl%id(family(:)) + ! Collisional fragments will be uniformly distributed around the pre-impact barycenter nfrag = NFRAG_DISRUPT allocate(m_frag(nfrag)) @@ -117,6 +119,8 @@ module function symba_fragmentation_casehitandrun(system, param, family, x, v, m logical :: lpure logical, dimension(system%pl%nbody) :: lmask + write(*, '("Hit and run between bodies ",I8,99(:,",",I8))') system%pl%id(family(:)) + mtot = sum(mass(:)) xcom(:) = (mass(1) * x(:,1) + mass(2) * x(:,2)) / mtot vcom(:) = (mass(1) * v(:,1) + mass(2) * v(:,2)) / mtot @@ -175,8 +179,14 @@ module function symba_fragmentation_casehitandrun(system, param, family, x, v, m write(*,'("Generating ",I2.0," fragments")') nfrag end if end if - if (lpure) then + if (lpure) then ! Reset these bodies back to being active so that nothing further is done to them status = ACTIVE + select type(pl => system%pl) + class is (symba_pl) + pl%status(family(:)) = status + pl%ldiscard(family(:)) = .false. + pl%lcollision(family(:)) = .false. + end select else status = HIT_AND_RUN call symba_fragmentation_mergeaddsub(system, param, family, id_frag, Ip_frag, m_frag, rad_frag, xb_frag, vb_frag, rot_frag, status) @@ -217,7 +227,7 @@ module function symba_fragmentation_casemerge(system, param, family, x, v, mass, select type(pl => system%pl) class is (symba_pl) - write(*, '("Merging bodies ",99(I8,",",:))') pl%id(family(:)) + write(*, '("Merging bodies ",I8,99(:,",",I8))') pl%id(family(:)) ibiggest = maxloc(pl%Gmass(family(:)), dim=1) id_frag(1) = pl%id(family(ibiggest)) @@ -304,6 +314,8 @@ module function symba_fragmentation_casesupercatastrophic(system, param, family, logical :: lfailure logical, dimension(system%pl%nbody) :: lmask + write(*, '("Supercatastrophic disruption between bodies ",I8,99(:,",",I8))') system%pl%id(family(:)) + ! Collisional fragments will be uniformly distributed around the pre-impact barycenter nfrag = NFRAG_SUPERCAT allocate(m_frag(nfrag)) diff --git a/src/symba/symba_io.f90 b/src/symba/symba_io.f90 index 97e5798c2..73027c351 100644 --- a/src/symba/symba_io.f90 +++ b/src/symba/symba_io.f90 @@ -344,6 +344,8 @@ module subroutine symba_io_write_discard(self, param) character(*), parameter :: PLNAMEFMT = '(I8, 2(1X, E23.16))' class(swiftest_body), allocatable :: pltemp + if (param%discard_out == "") return + associate(pl => self%pl, npl => self%pl%nbody, pl_adds => self%pl_adds) if (self%tp_discards%nbody > 0) call io_write_discard(self, param) select type(pl_discards => self%pl_discards) diff --git a/src/symba/symba_kick.f90 b/src/symba/symba_kick.f90 index 8625b3d81..7a98c2f69 100644 --- a/src/symba/symba_kick.f90 +++ b/src/symba/symba_kick.f90 @@ -128,18 +128,21 @@ module subroutine symba_kick_pltpenc(self, system, dt, irec, sgn) if (tp%nbody > 0) tp%lmask(:) = tp%status(:) /= INACTIVE irm1 = irec - 1 + if (sgn < 0) then irecl = irec - 1 else irecl = irec end if + + if (isplpl) then + pl%ah(:,ind1(1:self%nenc)) = 0.0_DP + pl%ah(:,ind2(1:self%nenc)) = 0.0_DP + else + tp%ah(:,ind2(1:self%nenc)) = 0.0_DP + end if + do k = 1, self%nenc - if (isplpl) then - pl%ah(:,ind1(k)) = 0.0_DP - pl%ah(:,ind2(k)) = 0.0_DP - else - tp%ah(:,ind2(k)) = 0.0_DP - end if if (isplpl) then lgoodlevel = (pl%levelg(ind1(k)) >= irm1) .and. (pl%levelg(ind2(k)) >= irm1) else @@ -182,7 +185,7 @@ module subroutine symba_kick_pltpenc(self, system, dt, irec, sgn) pl%vb(:,ind1(k)) = pl%vb(:,ind1(k)) + sgn * dt * pl%ah(:,ind1(k)) pl%vb(:,ind2(k)) = pl%vb(:,ind2(k)) + sgn * dt * pl%ah(:,ind2(k)) pl%ah(:,ind1(k)) = 0.0_DP - pl%ah(:,ind1(k)) = 0.0_DP + pl%ah(:,ind2(k)) = 0.0_DP end do else do k = 1, self%nenc diff --git a/src/symba/symba_setup.f90 b/src/symba/symba_setup.f90 index e06fb20b5..e4644f25e 100644 --- a/src/symba/symba_setup.f90 +++ b/src/symba/symba_setup.f90 @@ -76,6 +76,7 @@ module subroutine symba_setup_initialize_system(self, param) call whm_setup_initialize_system(system, param) call system%pltpenc_list%setup(0) call system%plplenc_list%setup(0) + call system%plplcollision_list%setup(0) select type(pl => system%pl) class is (symba_pl) call pl%sort("mass", ascending=.false.) diff --git a/src/symba/symba_step.f90 b/src/symba/symba_step.f90 index 7065625b4..847d6e4ee 100644 --- a/src/symba/symba_step.f90 +++ b/src/symba/symba_step.f90 @@ -237,7 +237,7 @@ module subroutine symba_step_reset_system(self) ! Internals integer(I4B) :: i - associate(system => self, pltpenc_list => self%pltpenc_list, plplenc_list => self%plplenc_list, pl_adds => self%pl_adds, pl_discards => self%pl_discards) + associate(system => self) select type(pl => system%pl) class is (symba_pl) select type(tp => system%tp) @@ -255,18 +255,19 @@ module subroutine symba_step_reset_system(self) pl%levelm(:) = 0 pl%lencounter = .false. pl%lcollision = .false. - plplenc_list%nenc = 0 + system%plplenc_list%nenc = 0 + system%plplcollision_list%nenc = 0 end if if (tp%nbody > 0) then tp%nplenc(:) = 0 tp%levelg(:) = 0 tp%levelm(:) = 0 - pltpenc_list%nenc = 0 + system%pltpenc_list%nenc = 0 end if - call pl_adds%resize(0) - call pl_discards%resize(0) + call system%pl_adds%resize(0) + call system%pl_discards%resize(0) end select end select end associate diff --git a/src/symba/symba_util.f90 b/src/symba/symba_util.f90 index 4c0f256e3..1335272fe 100644 --- a/src/symba/symba_util.f90 +++ b/src/symba/symba_util.f90 @@ -704,14 +704,25 @@ module subroutine symba_util_spill_arr_info(keeps, discards, lspill_list, ldestr type(symba_particle_info), dimension(:), allocatable, intent(inout) :: discards !! Array of discards logical, dimension(:), intent(in) :: lspill_list !! Logical array of bodies to spill into the discardss logical, intent(in) :: ldestructive !! Logical flag indicating whether or not this operation should alter the keeps array or not + ! Internals + integer(I4B) :: nspill, nkeep, nlist + + nkeep = count(.not.lspill_list(:)) + nspill = count(lspill_list(:)) + nlist = size(lspill_list(:)) + + if (.not.allocated(keeps) .or. nspill == 0) return + if (.not.allocated(discards)) then + allocate(discards(nspill)) + else if (size(discards) /= nspill) then + deallocate(discards) + allocate(discards(nspill)) + end if - if (.not.allocated(keeps) .or. count(lspill_list(:)) == 0) return - if (.not.allocated(discards)) allocate(discards(count(lspill_list(:)))) - - discards(:) = pack(keeps(:), lspill_list(:)) + discards(:) = pack(keeps(1:nlist), lspill_list(1:nlist)) if (ldestructive) then - if (count(.not.lspill_list(:)) > 0) then - keeps(:) = pack(keeps(:), .not. lspill_list(:)) + if (nkeep > 0) then + keeps(:) = pack(keeps(1:nlist), .not. lspill_list(1:nlist)) else deallocate(keeps) end if @@ -732,14 +743,25 @@ module subroutine symba_util_spill_arr_kin(keeps, discards, lspill_list, ldestru type(symba_kinship), dimension(:), allocatable, intent(inout) :: discards !! Array of discards logical, dimension(:), intent(in) :: lspill_list !! Logical array of bodies to spill into the discardss logical, intent(in) :: ldestructive !! Logical flag indicating whether or not this operation should alter the keeps array or not + ! Internals + integer(I4B) :: nspill, nkeep, nlist + + nkeep = count(.not.lspill_list(:)) + nspill = count(lspill_list(:)) + nlist = size(lspill_list(:)) + + if (.not.allocated(keeps) .or. nspill == 0) return + if (.not.allocated(discards)) then + allocate(discards(nspill)) + else if (size(discards) /= nspill) then + deallocate(discards) + allocate(discards(nspill)) + end if - if (.not.allocated(keeps) .or. count(lspill_list(:)) == 0) return - if (.not.allocated(discards)) allocate(discards(count(lspill_list(:)))) - - discards(:) = pack(keeps(:), lspill_list(:)) + discards(:) = pack(keeps(1:nlist), lspill_list(1:nlist)) if (ldestructive) then - if (count(.not.lspill_list(:)) > 0) then - keeps(:) = pack(keeps(:), .not. lspill_list(:)) + if (nkeep > 0) then + keeps(:) = pack(keeps(1:nlist), .not. lspill_list(1:nlist)) else deallocate(keeps) end if @@ -806,7 +828,7 @@ module subroutine symba_util_spill_pltpenc(self, discards, lspill_list, ldestruc ! Internals integer(I4B) :: i - associate(keeps => self) + associate(keeps => self, nenc => self%nenc) select type(discards) class is (symba_pltpenc) call util_spill(keeps%level, discards%level, lspill_list, ldestructive) diff --git a/src/util/util_copy.f90 b/src/util/util_copy.f90 index f44777eec..87634a419 100644 --- a/src/util/util_copy.f90 +++ b/src/util/util_copy.f90 @@ -21,6 +21,7 @@ module subroutine util_copy_encounter(self, source) self%x2(:,1:n) = source%x2(:,1:n) self%v1(:,1:n) = source%v1(:,1:n) self%v2(:,1:n) = source%v2(:,1:n) + self%t(1:n) = source%t(1:n) end associate return diff --git a/src/util/util_rescale.f90 b/src/util/util_rescale.f90 index 62a9409ec..50ebc088c 100644 --- a/src/util/util_rescale.f90 +++ b/src/util/util_rescale.f90 @@ -20,9 +20,11 @@ module subroutine util_rescale_system(self, param, mscale, dscale, tscale) ! Calculate the G for the system units param%GU = GC / (param%DU2M**3 / (param%MU2KG * param%TU2S**2)) - ! Calculate the inverse speed of light in the system units - param%inv_c2 = einsteinC * param%TU2S / param%DU2M - param%inv_c2 = (param%inv_c2)**(-2) + if (param%lgr) then + ! Calculate the inverse speed of light in the system units + param%inv_c2 = einsteinC * param%TU2S / param%DU2M + param%inv_c2 = (param%inv_c2)**(-2) + end if vscale = dscale / tscale diff --git a/src/util/util_spill.f90 b/src/util/util_spill.f90 index 9c76ff5e9..a26fbfad7 100644 --- a/src/util/util_spill.f90 +++ b/src/util/util_spill.f90 @@ -13,14 +13,25 @@ module subroutine util_spill_arr_char_string(keeps, discards, lspill_list, ldest character(len=STRMAX), dimension(:), allocatable, intent(inout) :: discards !! Array of discards logical, dimension(:), intent(in) :: lspill_list !! Logical array of bodies to spill into the discards logical, intent(in) :: ldestructive !! Logical flag indicating whether or not this operation should alter the keeps array or not + ! Internals + integer(I4B) :: nspill, nkeep, nlist + + nkeep = count(.not.lspill_list(:)) + nspill = count(lspill_list(:)) + nlist = size(lspill_list(:)) + + if (.not.allocated(keeps) .or. nspill == 0) return + if (.not.allocated(discards)) then + allocate(discards(nspill)) + else if (size(discards) /= nspill) then + deallocate(discards) + allocate(discards(nspill)) + end if - if (.not.allocated(keeps) .or. count(lspill_list(:)) == 0) return - if (.not.allocated(discards)) allocate(discards(count(lspill_list(:)))) - - discards(:) = pack(keeps(:), lspill_list(:)) + discards(:) = pack(keeps(1:nlist), lspill_list(1:nlist)) if (ldestructive) then - if (count(.not.lspill_list(:)) > 0) then - keeps(:) = pack(keeps(:), .not. lspill_list(:)) + if (nkeep > 0) then + keeps(:) = pack(keeps(1:nlist), .not. lspill_list(1:nlist)) else deallocate(keeps) end if @@ -40,14 +51,25 @@ module subroutine util_spill_arr_DP(keeps, discards, lspill_list, ldestructive) real(DP), dimension(:), allocatable, intent(inout) :: discards !! Array of discards logical, dimension(:), intent(in) :: lspill_list !! Logical array of bodies to spill into the discardss logical, intent(in) :: ldestructive !! Logical flag indicating whether or not this operation should alter the keeps array or not + ! Internals + integer(I4B) :: nspill, nkeep, nlist + + nkeep = count(.not.lspill_list(:)) + nspill = count(lspill_list(:)) + nlist = size(lspill_list(:)) + + if (.not.allocated(keeps) .or. nspill == 0) return + if (.not.allocated(discards)) then + allocate(discards(nspill)) + else if (size(discards) /= nspill) then + deallocate(discards) + allocate(discards(nspill)) + end if - if (.not.allocated(keeps) .or. count(lspill_list(:)) == 0) return - if (.not.allocated(discards)) allocate(discards(count(lspill_list(:)))) - - discards(:) = pack(keeps(:), lspill_list(:)) + discards(:) = pack(keeps(1:nlist), lspill_list(1:nlist)) if (ldestructive) then - if (count(.not.lspill_list(:)) > 0) then - keeps(:) = pack(keeps(:), .not. lspill_list(:)) + if (nkeep > 0) then + keeps(:) = pack(keeps(1:nlist), .not. lspill_list(1:nlist)) else deallocate(keeps) end if @@ -68,19 +90,30 @@ module subroutine util_spill_arr_DPvec(keeps, discards, lspill_list, ldestructiv logical, dimension(:), intent(in) :: lspill_list !! Logical array of bodies to spill into the discards logical, intent(in) :: ldestructive !! Logical flag indicating whether or not this operation should alter the keeps array or not ! Internals - integer(I4B) :: i - - if (.not.allocated(keeps) .or. count(lspill_list(:)) == 0) return - if (.not.allocated(discards)) allocate(discards(NDIM, count(lspill_list(:)))) + integer(I4B) :: i, nspill, nkeep, nlist + + nkeep = count(.not.lspill_list(:)) + nspill = count(lspill_list(:)) + nlist = size(lspill_list(:)) + + if (.not.allocated(keeps) .or. nspill == 0) return + if (.not.allocated(discards)) then + allocate(discards(NDIM, nspill)) + else if (size(discards, dim=2) /= nspill) then + deallocate(discards) + allocate(discards(NDIM, nspill)) + end if do i = 1, NDIM - discards(i,:) = pack(keeps(i,:), lspill_list(:)) + discards(i,:) = pack(keeps(i,1:nlist), lspill_list(1:nlist)) end do if (ldestructive) then - if (count(.not.lspill_list(:)) > 0) then + if (nkeep > 0) then do i = 1, NDIM - keeps(i,:) = pack(keeps(i,:), .not. lspill_list(:)) + keeps(i,:) = pack(keeps(i,1:nlist), .not. lspill_list(1:nlist)) end do + else + deallocate(keeps) end if end if @@ -98,14 +131,25 @@ module subroutine util_spill_arr_I4B(keeps, discards, lspill_list, ldestructive) integer(I4B), dimension(:), allocatable, intent(inout) :: discards !! Array of discards logical, dimension(:), intent(in) :: lspill_list !! Logical array of bodies to spill into the discards logical, intent(in) :: ldestructive !! Logical flag indicating whether or not this operation should alter the keeps array or not + ! Internals + integer(I4B) :: nspill, nkeep, nlist + + nkeep = count(.not.lspill_list(:)) + nspill = count(lspill_list(:)) + nlist = size(lspill_list(:)) + + if (.not.allocated(keeps) .or. nspill == 0) return + if (.not.allocated(discards)) then + allocate(discards(nspill)) + else if (size(discards) /= nspill) then + deallocate(discards) + allocate(discards(nspill)) + end if - if (.not.allocated(keeps) .or. count(lspill_list(:)) == 0) return - if (.not.allocated(discards)) allocate(discards(count(lspill_list(:)))) - - discards(:) = pack(keeps(:), lspill_list(:)) + discards(:) = pack(keeps(1:nlist), lspill_list(1:nlist)) if (ldestructive) then - if (count(.not.lspill_list(:)) > 0) then - keeps(:) = pack(keeps(:), .not. lspill_list(:)) + if (nkeep > 0) then + keeps(:) = pack(keeps(1:nlist), .not. lspill_list(1:nlist)) else deallocate(keeps) end if @@ -125,14 +169,25 @@ module subroutine util_spill_arr_logical(keeps, discards, lspill_list, ldestruct logical, dimension(:), allocatable, intent(inout) :: discards !! Array of discards logical, dimension(:), intent(in) :: lspill_list !! Logical array of bodies to spill into the discards logical, intent(in) :: ldestructive !! Logical flag indicating whether or not this operation should alter the keeps array or no + ! Internals + integer(I4B) :: nspill, nkeep, nlist + + nkeep = count(.not.lspill_list(:)) + nspill = count(lspill_list(:)) + nlist = size(lspill_list(:)) + + if (.not.allocated(keeps) .or. nspill == 0) return + if (.not.allocated(discards)) then + allocate(discards(nspill)) + else if (size(discards) /= nspill) then + deallocate(discards) + allocate(discards(nspill)) + end if - if (.not.allocated(keeps) .or. count(lspill_list(:)) == 0) return - if (.not.allocated(discards)) allocate(discards(count(lspill_list(:)))) - - discards(:) = pack(keeps(:), lspill_list(:)) + discards(:) = pack(keeps(1:nlist), lspill_list(1:nlist)) if (ldestructive) then - if (count(.not.lspill_list(:)) > 0) then - keeps(:) = pack(keeps(:), .not. lspill_list(:)) + if (nkeep > 0) then + keeps(:) = pack(keeps(1:nlist), .not. lspill_list(1:nlist)) else deallocate(keeps) end if @@ -154,7 +209,6 @@ module subroutine util_spill_body(self, discards, lspill_list, ldestructive) logical, dimension(:), intent(in) :: lspill_list !! Logical array of bodies to spill into the discards logical, intent(in) :: ldestructive !! Logical flag indicating whether or not this operation should alter body by removing the discard list ! Internals - integer(I4B) :: i ! For each component, pack the discarded bodies into the discard object and do the inverse with the keeps !> Spill all the common components @@ -185,7 +239,6 @@ module subroutine util_spill_body(self, discards, lspill_list, ldestructive) discards%nbody = count(lspill_list(:)) keeps%nbody = keeps%nbody - discards%nbody if (keeps%nbody > size(keeps%status)) keeps%status(keeps%nbody+1:size(keeps%status)) = INACTIVE - end associate return @@ -201,11 +254,8 @@ module subroutine util_spill_encounter(self, discards, lspill_list, ldestructive class(swiftest_encounter), intent(inout) :: discards !! Discarded object logical, dimension(:), intent(in) :: lspill_list !! Logical array of bodies to spill into the discards logical, intent(in) :: ldestructive !! Logical flag indicating whether or not this operation should alter body by removing the discard list - ! Internals - integer(I4B) :: i associate(keeps => self) - call util_spill(keeps%lvdotr, discards%lvdotr, lspill_list, ldestructive) call util_spill(keeps%status, discards%status, lspill_list, ldestructive) call util_spill(keeps%index1, discards%index1, lspill_list, ldestructive) @@ -214,6 +264,7 @@ module subroutine util_spill_encounter(self, discards, lspill_list, ldestructive call util_spill(keeps%x2, discards%x2, lspill_list, ldestructive) call util_spill(keeps%v1, discards%v1, lspill_list, ldestructive) call util_spill(keeps%v2, discards%v2, lspill_list, ldestructive) + call util_spill(keeps%t, discards%t, lspill_list, ldestructive) ! This is the base class, so will be the last to be called in the cascade. ! Therefore we need to set the nenc values for both the keeps and discareds @@ -237,11 +288,8 @@ module subroutine util_spill_pl(self, discards, lspill_list, ldestructive) class(swiftest_body), intent(inout) :: discards !! Discarded object logical, dimension(:), intent(in) :: lspill_list !! Logical array of bodies to spill into the discards logical, intent(in) :: ldestructive !! Logical flag indicating whether or not this operation should alter body by removing the discard list - ! Internals - integer(I4B) :: i associate(keeps => self) - select type (discards) ! The standard requires us to select the type of both arguments in order to access all the components class is (swiftest_pl) !> Spill components specific to the massive body class