1

I'm trying to simulate a model very used with analyses of magnetic materials. Initially, I wrote this code in python but did not obtain results; then, I decided to use Mathematica to make the same simulation [1] and It seems to work properly, at least It shows a curve similar to expected, but when It's translated to python (using the lmfit package and scipy.integrate.quad) it only shows a line. What else should I consider to make it work properly in python?

Mathematica code

Here what I have:

from lmfit import minimize, Minimizer, Parameters, Parameter, report_fit
from numpy import loadtxt, vectorize, sqrt, log, inf, exp, pi, tanh,     linspace
from scipy.integrate import quad
from scipy.constants import codata
import matplotlib.pyplot as plt
from mpmath import coth


T1 = linspace(1, 300, 100)

def PDF(Dmag, v):
#    v = params.valuesdict()
    pdf=(1/((Dmag*1e-9)*v['w']*sqrt(2*pi)))*exp(-((log((Dmag*1e-9)/Dmag))**2)/(2*v['w']**2))*1e-9
    return pdf

def mZFC(Dmag, args):
    T1, params = args
    v  = params.valuesdict()
    V  = ((pi/6)*(Dmag*1e-9)**3)
    nu = v['v0']*exp(-(V*v['Kv'])/(v['kB']*T1))
    dt = v['kB']*T1**2/((V*v['Kv'])*v['vT']);
    return (v['B']*(v['Ms']*V)**2/(3*(V*v['Kv'])))*(exp(-nu*dt) +  ((V*v['Kv'])/(v['kB']*T1))*(1 - exp(-nu*dt)))*PDF(Dmag, v)

def curveZFC(T1, params):
    return params['NT']*quad(mZFC, 0.0, inf, [T1, params])[0]

vcurveZFC = vectorize(curveZFC, excluded=set([1]))


kB = codata.value('Boltzmann constant')

params = Parameters()
params.add('Dm' , value= 3.2e-9 , vary= True)
params.add('w'  , value= 0.26    , vary= True)
params.add('Ms' , value= 1350000   , vary= True)
params.add('T'  , value= 300     , vary= False)
params.add('kB' , value= kB      , vary= False)
params.add('NT' , value= 1.7e12  , vary= True)
#m_ZFC params
params.add('B'   , value= 0.0050, vary= False)
params.add('v0'  , value= 1e10 , vary= False)
params.add('Kv'  , value= 178000  , vary= False)
params.add('vT'  , value= 0.05 , vary= False)


y2init = vcurveZFC(T1, params)

plt.plot(T1, y2init, 'r*')
plt.show()

Any help would be appreciated. Thanks!

Ernesto
  • 21
  • 4
  • Your dropbox link is private, you have to share it if we are going to be able to view the code. – Daniel Jun 25 '17 at 19:24
  • Which version of python are you using? – Frank Wilson Jun 25 '17 at 19:27
  • I changed the link to download and my python version is 2.7.12 – Ernesto Jun 25 '17 at 19:44
  • isn't this the same question (or very similar code) as https://stackoverflow.com/questions/44464981 ? I think the problem earlier was with your `vcurveZFC` function not working as you expect. That is, you already *know* where the problem lies, so can ask a more specific question. First, does `quad` do what you expect? Does `curveZFC` do what you expect, and so forth. – M Newville Jun 27 '17 at 12:36
  • Dear @MNewville, The other code only worked when rescaling one of the variables by 1e-9. It was not necessary to make any additional changes. It seems that the power of kB introduces some difficulties in the calculation. Now, I'm working in the second part of the overall code. I only included the mZFC function, which works properly in Mathematica, but in python, It does not! I tried to rewrite the equations and other solutions but none has worked so far. – Ernesto Jun 27 '17 at 14:36

0 Answers0