r_soviet_overpressure throws error when psi between 1.0-1.5

[lwillard]

e.g., s = r_soviet_overpressure(float(20000), 1.0, float(0), False, "kT", "m", "psi")

Error does not appear to be yield or HOB dependent, purely overpressure below a certain threshold.

File "D:\data\SIOP\ottawa.py", line 19, in
s = r_soviet_overpressure(float(20000), 1.2, float(0), False, "kT", "m", "psi")
File "C:\Users\xxxx\AppData\Local\Programs\Python\Python310\lib\site-packages\glasstone-0.0.1-py3.10.egg\glasstone\overpressure.py", line 869, in r_soviet_overpressure
File "C:\Users\xxxx\AppData\Local\Programs\Python\Python310\lib\site-packages\glasstone-0.0.1-py3.10.egg\glasstone\overpressure.py", line 826, in _rsovietmach
File "C:\Users\xxxx\AppData\Local\Programs\Python\Python310\lib\site-packages\glasstone-0.0.1-py3.10.egg\glasstone\overpressure.py", line 771, in _rsoviet_mach_sh7
glasstone.utilities.ValueOutsideGraphError: -1.0738201869691235

[GOFAI]

I know it's frustrating, but this is the intended behavior. The r_soviet_overpressure function is based upon digitized data from a chart in a 1987 Soviet military nuclear weapons effects manual. In order to avoid giving misleading results, I designed all these functions to throw the ValueOutsideGraphError when they are called with values that aren't included in the original Soviet chart.

One of the deficiencies of the original Soviet charts that overpressure function is based upon is that they do not extend into the very high overpressure range. Since you seem to be interested in very high overpressure (20,000psi) in a case without the thermal layer, I suggest using the Brode overpressure function to plot your case of interest to estimate the ranges at which various overpressures will occur. The Brode equation is supposed to work well into the high overpressure range for cases without the thermal layer.

[lwillard]

I'm actually not that interested in high overpressure ranges (over 20 PSI, things start to matter less and less). What I'm doing it taking David Teter's OPEN-RISOP laydown files ( https://github.com/davidteter/OPEN-RISOP), and using that to generate KML files with 1.0, 2.0, and 5.0 PSI blast rings, so I was using the r_soviet_overpressure function because it gives me the radius for a given overpressure. Except it doesn't for certain combinations of height of burst and low overpressure. If there's a better way of doing it using Glasstone, I would love to know, but as it stands now I'm looking at brute forcing it using curve_fit and overpressures calculated at specific distances using the Brode equation. But all this, from a (personal) requirement perspective, needs to be automated. I don't want to have to hand calculate overpressure distances every time yield or hob changes. Example output from earlier runs: [image: DFW Footprint.png]

…

On Fri, Apr 29, 2022 at 9:52 AM Edward Geist ***@***.***> wrote: I know it's frustrating, but this is the intended behavior. The r_soviet_overpressure function is based upon digitized data from a chart in a 1987 Soviet military nuclear weapons effects manual. In order to avoid giving misleading results, I designed all these functions to throw the ValueOutsideGraphError when they are called with values that aren't included in the original Soviet chart. One of the deficiencies of the original Soviet charts that overpressure function is based upon is that they do not extend into the very high overpressure range. Since you seem to be interested in very high overpressure (20,000psi) in a case without the thermal layer, I suggest using the Brode overpressure function to plot your case of interest to estimate the ranges at which various overpressures will occur. The Brode equation is supposed to work well into the high overpressure range for cases without the thermal layer. — Reply to this email directly, view it on GitHub <#7 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ABLDDKD27QWFJS3PVWJ2FDTVHPZRTANCNFSM5UWBT35A> . You are receiving this because you authored the thread.Message ID: ***@***.***>

[GOFAI]

How are you selecting heights-of-burst? Are you aiming to select an "optimal" HOB for each blast ring, and if so for what burst yields and detonation conditions (thermal precursors etc.)? Depending on quite what you need I may be able to suggest a more elegant solution.

[lwillard]

Sorry, just saw this buried in my email. I'm using the laydown file data in OPEN-RISOP by David Teeter ( https://github.com/davidteter/OPEN-RISOP). IIRC, he pre-selected HOB to maximize 5 PSI radius for airbursts depending on yield.

…

-Lane

On Sun, May 1, 2022 at 1:42 AM Edward Geist ***@***.***> wrote: How are you selecting heights-of-burst? Are you aiming to select an "optimal" HOB for each blast ring, and if so for what burst yields and detonation conditions (thermal precursors etc.)? Depending on quite what you need I may be able to suggest a more elegant solution. — Reply to this email directly, view it on GitHub <#7 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ABLDDKA76OJOYX33ZMSPDELVHYRVDANCNFSM5UWBT35A> . You are receiving this because you authored the thread.Message ID: ***@***.***>

[GOFAI]

So I've checked and the fixed_point function in scipy.optimize should do what you need, e.g.:

optimize.fixed_point(lambda r: brode_overpressure(teters_yield, r, teters_height, opunits='psi'), desired_overpressure)

[lwillard]

Mega delayed response, but I finally had time to work on this problem. So, optimize.fixed_point didn't work for me, but I ultimately got a working brute force solution. Definitely sub-optimal from a performance perspective, but seems to be accurate and it's fast enough for government work. Obviously, this is derivative work from (and uses) your glasstone library you developed and uploaded to github. Since I've cleaned up the code enough to where I wouldn't be horribly embarrassed if some of my co-workers saw it, I'll be uploading to github soon (need to create some test harnesses/examples for people if they are interested in using it). I'll shoot you over a link when it's uploaded (mainly to review that proper attribution was made in regards to your work). # use this class to calculate the altitude to detonate a nuke in order # to optimize the area under a particular PSI footprint radius class YieldAltitudeOptimizer: def __init__(self, max_radius=25000, max_height=5000, radius_interval= 100, height_interval=10): self.max_radius = max_radius self.max_height = max_height self.radius_interval = radius_interval self.height_interval = height_interval # distance/radius intervals self.t = np.arange(0.01, self.max_radius, radius_interval) # height of detonation intervals self.u = np.arange(0, self.max_height, height_interval) # for a given bomb yield and desired over-pressure, calculates and returns the optimal detonation height to maximize the radius of the given overpressure def optimize_for_overpressure(self, bomb_yield, overpressure): arrs = [] for height_of_burst in self.u: ovps = brode_overpressure(bomb_yield, self.t, height_of_burst, 'kT', dunits='m', opunits='psi') # this defines a function that calculates/interpolates the extent of the radius where an overpressure is found for a given yield and height of burst f = interpolate.interp1d(ovps, self.t, bounds_error=False, fill_value=0.0) # append to the array arrs.append(f(overpressure)) # find the max radius of the overpressure across all of the heights max_radius = max(arrs) # find the index of the max radius max_index = arrs.index(max_radius) # return the value for the index of the height array used to generate the max radius return self.u[max_index]

…

On Sun, May 22, 2022 at 1:05 AM Edward Geist ***@***.***> wrote: So I've checked and the fixed_point function in scipy.optimize <https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.fixed_point.html> should do what you need, e.g.: optimize.fixed_point(lambda r: brode_overpressure(teters_yield, r, teters_height, opunits='psi'), desired_overpressure) — Reply to this email directly, view it on GitHub <#7 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ABLDDKFMCY4UZH6QOVNQ6A3VLHFA3ANCNFSM5UWBT35A> . You are receiving this because you authored the thread.Message ID: ***@***.***>

[GOFAI]

Can you say more about why optimize.fixed_point didn't work for you? I checked just now using Python 2 and this worked:

>>> from glasstone.overpressure import brode_overpressure
>>> from scipy import optimize
>>> optimize.fixed_point(lambda r: brode_overpressure(750, r, 2000, opunits='psi'), 5.0)
array(35.0354849668831)
>>> optimize.fixed_point(lambda r: brode_overpressure(750, r, 200, opunits='psi'), 5.0)
array(530.7211266613074)

[lwillard]

Oh geez, I just realized I had conflated two different problems. One, I told you about, and one I didn't. Problem #2 (the "secret one"): given a weapon yield and a desired PSI footprint, what is the detonation height that will maximize the radius exposed to the given PSI (or greater)? That's the code snippet I excitedly sent you this morning and which probably confused you quite a bit since all of the context was missing. Problem #1 was: given a weapon yield and detonation altitude / height of burst, what is the radius of the ring containing a particular PSI. That's the one we corresponded on, but I probably described it poorly at that point. For problem #1, a 750 kiloton explosion @2000 meters will have a 5 PSI ring out to 5821 meters using the Glasstone / interpolate.interp1d method. Nukemap.org yields 6km as the result for 750kt and 2000m HOB. You can visually confirm these results are in the right ballpark using your example interactive_overpressure.py. optimize.fixed_point yields 35.0354xxx as you pointed out, which does not appear to be correct. attaching the code snippet for the problem #1 generic solution: # use this class to calculate overpressure ring radii given a yield (in kt) and altitude of # detonation (in meters) class OverpressureRingCalculator: def __init__(self, max_radius=25000, max_height=5000, radius_interval= 100, height_interval=10): self.max_radius = max_radius self.max_height = max_height self.radius_interval = radius_interval self.height_interval = height_interval # distance/radius intervals self.t = np.arange(0.01, self.max_radius, radius_interval) # height of detonation intervals self.u = np.arange(0, self.max_height, height_interval) def GetOverpressureRadii(self, bomb_yield, height_of_burst, overpressures): ovps = brode_overpressure(bomb_yield, self.t, height_of_burst, 'kT', dunits='m', opunits='psi') f = interpolate.interp1d(ovps, self.t, bounds_error=False, fill_value=0.0) return f(overpressures)

…

On Sun, Dec 18, 2022 at 11:25 AM Edward Geist ***@***.***> wrote: Can you say more about why optimize.fixed_point didn't work for you? I checked just now using Python 2 and this worked: >>> from glasstone.overpressure import brode_overpressure >>> from scipy import optimize >>> optimize.fixed_point(lambda r: brode_overpressure(750, r, 2000, opunits='psi'), 5.0) array(35.0354849668831) >>> optimize.fixed_point(lambda r: brode_overpressure(750, r, 200, opunits='psi'), 5.0) array(530.7211266613074) — Reply to this email directly, view it on GitHub <#7 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ABLDDKELKGKGDRFS35D4FXDWN5CJNANCNFSM5UWBT35A> . You are receiving this because you authored the thread.Message ID: ***@***.***>

[GOFAI]

Wow, now I feel foolish--I had totally misunderstood what optimize.fixed_point does. But scipy.optimize does contain root finding functions that do what you want:

>>> brode_overpressure(750, 5821, 2000, opunits = 'psi')
4.999922789730555
>>> optimize.bisect(lambda r: brode_overpressure(750, r, 2000, opunits='psi') - 5.0, 1000.0, 10000.0)
5820.946857823928

I believe that you can combine this with the minimization techniques scipy provides to solve your second problem. But these solvers are pretty expensive, so it may be faster in practice (depending on what you're trying to do) to just generate tables and interpolate inside them.

[lwillard]

Quick question since I'm pretty much finished with the blast effects part of the work I'm doing and moving onto the fallout portion. Judging by wseg10.py, it looks like you are creating a 11 mile x 11 mile grid (indexed from -1,-1, so ground zero = 0,0), subdivided into 0.1 mile bins. WSEG creates a dose (stored in Z) for each X,Y grid element, and then mapplotlib contour auto-contours...is that basically correct? And lastly, have you done any modeling with Hysplit / pysplit on fallout? BTW, I uploaded my code to: https://github.com/lwillard/pythonnuketools IANAL, but in the readme I attributed it as derived from/using your work. Let me know if the verbiage is satisfactory or if you would want to remove or change anything. Thanks!

…

-Lane

On Sun, Dec 18, 2022 at 1:39 PM Edward Geist ***@***.***> wrote: Wow, now I feel foolish--I had totally misunderstood what optimize.fixed_point does. But scipy.optimize *does* contain root finding functions that do what you want: >>> brode_overpressure(750, 5821, 2000, opunits = 'psi') 4.999922789730555 >>> optimize.bisect(lambda r: brode_overpressure(750, r, 2000, opunits='psi') - 5.0, 1000.0, 10000.0) 5820.946857823928 I believe that you can combine this with the minimization techniques scipy provides to solve your second problem. But these solvers are pretty expensive, so it may be faster in practice (depending on what you're trying to do) to just generate tables and interpolating inside them. — Reply to this email directly, view it on GitHub <#7 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ABLDDKEDTSV3U4LGDK2UATTWN5R7TANCNFSM5UWBT35A> . You are receiving this because you authored the thread.Message ID: ***@***.***>

[GOFAI]

Last year Ivan Stepanov found an old code listing that suggests there's an error in the WSEG-10 model specification I was working from (in the Hanifen AFIT thesis), as well as spotting a couple of other goofs I'd made. He has a version in his fork of the Python library that addresses those issues:
https://github.com/Ivan-Stepanov/glasstone/blob/master/glasstone/fallout.py

And in hindsight the way I implemented WSEG-10 is wildly inefficient. The way I should've done it given that the fallout field has bilateral symmetry would've been to create a spline/interpolation array and then use affine to read from that.

These old "smearing" fallout models and WSEG-10 in particular are of considerable historical interest because they were used in the 1960s and 1970s to evaluate things like McNamara's "assured destruction" criterion, but they aren't the right tool for predicting fallout if you can provide something better. Drawing upon modern transport models like Hysplit or FLEXPART is the right way to go if you have the compute budget. I've looked into it a little bit but I'm not the real expert on it--you should reach out to Stepanov or Sebastien Phillippe, who was doing work on modeling fallout with Hysplit while he was a postdoc at Princeton:
https://sgs.princeton.edu/team/sebastien-philippe

Thanks for acknowledging me in your Python library!

[lwillard]

Small world. I've corresponded a little bit with Ivan on his discord for Nuclear War Simulator (where I am a member of the beta), as well as David Teeter with his work on OPEN-RISOP. Speaking of which, I need to let David know I published my code to a repo, I automated some of the stuff he was doing by hand on his maps. I reached out earlier today to one of the co-authors of the "Predictions of dispersion and deposition of fallout from nuclear testing using the NOAA-HYSPLIT meteorological model" paper to see if there is any insight on how they approached their Hysplit models. I've played around with Hysplit enough to be able to get a puff particle tracking a path against wind speeds and velocity, and pretty confident I can simulate a fair approximation of the fallout pattern, but where I fall short in my technical knowledge is deposition rates and total dose exposure and rates. Compute time isn't a huge issue, this isn't a real time computation need (like NWS), so multi-hour processing jobs are ok if necessary. I loaded the fine grain pop data here: https://sedac.ciesin.columbia.edu/data/set/gpw-v4-admin-unit-center-points-population-estimates-rev11/docs into a postgres db for casualty calculations...and,well, it's not very fast. If push comes to shove and I can't take the results delay, I'll figure out how to do the hard number crunching in the gpu and go buy a 4090 or push the load to AWS.

…

On Sun, Dec 18, 2022 at 5:11 PM Edward Geist ***@***.***> wrote: Last year Ivan Stepanov found an old code listing that suggests there's an error in the WSEG-10 model specification I was working from (in the Hanifen AFIT thesis), as well as spotting a couple of other goofs I'd made. He has a version in his fork of the Python library that addresses those issues: https://github.com/Ivan-Stepanov/glasstone/blob/master/glasstone/fallout.py And in hindsight the way I implemented WSEG-10 is wildly inefficient. The way I should've done it given that the fallout field has bilateral symmetry would've been to create a spline/interpolation array and then use affine to read from that. These old "smearing" fallout models and WSEG-10 in particular are of considerable historical interest because they were used in the 1960s and 1970s to evaluate things like McNamara's "assured destruction" criterion, but they aren't the right tool for predicting fallout if you can provide something better. Drawing upon modern transport models like Hysplit or FLEXPART is the right way to go if you have the compute budget. I've looked into it a little bit but I'm not the real expert on it--you should reach out to Stepanov or Sebastien Phillippe, who was doing work on modeling fallout with Hysplit while he was a postdoc at Princeton: https://sgs.princeton.edu/team/sebastien-philippe Thanks for acknowledging me in your Python library! — Reply to this email directly, view it on GitHub <#7 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ABLDDKADAHWNWC2WIVLA2MDWN6KZPANCNFSM5UWBT35A> . You are receiving this because you authored the thread.Message ID: ***@***.***>