```python
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
```

## Curve fitting

\begin{equation}
RMSE=\sqrt{\frac{1}{n}\sum{(y-\hat{y})^2}}\\
R^2=1-\frac{\sum{(y-\hat{y})^2}}{\sum{(y-\bar{y})^2}}
\end{equation}

## Example:

### 1.Fit gas equation
\begin{equation}
PV^\gamma=constant(K)\\
logP=logK-\gamma logV\\
y=a+bx
\end{equation}


```python
df = pd.read_excel('data/gas.xlsx')
df
```




<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }

    .dataframe tbody tr th {
        vertical-align: top;
    }

    .dataframe thead th {
        text-align: right;
    }
</style>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>volume(cm3)</th>
      <th>pressure(gm/cm2)</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>0</th>
      <td>54.3</td>
      <td>61.2</td>
    </tr>
    <tr>
      <th>1</th>
      <td>61.8</td>
      <td>49.2</td>
    </tr>
    <tr>
      <th>2</th>
      <td>72.4</td>
      <td>37.6</td>
    </tr>
    <tr>
      <th>3</th>
      <td>88.7</td>
      <td>28.4</td>
    </tr>
    <tr>
      <th>4</th>
      <td>118.6</td>
      <td>19.2</td>
    </tr>
    <tr>
      <th>5</th>
      <td>194.0</td>
      <td>10.1</td>
    </tr>
  </tbody>
</table>
</div>




```python
plt.figure(figsize = [6,6])
plt.scatter(df['volume(cm3)'],df['pressure(gm/cm2)'])
plt.xlabel('volume(cm3)')
plt.ylabel('pressure(gm/cm2)')
plt.title('Gas')
plt.show()
```


```python
n=6
x=np.log(df['volume(cm3)'])
y=np.log(df['pressure(gm/cm2)'])
z=np.polyfit(x, y, 1) # linear  fit =y=a+bx
p=np.poly1d(z)  #  [a,b]
ey=p(x) # expected value
a="{:.3}".format(p[0])
b="{:.3}".format(p[1])
K=np.exp(p[0])
gamma=p[1]

RMSE=np.sqrt(np.sum((y-ey)**2)/n)
R_square=1-(np.sum(np.sum(np.sum((y-ey)**2))/np.sum((y-np.mean(y))**2)))
rmse="{:.2}".format(RMSE)
r2="{:.2}".format(R_square)


plt.figure(figsize = [6,6])
plt.plot(x, y, 'ro',label='observed')
plt.plot(x, ey, 'b-',label='fitted, a= '+str(a)+', b='+str(b))
plt.xlabel('log(V)')
plt.ylabel('log(P)')
plt.figtext(0.2, 0.3, 'RMSE = '+str(rmse)+',r2='+str(r2))
plt.legend()
plt.title('Fitted curve')
plt.show()
print('gamma =',"{:.2}".format(-gamma))
print('K=',"{:.2}".format(K))

```

### 2.Fit 
\begin{equation}
f(\theta)=A+B\cos\theta+C\cos^2\theta\\
y=A+Bx+Cx^2
\end{equation}


```python
df1 = pd.read_excel('data/scattering.xlsx')
df1
```




<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }

    .dataframe tbody tr th {
        vertical-align: top;
    }

    .dataframe thead th {
        text-align: right;
    }
</style>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>theta(degree)</th>
      <th>f(theta)</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>0</th>
      <td>30</td>
      <td>11</td>
    </tr>
    <tr>
      <th>1</th>
      <td>45</td>
      <td>13</td>
    </tr>
    <tr>
      <th>2</th>
      <td>90</td>
      <td>16</td>
    </tr>
    <tr>
      <th>3</th>
      <td>120</td>
      <td>17</td>
    </tr>
    <tr>
      <th>4</th>
      <td>150</td>
      <td>14</td>
    </tr>
  </tbody>
</table>
</div>




```python
plt.figure(figsize = [6,6])
plt.scatter(df1['theta(degree)'],df1['f(theta)'])
plt.xlabel(r'$\theta$(degree)')
plt.ylabel(r'f($\theta$)')
plt.title('Scattering')
plt.show()
```


```python
n1=5
x1=np.cos(df1['theta(degree)']*np.pi/180)
y1=df1['f(theta)']
z1=np.polyfit(x1, y1, 2) # second order  fit y=a+bx+cx2
p1=np.poly1d(z1)  #  [c,b,a]
ey1=p1(x1) #expected value
A="{:.3}".format(p1[0])
B="{:.3}".format(p1[1])
C="{:.3}".format(p1[2])


RMSE1=np.sqrt(np.sum((y1-ey1)**2)/n1)
R_square1=1-(np.sum(np.sum(np.sum((y1-ey1)**2))/np.sum((y1-np.mean(y1))**2)))
rmse1="{:.2}".format(RMSE1)
r21="{:.2}".format(R_square1)


plt.figure(figsize = [6,6])
plt.plot(x1, y1, 'ro',label='observed')
plt.plot(x1, ey1, 'b-',label='fitted, A= '+str(A)+', B='+str(B)+',C='+str(C))
plt.xlabel(r'cos($\theta$)')
plt.ylabel(r'f($\theta$)')
plt.figtext(0.2, 0.3, 'RMSE = '+str(rmse1)+',r2='+str(r21))
plt.legend()
plt.title('Fitted curve')
plt.show()

```

### 3.Determination of decay conctant and half life
\begin{equation}
N=N_0~e^{-\lambda~t}\\
logN=logN_0-\lambda~t\\
y=A+BX
\end{equation}


```python
df2 = pd.read_excel('data/env-dec.xlsx')
df2
```




<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }

    .dataframe tbody tr th {
        vertical-align: top;
    }

    .dataframe thead th {
        text-align: right;
    }
</style>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>t [s]</th>
      <th>N [s-1]</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>0</th>
      <td>0</td>
      <td>32</td>
    </tr>
    <tr>
      <th>1</th>
      <td>5</td>
      <td>28</td>
    </tr>
    <tr>
      <th>2</th>
      <td>10</td>
      <td>29</td>
    </tr>
    <tr>
      <th>3</th>
      <td>15</td>
      <td>28</td>
    </tr>
    <tr>
      <th>4</th>
      <td>20</td>
      <td>25</td>
    </tr>
    <tr>
      <th>...</th>
      <td>...</td>
      <td>...</td>
    </tr>
    <tr>
      <th>78</th>
      <td>390</td>
      <td>1</td>
    </tr>
    <tr>
      <th>79</th>
      <td>395</td>
      <td>1</td>
    </tr>
    <tr>
      <th>80</th>
      <td>400</td>
      <td>2</td>
    </tr>
    <tr>
      <th>81</th>
      <td>405</td>
      <td>2</td>
    </tr>
    <tr>
      <th>82</th>
      <td>410</td>
      <td>1</td>
    </tr>
  </tbody>
</table>
<p>83 rows × 2 columns</p>
</div>




```python
plt.figure(figsize = [6,6])
plt.scatter(df2['t [s]'],df2['N [s-1]'])
plt.xlabel('time(s)')
plt.ylabel('Count($s^{-1}$)')
plt.title('Decay curve')
plt.show()
```


```python
n2=83
x2=df2['t [s]']
y2=np.log(df2['N [s-1]'])
z2=np.polyfit(x2, y2, 1) # second order  fit y=a+bx
p2=np.poly1d(z2)  #  [c,b,a]
ey2=p2(x2) #expected value
A1="{:.3}".format(p2[0])
B1="{:.3}".format(-p2[1])


RMSE2=np.sqrt(np.sum((y2-ey2)**2)/n2)
R_square2=1-(np.sum(np.sum(np.sum((y2-ey2)**2))/np.sum((y2-np.mean(y2))**2)))
rmse2="{:.2}".format(RMSE2)
r22="{:.2}".format(R_square2)


plt.figure(figsize = [6,6])
plt.plot(x2, y2, 'ro',label='observed')
plt.plot(x2, ey2, 'b-',label='fitted, A= '+str(A1)+', B='+str(B1))
plt.xlabel('time(sec)')
plt.ylabel('ln(A)')
plt.figtext(0.2, 0.3, 'RMSE = '+str(rmse2)+',r2='+str(r22))
plt.legend()
plt.title('Fitted curve')
plt.show()
```


```python
print('decay constant =%0.6f /sec'%-p2[1])
```

    decay constant =0.008246 /sec



```python
T=-0.693/p2[1]
print('half life = %0.2f sec'%T)
```

    half life = 84.04 sec


### 4.Absorption of coefficient
\begin{equation}
N=N_0~e^{-\mu~x}\\
logN=logN_0-\mu~x\\
y=A+BX
\end{equation}


```python
df3 = pd.read_excel('data/env-beta.xlsx')
df3
```




<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }

    .dataframe tbody tr th {
        vertical-align: top;
    }

    .dataframe thead th {
        text-align: right;
    }
</style>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>n</th>
      <th>N [s-1]</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>0</th>
      <td>0</td>
      <td>8.0</td>
    </tr>
    <tr>
      <th>1</th>
      <td>1</td>
      <td>6.0</td>
    </tr>
    <tr>
      <th>2</th>
      <td>2</td>
      <td>5.0</td>
    </tr>
    <tr>
      <th>3</th>
      <td>3</td>
      <td>4.0</td>
    </tr>
    <tr>
      <th>4</th>
      <td>4</td>
      <td>3.0</td>
    </tr>
    <tr>
      <th>5</th>
      <td>5</td>
      <td>2.0</td>
    </tr>
    <tr>
      <th>6</th>
      <td>6</td>
      <td>2.0</td>
    </tr>
    <tr>
      <th>7</th>
      <td>7</td>
      <td>2.0</td>
    </tr>
    <tr>
      <th>8</th>
      <td>8</td>
      <td>2.0</td>
    </tr>
    <tr>
      <th>9</th>
      <td>9</td>
      <td>1.0</td>
    </tr>
    <tr>
      <th>10</th>
      <td>10</td>
      <td>1.0</td>
    </tr>
    <tr>
      <th>11</th>
      <td>11</td>
      <td>1.0</td>
    </tr>
    <tr>
      <th>12</th>
      <td>12</td>
      <td>0.5</td>
    </tr>
    <tr>
      <th>13</th>
      <td>13</td>
      <td>0.5</td>
    </tr>
    <tr>
      <th>14</th>
      <td>14</td>
      <td>0.5</td>
    </tr>
    <tr>
      <th>15</th>
      <td>15</td>
      <td>0.5</td>
    </tr>
    <tr>
      <th>16</th>
      <td>16</td>
      <td>1.0</td>
    </tr>
    <tr>
      <th>17</th>
      <td>17</td>
      <td>0.5</td>
    </tr>
    <tr>
      <th>18</th>
      <td>18</td>
      <td>0.4</td>
    </tr>
    <tr>
      <th>19</th>
      <td>19</td>
      <td>0.3</td>
    </tr>
    <tr>
      <th>20</th>
      <td>20</td>
      <td>0.1</td>
    </tr>
  </tbody>
</table>
</div>




```python
plt.figure(figsize = [6,6])
plt.scatter(0.03*df3['n'],df3['N [s-1]'])# one sheet of Al is 0.03 mm thick
plt.xlabel('Thickness(mm)')
plt.ylabel('Count($s^{-1}$)')
plt.title('Absorption curve')
plt.show()
```


```python
n3=20
x3=0.03*df3['n']
y3=np.log(df3['N [s-1]'])
z3=np.polyfit(x3, y3, 1) # second order  fit y=a+bx
p3=np.poly1d(z3)  #  [c,b,a]
ey3=p3(x3) #expected value
A2="{:.3}".format(p3[0])
B2="{:.3}".format(-p3[1])


RMSE3=np.sqrt(np.sum((y3-ey3)**2)/n3)
R_square3=1-(np.sum(np.sum(np.sum((y3-ey3)**2))/np.sum((y3-np.mean(y2))**2)))
rmse3="{:.2}".format(RMSE3)
r23="{:.2}".format(R_square3)


plt.figure(figsize = [6,6])
plt.plot(x3, y3, 'ro',label='observed')
plt.plot(x3, ey3, 'b-',label='fitted, A= '+str(A2)+', B='+str(B2))
plt.xlabel('Thickness(mm)')
plt.ylabel('ln(N)')
plt.figtext(0.2, 0.3, 'RMSE = '+str(rmse3)+',r2='+str(r23))
plt.legend()
plt.title('Fitted curve')
plt.show()
```


```python
print('absorption coefficient =%0.3f /mm'%-p3[1])
```

    absorption coefficient =5.695 /mm



```python
density=2.7 # in cgs
mu_m=-10*p3[1]/density # in cgs
print('mass absorption coefficient =%0.3f sq.cm/gm'%mu_m)
```

    mass absorption coefficient =21.091 sq.cm/gm



```python

```

104efcbb3a2897b7cd178c9873ddd8077183f565
ipynb
Jupyter Notebook
PMS_curve.ipynb
joshidot/NPS
0b5b7dde9b5a9769c8a437d193b210545f9344ca
MIT
PMS_curve.ipynb
joshidot/NPS
0b5b7dde9b5a9769c8a437d193b210545f9344ca
MIT
PMS_curve.ipynb
joshidot/NPS
0b5b7dde9b5a9769c8a437d193b210545f9344ca
MIT
Qwen/Qwen-72B
1. YES
2. YES