{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "48023aee",
   "metadata": {},
   "outputs": [],
   "source": [
    "from tadatakit.classes import Experiment\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import scipy.integrate as integrate\n",
    "from lmfit.models import Pearson4Model\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "20f04ced",
   "metadata": {},
   "outputs": [],
   "source": [
    "experiment =  Experiment.from_json(\"hdpe_pp_50_50.json\")\n",
    "df = experiment.get_dataframe('processed')\n",
    "\n",
    "step_df = df[df['Procedure Step Id']==experiment.procedure.steps[1].id]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "101cb495",
   "metadata": {},
   "outputs": [],
   "source": [
    "peaks_df = step_df[step_df['Temperature / °C'].between(0, 185)].copy()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "a4647dda",
   "metadata": {},
   "outputs": [],
   "source": [
    "def calculate_linear_baseline(peak_data):\n",
    "    temperature = peak_data[\"Temperature / °C\"]\n",
    "    heat_flow = peak_data[\"Heat Flow (Normalized) / W/g\"]\n",
    "    baseline = np.interp(temperature, [temperature.iloc[0], temperature.iloc[-1]], [heat_flow.iloc[0], heat_flow.iloc[-1]])\n",
    "    return baseline\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "7993b335",
   "metadata": {},
   "outputs": [],
   "source": [
    "peaks_df['Baseline Heat Flow (Normalized) / W/g'] = calculate_linear_baseline(peaks_df)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "6a852875",
   "metadata": {},
   "outputs": [],
   "source": [
    "peaks_df['Baseline Subtracted Heat Flow (Normalized) / W/g'] = peaks_df['Heat Flow (Normalized) / W/g'] - peaks_df['Baseline Heat Flow (Normalized) / W/g'] \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "f0918981",
   "metadata": {},
   "outputs": [],
   "source": [
    "p41 = Pearson4Model(prefix='p41_')\n",
    "p42 = Pearson4Model(prefix='p42_')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "428933af",
   "metadata": {},
   "outputs": [],
   "source": [
    "p1=p41.make_params()\n",
    "p1['p41_center'].set(133)\n",
    "p1['p41_amplitude'].set(-1.2)\n",
    "p2=p42.make_params()\n",
    "p2['p42_center'].set(162)\n",
    "p2['p42_amplitude'].set(-0.8)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "2ff6d2bb",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table class=\"jp-toc-ignore\"><caption>Parameters</caption><tr><th style='text-align:left'>name</th><th style='text-align:left'>value</th><th style='text-align:left'>initial value</th><th style='text-align:left'>min</th><th style='text-align:left'>max</th><th style='text-align:left'>vary</th><th style='text-align:right'>expression</th></tr><tr><td style='text-align:left'>p41_amplitude</td><td style='text-align:left'>-1.20000000</td><td style='text-align:left'>-1.2</td><td style='text-align:left'>       -inf</td><td style='text-align:left'>        inf</td><td style='text-align:left'>True</td><td style='text-align:right'></td></tr><tr><td style='text-align:left'>p41_center</td><td style='text-align:left'> 133.000000</td><td style='text-align:left'>133</td><td style='text-align:left'>       -inf</td><td style='text-align:left'>        inf</td><td style='text-align:left'>True</td><td style='text-align:right'></td></tr><tr><td style='text-align:left'>p41_sigma</td><td style='text-align:left'> 1.00000000</td><td style='text-align:left'>None</td><td style='text-align:left'>       -inf</td><td style='text-align:left'>        inf</td><td style='text-align:left'>True</td><td style='text-align:right'></td></tr><tr><td style='text-align:left'>p41_expon</td><td style='text-align:left'> 1.50000000</td><td style='text-align:left'>None</td><td style='text-align:left'> 0.50000000</td><td style='text-align:left'> 1000.00000</td><td style='text-align:left'>True</td><td style='text-align:right'></td></tr><tr><td style='text-align:left'>p41_skew</td><td style='text-align:left'> 0.00000000</td><td style='text-align:left'>None</td><td style='text-align:left'>-1000.00000</td><td style='text-align:left'> 1000.00000</td><td style='text-align:left'>True</td><td style='text-align:right'></td></tr><tr><td style='text-align:left'>p41_fwhm</td><td style='text-align:left'> 1.53284187</td><td style='text-align:left'>None</td><td style='text-align:left'>       -inf</td><td style='text-align:left'>        inf</td><td style='text-align:left'>False</td><td style='text-align:right'>p41_sigma*sqrt(2**(1/p41_expon)-1)*pi/arctan2(exp(1)*p41_expon, p41_skew)</td></tr><tr><td style='text-align:left'>p41_height</td><td style='text-align:left'>-0.60000000</td><td style='text-align:left'>None</td><td style='text-align:left'>       -inf</td><td style='text-align:left'>        inf</td><td style='text-align:left'>False</td><td style='text-align:right'>(p41_amplitude / p41_sigma) * exp(2 * (real(loggammafcn(p41_expon + p41_skew * 0.5j)) - loggammafcn(p41_expon)) - betalnfnc(p41_expon-0.5, 0.5) - p41_expon * log1p(square(p41_skew/(2*p41_expon))) - p41_skew * arctan(-p41_skew/(2*p41_expon)))</td></tr><tr><td style='text-align:left'>p41_position</td><td style='text-align:left'> 133.000000</td><td style='text-align:left'>None</td><td style='text-align:left'>       -inf</td><td style='text-align:left'>        inf</td><td style='text-align:left'>False</td><td style='text-align:right'>p41_center-p41_sigma*p41_skew/(2*p41_expon)</td></tr><tr><td style='text-align:left'>p42_amplitude</td><td style='text-align:left'>-0.80000000</td><td style='text-align:left'>-0.8</td><td style='text-align:left'>       -inf</td><td style='text-align:left'>        inf</td><td style='text-align:left'>True</td><td style='text-align:right'></td></tr><tr><td style='text-align:left'>p42_center</td><td style='text-align:left'> 162.000000</td><td style='text-align:left'>162</td><td style='text-align:left'>       -inf</td><td style='text-align:left'>        inf</td><td style='text-align:left'>True</td><td style='text-align:right'></td></tr><tr><td style='text-align:left'>p42_sigma</td><td style='text-align:left'> 1.00000000</td><td style='text-align:left'>None</td><td style='text-align:left'>       -inf</td><td style='text-align:left'>        inf</td><td style='text-align:left'>True</td><td style='text-align:right'></td></tr><tr><td style='text-align:left'>p42_expon</td><td style='text-align:left'> 1.50000000</td><td style='text-align:left'>None</td><td style='text-align:left'> 0.50000000</td><td style='text-align:left'> 1000.00000</td><td style='text-align:left'>True</td><td style='text-align:right'></td></tr><tr><td style='text-align:left'>p42_skew</td><td style='text-align:left'> 0.00000000</td><td style='text-align:left'>None</td><td style='text-align:left'>-1000.00000</td><td style='text-align:left'> 1000.00000</td><td style='text-align:left'>True</td><td style='text-align:right'></td></tr><tr><td style='text-align:left'>p42_fwhm</td><td style='text-align:left'> 1.53284187</td><td style='text-align:left'>None</td><td style='text-align:left'>       -inf</td><td style='text-align:left'>        inf</td><td style='text-align:left'>False</td><td style='text-align:right'>p42_sigma*sqrt(2**(1/p42_expon)-1)*pi/arctan2(exp(1)*p42_expon, p42_skew)</td></tr><tr><td style='text-align:left'>p42_height</td><td style='text-align:left'>-0.40000000</td><td style='text-align:left'>None</td><td style='text-align:left'>       -inf</td><td style='text-align:left'>        inf</td><td style='text-align:left'>False</td><td style='text-align:right'>(p42_amplitude / p42_sigma) * exp(2 * (real(loggammafcn(p42_expon + p42_skew * 0.5j)) - loggammafcn(p42_expon)) - betalnfnc(p42_expon-0.5, 0.5) - p42_expon * log1p(square(p42_skew/(2*p42_expon))) - p42_skew * arctan(-p42_skew/(2*p42_expon)))</td></tr><tr><td style='text-align:left'>p42_position</td><td style='text-align:left'> 162.000000</td><td style='text-align:left'>None</td><td style='text-align:left'>       -inf</td><td style='text-align:left'>        inf</td><td style='text-align:left'>False</td><td style='text-align:right'>p42_center-p42_sigma*p42_skew/(2*p42_expon)</td></tr></table>"
      ],
      "text/plain": [
       "Parameters([('p41_amplitude', <Parameter 'p41_amplitude', value=-1.2, bounds=[-inf:inf]>), ('p41_center', <Parameter 'p41_center', value=133, bounds=[-inf:inf]>), ('p41_sigma', <Parameter 'p41_sigma', value=1.0, bounds=[-inf:inf]>), ('p41_expon', <Parameter 'p41_expon', value=1.5, bounds=[0.500000000000001:1000]>), ('p41_skew', <Parameter 'p41_skew', value=0.0, bounds=[-1000:1000]>), ('p41_fwhm', <Parameter 'p41_fwhm', value=1.5328418730817595, bounds=[-inf:inf], expr='p41_sigma*sqrt(2**(1/p41_expon)-1)*pi/arctan2(exp(1)*p41_expon, p41_skew)'>), ('p41_height', <Parameter 'p41_height', value=-0.5999999999999982, bounds=[-inf:inf], expr='(p41_amplitude / p41_sigma) * exp(2 * (real(loggammafcn(p41_expon + p41_skew * 0.5j)) - loggammafcn(p41_expon)) - betalnfnc(p41_expon-0.5, 0.5) - p41_expon * log1p(square(p41_skew/(2*p41_expon))) - p41_skew * arctan(-p41_skew/(2*p41_expon)))'>), ('p41_position', <Parameter 'p41_position', value=133.0, bounds=[-inf:inf], expr='p41_center-p41_sigma*p41_skew/(2*p41_expon)'>), ('p42_amplitude', <Parameter 'p42_amplitude', value=-0.8, bounds=[-inf:inf]>), ('p42_center', <Parameter 'p42_center', value=162, bounds=[-inf:inf]>), ('p42_sigma', <Parameter 'p42_sigma', value=1.0, bounds=[-inf:inf]>), ('p42_expon', <Parameter 'p42_expon', value=1.5, bounds=[0.500000000000001:1000]>), ('p42_skew', <Parameter 'p42_skew', value=0.0, bounds=[-1000:1000]>), ('p42_fwhm', <Parameter 'p42_fwhm', value=1.5328418730817595, bounds=[-inf:inf], expr='p42_sigma*sqrt(2**(1/p42_expon)-1)*pi/arctan2(exp(1)*p42_expon, p42_skew)'>), ('p42_height', <Parameter 'p42_height', value=-0.39999999999999886, bounds=[-inf:inf], expr='(p42_amplitude / p42_sigma) * exp(2 * (real(loggammafcn(p42_expon + p42_skew * 0.5j)) - loggammafcn(p42_expon)) - betalnfnc(p42_expon-0.5, 0.5) - p42_expon * log1p(square(p42_skew/(2*p42_expon))) - p42_skew * arctan(-p42_skew/(2*p42_expon)))'>), ('p42_position', <Parameter 'p42_position', value=162.0, bounds=[-inf:inf], expr='p42_center-p42_sigma*p42_skew/(2*p42_expon)'>)])"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model = p41 + p42\n",
    "params=model.make_params()\n",
    "params.update(p1)\n",
    "params.update(p2)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "1e45424e",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 50168 characters]",
      "text/plain": [
       "<Figure size 1000x600 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "init=model.eval(params=params,x=peaks_df[\"Temperature / °C\"])\n",
    "peak1=p41.eval(params=params,x=peaks_df[\"Temperature / °C\"])\n",
    "peak2=p42.eval(params=params,x=peaks_df[\"Temperature / °C\"])\n",
    "\n",
    "fig, ax1 = plt.subplots(figsize=(10,6))\n",
    " \n",
    "ax1.plot(peaks_df[\"Temperature / °C\"], peaks_df['Baseline Subtracted Heat Flow (Normalized) / W/g'], label = 'Original Data')\n",
    "ax1.plot(peaks_df[\"Temperature / °C\"], peak1, label = 'Peak 1 Fit (Initial Guess)')\n",
    "ax1.plot(peaks_df[\"Temperature / °C\"], peak2, label = 'Peak 2 Fit (Initial Guess)')\n",
    " \n",
    "plt.legend()\n",
    "plt.show() \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "e4db9c51",
   "metadata": {},
   "outputs": [],
   "source": [
    "result = model.fit(data=peaks_df['Baseline Subtracted Heat Flow (Normalized) / W/g'], params=params,x=peaks_df[\"Temperature / °C\"])\n",
    "comps = result.eval_components()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "352a363b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 41712 characters]",
      "text/plain": [
       "<Figure size 1000x600 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax1 = plt.subplots(figsize=(10,6))\n",
    "\n",
    "ax1.plot(peaks_df[\"Temperature / °C\"], peaks_df['Baseline Subtracted Heat Flow (Normalized) / W/g'], label = 'Original Data')\n",
    "ax1.plot(peaks_df[\"Temperature / °C\"], result.best_fit, label = 'Fit')\n",
    "\n",
    "plt.legend()\n",
    "plt.show() \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "beba99d1",
   "metadata": {},
   "outputs": [],
   "source": [
    "final_df = peaks_df[['Step time / min', \n",
    "                     'Temperature / °C', \n",
    "                     'Heat Flow (Normalized) / W/g', \n",
    "                     'Baseline Subtracted Heat Flow (Normalized) / W/g', \n",
    "                     'Baseline Heat Flow (Normalized) / W/g']\n",
    "                     ].copy()\n",
    "\n",
    "final_df = final_df.reset_index(drop=True)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "bdc1ebe2",
   "metadata": {},
   "outputs": [],
   "source": [
    "melt1_df = pd.DataFrame.from_records(comps)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "6a640335",
   "metadata": {},
   "outputs": [],
   "source": [
    "final_df['Fit / W/g'] = result.best_fit + final_df['Baseline Heat Flow (Normalized) / W/g']\n",
    "final_df['Melt 1 / W/g'] = melt1_df['p41_'] + final_df['Baseline Heat Flow (Normalized) / W/g']\n",
    "final_df['Melt 2 / W/g'] = melt1_df['p42_'] + final_df['Baseline Heat Flow (Normalized) / W/g']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "a2734367",
   "metadata": {},
   "outputs": [],
   "source": [
    "x = final_df['Step time / min']\n",
    "y1 = final_df['Melt 1 / W/g']\n",
    "y2 = final_df['Melt 2 / W/g']\n",
    " \n",
    "baseline1 = (np.interp(x, [x.iloc[0], x.iloc[-1]],[y1.iloc[0], y1.iloc[-1]]))\n",
    "baseline2 = (np.interp(x, [x.iloc[0], x.iloc[-1]],[y2.iloc[0], y2.iloc[-1]]))\n",
    "Enthalpy1 = (integrate.trapezoid(y1 - baseline1, x))*60\n",
    "Enthalpy2 = (integrate.trapezoid(y2 - baseline2, x))*60 \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "24ebc392",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 58912 characters]",
      "text/plain": [
       "<Figure size 1000x600 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax1 = plt.subplots(figsize=(10,6))\n",
    " \n",
    "ax1.plot(final_df[\"Temperature / °C\"],final_df['Heat Flow (Normalized) / W/g'], label='Original Data')\n",
    "ax1.plot(final_df[\"Temperature / °C\"],final_df['Melt 1 / W/g'], label=\"Melt 1 \" f\"Enthalpy = {Enthalpy1:.2f} J/g\")\n",
    "ax1.plot(final_df[\"Temperature / °C\"],final_df['Melt 2 / W/g'], label=\"Melt 2 \" f\"Enthalpy = {Enthalpy2:.2f} J/g\")\n",
    "\n",
    "plt.legend()\n",
    "plt.show() \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0e6bcae8",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.12.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
