#1073 closed defect (invalid)
Discrepancy in code for polymer_excl_volume
Reported by: | smk78 | Owned by: | |
---|---|---|---|
Priority: | major | Milestone: | sasmodels 1.0 |
Component: | sasmodels | Keywords: | |
Cc: | Work Package: | SasView Bug Fixing |
Description ¶
Reported by pkienzle:
Just went through the code for polymer_excl_volume. The result should match older versions of SasView, but this does not match the documentation, nor does it match Igor. For some reason, when Jae He added the model he scaled the gammainc term by gamma. Maybe Paul or Mat has some memory of why?
Here the code that Jae He added:
$ git show -p ec1c5540 2010-10-04 ec1c5540 added polymerexclVolume model and its test Jae Cho diff --git a/sansmodels/src/sans/models/PolymerExclVolume.py b/sansmodels/src/sans/models/PolymerExclVolume.py ... + nu = 1.0 / mm + + Xx = x * x * rg * rg *(2.0 * nu + 1.0) * (2.0 * nu + 2.0) / 6.0 + onu = 1.0 / nu + o2nu = 1.0 /(2.0 * nu) + Ps =(1.0 / (nu * pow(Xx,o2nu))) * (gamma(o2nu)*gammainc(o2nu,Xx) - \ + 1.0 / pow(Xx,o2nu) * gamma(onu)*gammainc(onu,Xx)) + + if x == 0: + Ps = 1.0 + + return (sc * Ps + bg);
Here's what is in Igor 7.12, Analysis/Models/NewModels_2010/PolymerExcludVol_v40.ipf:
nu=1/m Xx=qval^2*Rg^2*(2*nu+1)*(2*nu+2)/6 onu=1.0/(nu) o2nu=1.0/(2*nu) Ps=(1/(nu*Xx^o2nu))*(gammaInc(o2nu,Xx,0)-(1/Xx^o2nu)*gammaInc(onu,Xx,0)) Debye=(2/Xx^2)*(exp(-Xx)-1+Xx) if(qval == 0) Ps = 1 endif inten = I0*Ps + bgd Igor gammaInc(a, x, 0) = gamma(a) - gammaInc(a, x) = int_0^x e^-t t^(a-1) dt = scipy.special.gammainc(a, x).
Note: The current version looks different from Jae He's because I updated it during one of my cleanup passes. I'm pretty sure that I didn't change the result, but probably should revert it so that the equations in the code better match the equations in the documentation.
Change History (11)
Changed 7 years ago by smk78
comment:1 Changed 7 years ago by smk78
SasView documentation for polymer_excl_volume model references:
H Benoit, Comptes Rendus, 245 (1957) 2244-2247
B Hammouda, SANS from Homogeneous Polymer Mixtures - A Unified Overview, Advances in Polym. Sci. 106(1993) 87-133
Second of these is attached to this ticket.
comment:2 Changed 7 years ago by krzywon
From Mike Hore:
- The reason that gammainc is scaled by gamma in the SasView function is that gammainc in scipy is normalized by gamma, so that is 100% correct to do.
- Looking at the function in the version of SasView (4.1.2) that I have installed currently, the issue is in the following block of code:
usub = (q*rg)**2 * (2.0/porod_exp + 1.0) * (2.0/porod_exp + 2.0)/6.0 with errstate(divide='ignore', invalid='ignore'): upow = power(usub, -0.5*porod_exp) result= (porod_exp*upow * (gamma(0.5*porod_exp)*gammainc(0.5*porod_exp, usub) - upow*gamma(porod_exp)*gammainc(porod_exp, usub)))
- I think it should be something like the following instead:
usub = (q*rg)**2 * (2.0/porod_exp + 1.0) * (2.0/porod_exp + 2.0)/6.0 with errstate(divide='ignore', invalid='ignore'): upow1 = power(usub, -0.5*porod_exp) upow2 = power(usub, -porod_exp) result= (porod_exp*upow1 * (gamma(0.5*porod_exp)*gammainc(0.5*porod_exp, usub)upow2*gamma(porod_exp)*gammainc(porod_exp, usub)))
- Because U is raised to a different power for the first and second terms in the correct version of the function (½*nu for the first, and 1/nu for the second).
A second email after he sent the above:
- I see where the disagreement is from. Boualem’s paper from 1993 has the form that’s currently in SasView, but more recent publications like https://www.ncnr.nist.gov/staff/hammouda/publications/2017_hammouda_kim_j_mol_liq.pdf and mine https://www.ncnr.nist.gov/staff/hammouda/publications/2013_hore_et_al_macromolecules.pdf have the form that I’m talking about. So I guess it’s a matter of which is correct!
comment:3 Changed 7 years ago by smk78
Interestingly, https://www.ncnr.nist.gov/staff/hammouda/publications/2013_hore_et_al_macromolecules.pdf references the polymer with excluded volume form factor with exactly the same two references that are in the SasView documentation!!!
comment:4 Changed 7 years ago by krzywon
I just spoke with Boualem and we went through the math together. In the 1993 paper, [1/vX^(1/2v)] has been factored out. When multiplied through, the resulting equation is identical to the equation in the 2013 and 2017 papers. The calculation currently within SasView matches that in the 1993 paper so I think the model might be correct.
comment:5 Changed 7 years ago by smk78
So can we close this ticket now, or is there still an issue?
comment:6 Changed 7 years ago by krzywon
- Resolution set to invalid
- Status changed from new to closed
I contacted Mike Hore and he agrees, the mathematics is correct. Yes, we can close this ticket.
comment:7 Changed 7 years ago by smk78
I have pushed some updates to the documentation to Master to reflect what has been learnt from this exchange.
comment:8 Changed 7 years ago by Paul Kienzle <pkienzle@…>
comment:9 Changed 7 years ago by Paul Kienzle <pkienzle@…>
comment:10 Changed 7 years ago by butler
- Component changed from SasView to sasmodels
- Milestone changed from SasView 4.3.0 to sasmodels 1.0
Hammouda reference for polymer_excl_volume model