update readme and add notebook

amanmdesai · Feb 7, 2023 · 7abc4a5 · 7abc4a5
1 parent 282048f
commit 7abc4a5
Show file tree

Hide file tree

Showing 2 changed files with 225 additions and 3 deletions.
diff --git a/README.md b/README.md
@@ -14,10 +14,94 @@ Aman Desai
 A package to read SPARC data for Galactic Rotation Curves
 
 ## Installation
+
 ```bash
-pip install git+https:https://github.com/amanmdesai/pygrc
+pip install pygrc
+```
+
+## Example Usage
+
+- 1. Import libraries
+
+```python
+import math
+import numpy as np
+import matplotlib.pyplot as plt
+import pygrc as gr
+```
+
+- 2. Read Sparc Data
+
+```python
+df = gr.Reader.read(filepath="/home/amdesai/Physics/data/NGC5055_rotmod.dat")
+gr.Plot().overlap(df,"Rad",["Vobs","Vgas","Vbul","Vdisk"],"observed velocity")
+df
+```
+
+- 3.  Some sample function
+```python
+# mond and mass function are based on Galaxies 2018, 6(3), 70; https://doi.org/10.3390/galaxies6030070
+
+def mass(r, M0, R0):
+    return M0*(1- (1+(r/R0))*np.exp(-r/R0))
+
+def mond(r, M0, rc, R0,b, beta):
+    a = 1.2e-10
+    G = 4.300e-6 #parsecs
+    m = mass(r,M0, R0)
+    f = (G*m/r)*(1 + b*(1+(r/rc)))#*10e-3
+    return np.sqrt(f)*10e-3
+
+def newton(r, M0, rc, R0,beta):
+    a = 1.2e-10
+    G = 4.30e-6 #parsecs
+    m = mass(r,M0, R0)
+    f = G*m/r
+    #f = (G*m/r)*(1/np.sqrt(2)) * np.sqrt(1 + np.sqrt(1 + (r**4) * (2*a/(G*m))**2))
+    return np.sqrt(f)*10e-3
+```
+
+- 4a. Perform fit  (uses iminuit in the background)
+
+```python
+r = np.linspace(1e-5,df['Rad'].max(),2000)
+
+m_1=gr.Fit(df['Rad'],df['Vobs'],1.,1.,3,.35,1.)
+
+m_2=gr.Fit(df['Rad'],df['Vobs'],1.,1.,3,1.)
+```
+
+- 4b. Continue:
+
+```python
+m1= m_1.fit_lsq(mond, [(1.,1e17),(1.,10.),(2.,5.),(0.1,2),(0.1,2)],df['errV'],[False,False,True,True,True])
+m1
 ```
 
-SPARC Data can be obtained from here: 
-http://astroweb.cwru.edu/SPARC/
+- 4c. Continue:
+
+```python
+m2 = m_2.fit_lsq(newton, [(1.,1e16),(1.,10.),(1,10),(0.1,2)],df['errV'],[False,True,True,True])
+m2
+```
+
+- 5a. draw plots
+```python
+m_2.get_profile(m2,'M0')
+m_1.draw_contours(m1,['M0','rc'])
+```
+
+- 6b.
+
+```python
+fig, ax =plt.subplots()
+gr.Plot().plot_grc(df,m1,mond,'MOND','NGC',ax)
+gr.Plot().plot_grc(df,m2,newton,'Newton','NGC',ax)
+plt.savefig('1.pdf')
+```
+
+
+References:
 
+- SPARC Data can be obtained from here: http://astroweb.cwru.edu/SPARC/
+- Mond and Mass function are based on Galaxies 2018, 6(3), 70; https://doi.org/10.3390/galaxies6030070
diff --git a/notebooks/example.ipynb b/notebooks/example.ipynb
@@ -0,0 +1,138 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "e7a1fc45",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import math\n",
+    "import numpy as np\n",
+    "import matplotlib.pyplot as plt\n",
+    "import pygrc as gr"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "7e27a301",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "df = gr.Reader.read(filepath=\"/home/amdesai/Physics/data/NGC5055_rotmod.dat\")\n",
+    "gr.Plot().overlap(df,\"Rad\",[\"Vobs\",\"Vgas\",\"Vbul\",\"Vdisk\"],\"observed velocity\")\n",
+    "df"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "642a96a4",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# some example functions based on Galaxies 2018, 6(3), 70; https://doi.org/10.3390/galaxies6030070\n",
+    "\n",
+    "def mass(r, M0, R0):\n",
+    "    return M0*(1- (1+(r/R0))*np.exp(-r/R0))\n",
+    "\n",
+    "def mond(r, M0, rc, R0,b, beta):\n",
+    "    a = 1.2e-10\n",
+    "    G = 4.300e-6 #parsecs \n",
+    "    m = mass(r,M0, R0)\n",
+    "    f = (G*m/r)*(1 + b*(1+(r/rc)))#*10e-3\n",
+    "    return np.sqrt(f)*10e-3\n",
+    "\n",
+    "def newton(r, M0, rc, R0,beta):\n",
+    "    a = 1.2e-10\n",
+    "    G = 4.30e-6 #parsecs \n",
+    "    m = mass(r,M0, R0)\n",
+    "    f = G*m/r\n",
+    "    #f = (G*m/r)*(1/np.sqrt(2)) * np.sqrt(1 + np.sqrt(1 + (r**4) * (2*a/(G*m))**2))\n",
+    "    return np.sqrt(f)*10e-3\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "6df730b3",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "r = np.linspace(1e-5,df['Rad'].max(),2000)\n",
+    "\n",
+    "m_1=gr.Fit(df['Rad'],df['Vobs'],1.,1.,3,.35,1.)\n",
+    "\n",
+    "m_2=gr.Fit(df['Rad'],df['Vobs'],1.,1.,3,1.)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "0deaae95",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "m1= m_1.fit_lsq(mond, [(1.,1e17),(1.,10.),(2.,5.),(0.1,2),(0.1,2)],df['errV'],[False,False,True,True,True])\n",
+    "m1"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "64e5e876",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "m2 = m_2.fit_lsq(newton, [(1.,1e16),(1.,10.),(1,10),(0.1,2)],df['errV'],[False,True,True,True])\n",
+    "m2"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "efb56302",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "m_2.get_profile(m2,'M0')\n",
+    "m_1.draw_contours(m1,['M0','rc'])"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "98ed80d2",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fig, ax =plt.subplots()\n",
+    "gr.Plot().plot_grc(df,m1,mond,'MOND','NGC',ax)\n",
+    "gr.Plot().plot_grc(df,m2,newton,'Newton','NGC',ax)\n",
+    "plt.savefig('1.pdf')"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3 (ipykernel)",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.10.6"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}