Updated the SFD code to pin the top of the SFD rather than the bottom…

…. This means that bodies smaller than the minimum can be made, but the minimum still sets the total number

examples/Fragmentation/Fragmentation_Movie.py

@@ -166,7 +166,7 @@ def encounter_combiner(sim):
     enc = enc.sel(time=tgood)
 
     # The following will combine the two datasets along the time dimension, sort the time dimension, and then fill in any time gaps with interpolation
-    ds = xr.combine_nested([data,enc],concat_dim='time').sortby("time").chunk(dict(time=-1)).interpolate_na(dim="time", method="akima")
+    ds = xr.combine_nested([data,enc],concat_dim='time').sortby("time").interpolate_na(dim="time", method="akima")
 
     # Rename the merged Target body so that their data can be combined
     tname=[n for n in ds['name'].data if names[0] in n]
@@ -359,7 +359,7 @@ def vec_props(self, c):
         minimum_fragment_gmass = 0.01 * body_Gmass[style][1] 
         gmtiny = 0.10 * body_Gmass[style][1] 
         sim.set_parameter(collision_model="fraggle", encounter_save="both", gmtiny=gmtiny, minimum_fragment_gmass=minimum_fragment_gmass, nfrag_reduction=nfrag_reduction[style])
-        sim.run(dt=5e-4, tstop=tstop[style], istep_out=1, dump_cadence=0, dask=False)
+        sim.run(dt=5e-4, tstop=tstop[style], istep_out=1, dump_cadence=0)
 
         print("Generating animation")
         anim = AnimatedScatter(sim,movie_filename,movie_titles[style],style,nskip=1)

src/fraggle/fraggle_util.f90

@@ -175,15 +175,16 @@ module subroutine fraggle_util_set_mass_dist(self, param)
 
             ! The mass will be distributed evenly in logspace between the second-largest remnant and the minimum mass
             ! Use Newton's method solver to get the logspace slope of the mass function
-            Mrat = mremaining / min_mfrag
+            Mrat = (mremaining + Mslr) / Mslr
             nrem = nfrag - 2
-            x0 = 1.0_DP + 100*epsilon(1.0_DP)
+            x0 = 1.0_DP - 100*epsilon(1.0_DP)
             x1 = Mrat**(1.0/nrem)
             do j = 1, MAXLOOP 
                y0 = Mrat - (1.0_DP - x0**nrem)/(1.0_DP - x0)
                y1 = Mrat - (1.0_DP - x1**nrem)/(1.0_DP - x1)
                if (y0*y1 < 0.0_DP) exit
                x1 = x1 * 1.6_DP
+               x0 = x0 / 1.6_DP
             end do
 
             ! Find the mass scaling factor with bisection
@@ -214,11 +215,11 @@ module subroutine fraggle_util_set_mass_dist(self, param)
 
             Mslr = impactors%mass_dist(iMslr)
             mass(iMslr) = Mslr
-            mfrag = min_mfrag
+            mfrag = Mslr 
             do i = iMrem,nfrag
+               mfrag = mfrag * mscale
                mass(i) = mfrag
                mremaining = mremaining - mfrag
-               mfrag = max(mfrag * mscale, min_mfrag)
             end do
 
             ! There may still be some small residual due to round-off error. If so, simply add it to the last bin of the mass distribution.