{"nbformat":4,"nbformat_minor":0,"metadata":{"colab":{"provenance":[]},"kernelspec":{"name":"python3","display_name":"Python 3"},"language_info":{"name":"python"}},"cells":[{"cell_type":"markdown","source":["Developed by Megan Renz, Rebeckah Fussell, and Natasha Holmes for Cornell Physics Labs."],"metadata":{"id":"eoi0PcNapPPW"}},{"cell_type":"markdown","source":["# Packages"],"metadata":{"id":"ZchL-AIepbqi"}},{"cell_type":"code","source":["import sys\n","\n","%matplotlib inline\n","#%matplotlib notebook\n","import matplotlib.pyplot as plt\n","import numpy as np\n","import csv\n","from scipy.optimize import minimize, rosen, rosen_der, curve_fit\n","from matplotlib.widgets import TextBox\n","import pandas as pd"],"metadata":{"id":"7Dv5gwG9pcth"},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":["# Fitting"],"metadata":{"id":"iDrCSP8mpQlf"}},{"cell_type":"code","source":["\n","# def chiSquared2(args,x, y, dy, f ):\n","#     '''Function Chi-Squared.\n","#     x, y and dy are numpy arrays, referring to x, y and the uncertainty in y respectively.\n","#     f is the function we are fitting.\n","#     args are the arguments of the function we have fit.\n","#     '''\n","#     return 1/(len(x)*np.sum((f(x, args)-y)**2/dy**2))\n","\n","def chiSquared(x, y, dy, f, args):\n","    '''Function Chi-Squared.\n","    x, y and dy are numpy arrays, referring to x, y and the uncertainty in y respectively.\n","    f is the function we are fitting.\n","    args are the arguments of the function we have fit.\n","    '''\n","    return 1/(len(x))*np.sum((f(x, args)-y)**2/dy**2)\n","\n","def linear(x,args):\n","    '''\n","    A special case of Poly.\n","    '''\n","    return x*args[0]+args[1]\n","\n","def linearFit(x, m,b):\n","    return m*x+b\n","\n","def poly(x, args):\n","    '''\n","    returns the value of the polynomial sum (x**i*args[i])\n","    '''\n","    total=x**0*args[0]\n","    for i in range(1,len(args)):\n","        total+=x**i*args[i]\n","    return total\n","\n","def polyFit(x, A, B, C):\n","  return A*(x**2)+B*x+C\n","\n","def autoFit(x=[], y=[], dy=[], title=\"Use title= in your call\", xaxis=\"Use xaxis= in your call\", yaxis=\"Use yaxis= in your call\"):\n","    x0=[0,0]\n","    x=np.array(x)\n","    y=np.array(y)\n","    dy=np.array(dy)\n","\n","    args,cov=curve_fit(linearFit, x, y, x0, dy)\n","\n","    m=args[0]\n","    b=args[1]\n","    dm = np.sqrt(cov[0][0])\n","    db = np.sqrt(cov[1][1])\n","\n","    fig, ax=plt.subplots(1,2, figsize=(12,6), dpi=75)\n","    line,=ax[0].plot(x,linear(x, [m,b]))\n","    data=ax[0].errorbar(x,y, dy, fmt='.k', capsize=5, markersize=10)\n","    ax[0].set_title(title)\n","    ax[0].set_xlabel(xaxis)\n","    ax[0].set_ylabel(yaxis)\n","    ax[1].set_title(\"Residuals\")\n","    ax[1].set_xlabel(xaxis)\n","    ax[1].set_ylabel(\"point -- line\")\n","    residuals=y-linear(x, [m,b])\n","    res=ax[1].errorbar(x,residuals, dy, fmt='.k', capsize=5, markersize=10)\n","    ax[1].grid(True, which='both')\n","    plt.show()\n","\n","    print(\"chi-squared value =  {}\".format(chiSquared(x,y, dy, linear, [m, b]).round(3)))\n","    print(\"Best fit line :  \")\n","    print(\"y=\", round(args[0],4),\"x\", \"+\", round(args[1],4))\n","    print(\"m={}+\\-{}\".format(m.round(3),dm.round(3),))\n","    print(\"b={}+\\-{}\".format(b.round(3),db.round(3),))\n","\n","def manualFit(x=[], y=[], dy=[], m=[], b=[], title=\"Use title= in your call\", xaxis=\"Use xaxis= in your call\", yaxis=\"Use yaxis= in your call\"):\n","    x0=[0,0]\n","    x=np.array(x)\n","    y=np.array(y)\n","    dy=np.array(dy)\n","\n","    fig, ax=plt.subplots(1,2, figsize=(12,6), dpi=75)\n","    line,=ax[0].plot(x,linear(x, [m,b]))\n","    data=ax[0].errorbar(x,y, dy, fmt='.k')\n","    ax[0].set_title(title)\n","    ax[0].set_xlabel(xaxis)\n","    ax[0].set_ylabel(yaxis)\n","    ax[1].set_title(\"Residuals\")\n","    ax[1].set_xlabel(xaxis)\n","    ax[1].set_ylabel(\"point -- line\")\n","    residuals=y-linear(x, [m,b])\n","    res=ax[1].errorbar(x,residuals, dy, fmt='.k')\n","    ax[1].grid(True, which='both')\n","    plt.show()\n","\n","    print(\"chi-squared value =  {}\".format(chiSquared(x,y, dy, linear, [m, b]).round(3)))\n","    print(\"Best fit line :  \")\n","    print(\"y=\", m,\"x\", \"+\", b)\n","\n","\n","def autoFit2(x=[], y=[], dy=[], title=\"Use title= in your call\", xaxis=\"Use xaxis= in your call\", yaxis=\"Use yaxis= in your call\"):\n","    x0=[0,0]\n","    x=np.array(x)\n","    y=np.array(y)\n","    dy=np.array(dy)\n","\n","    args,cov=curve_fit(linearFit, x**2, y, x0, dy)\n","\n","    m=args[0]\n","    b=args[1]\n","    dm = np.sqrt(cov[0][0])\n","    db = np.sqrt(cov[1][1])\n","\n","    fig, ax=plt.subplots(1,2, figsize=(12,6), dpi=75)\n","    line,=ax[0].plot(x,linear(x**2, [m,b]))\n","    data=ax[0].errorbar(x,y, dy, fmt='.k')\n","    ax[0].set_title(title)\n","    ax[0].set_xlabel(xaxis)\n","    ax[0].set_ylabel(yaxis)\n","    ax[1].set_title(\"Residuals\")\n","    ax[1].set_xlabel(xaxis)\n","    ax[1].set_ylabel(\"point -- line\")\n","    residuals=y-linear(x**2, [m,b])\n","    res=ax[1].errorbar(x,residuals, dy, fmt='.k')\n","    ax[1].grid(True, which='both')\n","    plt.show()\n","\n","    print(\"chi-squared value =  {}\".format(chiSquared(x**2,y, dy, linear, [m, b]).round(3)))\n","    print(\"Best fit line :  \")\n","    print(\"y=\", round(args[0],4),\"x^2\", \"+\", round(args[1],4))\n","    print(\"m={}+\\-{}\".format(m.round(3),dm.round(3),))\n","    print(\"b={}+\\-{}\".format(b.round(3),db.round(3),))\n"],"metadata":{"id":"Wd11VBSYpRZf","executionInfo":{"status":"ok","timestamp":1757529562813,"user_tz":240,"elapsed":76,"user":{"displayName":"Benjamin Vaughan","userId":"12757759398363721665"}},"outputId":"46056028-f5b3-4925-93a0-52d2fc1402e1","colab":{"base_uri":"https://localhost:8080/"}},"execution_count":1,"outputs":[{"output_type":"stream","name":"stderr","text":["<>:68: SyntaxWarning: invalid escape sequence '\\-'\n","<>:69: SyntaxWarning: invalid escape sequence '\\-'\n","<>:126: SyntaxWarning: invalid escape sequence '\\-'\n","<>:127: SyntaxWarning: invalid escape sequence '\\-'\n","<>:68: SyntaxWarning: invalid escape sequence '\\-'\n","<>:69: SyntaxWarning: invalid escape sequence '\\-'\n","<>:126: SyntaxWarning: invalid escape sequence '\\-'\n","<>:127: SyntaxWarning: invalid escape sequence '\\-'\n","/tmp/ipython-input-2179927003.py:68: SyntaxWarning: invalid escape sequence '\\-'\n","  print(\"m={}+\\-{}\".format(m.round(3),dm.round(3),))\n","/tmp/ipython-input-2179927003.py:69: SyntaxWarning: invalid escape sequence '\\-'\n","  print(\"b={}+\\-{}\".format(b.round(3),db.round(3),))\n","/tmp/ipython-input-2179927003.py:126: SyntaxWarning: invalid escape sequence '\\-'\n","  print(\"m={}+\\-{}\".format(m.round(3),dm.round(3),))\n","/tmp/ipython-input-2179927003.py:127: SyntaxWarning: invalid escape sequence '\\-'\n","  print(\"b={}+\\-{}\".format(b.round(3),db.round(3),))\n"]}]},{"cell_type":"markdown","source":["# Uncertainty and $t^\\prime$"],"metadata":{"id":"UtqHE7iYpUaZ"}},{"cell_type":"code","source":["def standard_deviation(data):\n","    '''\n","    Takes in a one-dimensional numpy array argument, and returns the standard deviation of the dataset.\n","    '''\n","    N = len(data)\n","    return np.sqrt(np.sum((data - np.mean(data))**2)/(N-1))\n","\n","def standard_unc_of_mean(data):\n","    '''\n","    Takes in a one-dimensional numpy array argument, and returns the standard uncertainty of the mean of the dataset.\n","    '''\n","    N = len(data)\n","    return standard_deviation(data) / np.sqrt(N)\n","def t_prime(A, dA, B, dB):\n","    '''\n","    Returns the value of t prime. Requires 4 arguments in this order: measurement A, uncertainty of measurement A, measurement B, uncertainty of measurement B.\n","    '''\n","    return abs(A - B) / np.sqrt(dA**2 + dB**2)\n"],"metadata":{"id":"azT1FSIMpVs3"},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":["# Uncertainty propagation"],"metadata":{"id":"t8qUBUjLpXPL"}},{"cell_type":"code","source":["def dsum(x, dx, y, dy):\n","    return np.sqrt(dx**2+dy**2)\n","\n","def dsub(x, dx, y, dy):\n","    return np.sqrt(dx**2+dy**2)\n","\n","def dmult(x, dx, y, dy):\n","    return x*y*np.sqrt((dx/x)**2+(dy/y)**2)\n","\n","def ddiv(x, dx, y, dy):\n","    return x/y*np.sqrt((dx/x)**2+(dy/y)**2)\n","\n","def dpwr(x, dx, n):\n","    return n*dx*(x**n)/x\n","\n","def dln(x, dx):\n","    return dx/x"],"metadata":{"id":"TmIz9CTgpYzx"},"execution_count":null,"outputs":[]}]}