Skip to content

Instantly share code, notes, and snippets.

@pkaf
Created August 6, 2015 12:24
Show Gist options
  • Save pkaf/843b1e7dcfbbc7c2f050 to your computer and use it in GitHub Desktop.
Save pkaf/843b1e7dcfbbc7c2f050 to your computer and use it in GitHub Desktop.
Vectorization of scipy.integrate.quad function
Display the source blob
Display the rendered blob
Raw
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Zach Cano question in Python Astronomer Facebook page"
]
},
{
"cell_type": "code",
"execution_count": 61,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"array([ 1. , 1. , 1. , 1. , 1. ,\n",
" 1. , 1. , 1. , 1. , 1. ,\n",
" 1. , 1. , 1. , 1. , 1. ,\n",
" 1. , 1. , 1. , 1. , 1. ,\n",
" 1. , 1. , 1. , 1. , 1. ,\n",
" 1. , 1. , 1. , 0.99999999, 0.99999996,\n",
" 0.99999982, 0.99999914, 0.99999627, 0.99998505, 0.99994483,\n",
" 0.99981232, 0.99941155, 0.99829961, 0.9954717 , 0.9888859 ,\n",
" 0.97486006, 0.94759109, 0.89930752, 0.82170465, 0.70903954,\n",
" 0.56239836, 0.3934384 , 0.2251424 , 0.08773738, 0.01015117,\n",
" 0.01015117, 0.08773738, 0.2251424 , 0.3934384 , 0.56239836,\n",
" 0.70903954, 0.82170465, 0.89930752, 0.94759109, 0.97486006,\n",
" 0.9888859 , 0.9954717 , 0.99829961, 0.99941155, 0.99981232,\n",
" 0.99994483, 0.99998505, 0.99999627, 0.99999914, 0.99999982,\n",
" 0.99999996, 0.99999999, 1. , 1. , 1. ,\n",
" 1. , 1. , 1. , 1. , 1. ,\n",
" 1. , 1. , 1. , 1. , 1. ,\n",
" 1. , 1. , 1. , 1. , 1. ,\n",
" 1. , 1. , 1. , 1. , 1. ,\n",
" 1. , 1. , 1. , 1. , 1. ])"
]
},
"execution_count": 61,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"import numpy as np\n",
"from scipy.integrate import quad\n",
"\n",
"def A(z,tm,tni):\n",
" y=tm/(2*tni)\n",
" return 2*z*np.exp((-2*z*y)+(z**2))\n",
"\n",
"def B(z,tm,tni,tco):\n",
" y=tm/(2*tni)\n",
" s=(tm*(tco-tni))/(2*tco*tni)\n",
" return 2*z*np.exp(-2*z*y+2*z*s+z**2)\n",
"\n",
"def Arnett(t,tm,tni,tco,Mni,Eni,Eco): \n",
" x=t/tm\n",
" f,jk=quad(A,0, x, args=(tm,tni,)) #integral of A\n",
" h,jk=quad(B,0, x, args=(tm,tni,tco,)) #integral of B\n",
" g=np.exp((-(x/tm)**2))\n",
" return Mni*g*((Eni-Eco)*f+Eco*h)\n",
"\n",
"#MAGIC BIT\n",
"Arnt = np.vectorize(Arnett)\n",
"\n",
"#Fiducial values to test the code\n",
"t = np.linspace(-10,10,100)\n",
"z = np.linspace(-1,1,100)\n",
"quad( A, 0, 2, args=(1,1,))\n",
"\n",
"#Your final result still function of t\n",
"Arnt(t, 1, 1, 1, 1, 1, 1 )"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 2",
"language": "python",
"name": "python2"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 2
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython2",
"version": "2.7.9"
}
},
"nbformat": 4,
"nbformat_minor": 0
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment