ODEs-to-python.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "import re"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {},
   "outputs": [],
   "source": [
    "def prep_lines(lines):\n",
    "    '''\n",
    "    Remove Fortran-specific code, and remove \\n from reading in as strings. \n",
    "    '''\n",
    "    for i in range(len(lines)):\n",
    "        lines[i]=lines[i].replace('     * ', '')\n",
    "        lines[i]=lines[i].replace(' \\n', '')\n",
    "        lines[i]=lines[i].replace('\\n', '')\n",
    "    return lines\n",
    "\n",
    "\n",
    "def read_fortran_ODEs(file, cons = False):\n",
    "    '''\n",
    "    Read fortran ODEs and prep for translation.\\n\n",
    "        - file = file name\\n\n",
    "        - returns: dict\n",
    "    '''\n",
    "\n",
    "    ## read fortran ODEs from given file\n",
    "    with open(file, 'r') as f:\n",
    "        lines = f.readlines()\n",
    "\n",
    "    if cons == True:\n",
    "        lines = lines[16:93]\n",
    "    if cons == False:\n",
    "        ## select non-conserved species\n",
    "        lines = lines[123:]\n",
    "\n",
    "    ## remove unnecessary fortran code\n",
    "    lines = prep_lines(lines)\n",
    "\n",
    "    ODEs = stick_fortran(lines)\n",
    "\n",
    "    return ODEs\n",
    "\n",
    "\n",
    "def stick_fortran(lines):\n",
    "    '''\n",
    "    Make dictionary:\\n\n",
    "        - entry per species\\n\n",
    "        - stick together separate fortran strings (result: full ODE in 1 string) \\n\n",
    "        \\n\n",
    "    This is possible by searching for 'C' in Fortran, indicating commented code.\n",
    "    '''\n",
    "\n",
    "    ODEs = dict()\n",
    "\n",
    "    for i, line in enumerate(lines):\n",
    "        if len(line) != 0 and line[0] == 'C':\n",
    "            ODEs[line[2:]] = str('')            ## set species name (every that comes after the 'C')\n",
    "            print('i=',i, end = '\\r')\n",
    "            start = i + 1                       ## start searching where this ODE ends\n",
    "\n",
    "            for j in range(start,start+1000):\n",
    "                if j >= len(lines):\n",
    "                    break\n",
    "                if len(lines[j]) != 0 and lines[j][0] == 'C':   ## stop searching when the next comment 'C' is found\n",
    "                    stop = j\n",
    "                    break\n",
    "                \n",
    "            for j in range(i+1,stop):           ## stick together the different pieces of fortran ODE\n",
    "                ODEs[line[2:]] += lines[j]\n",
    "\n",
    "            i = stop + 1                        ## continue after the ODE\n",
    "\n",
    "    count = 0\n",
    "    for key in ODEs:            ## count empty keys in dict()\n",
    "        if key == '':\n",
    "            count += 1\n",
    "\n",
    "    while count > 0:            ## remove empty keys in dict()\n",
    "        del ODEs['']\n",
    "        count -= 1\n",
    "\n",
    "    return ODEs\n",
    "\n",
    "\n",
    "\n",
    "def fortran_to_python(ODEs, spec, nblen, nbs, cons = False, ODEtype = 'ipv'):\n",
    "    '''\n",
    "    Translate fortran ODE to python ODE. \\n\\n\n",
    "\n",
    "    Note: NOT fractional abundances, but in units of cm^-3.\n",
    "    '''\n",
    "    \n",
    "    print(spec)\n",
    "\n",
    "    if isinstance(ODEs[spec],str):\n",
    "        ## split string\n",
    "        ODEs[spec] = ODEs[spec].split(' ')\n",
    "\n",
    "    ## remove empty strings\n",
    "    to_remove = ''\n",
    "    while to_remove in ODEs[spec]:\n",
    "        ODEs[spec].remove('')\n",
    "\n",
    "    ## remove other fortran commands\n",
    "    while 'END' in ODEs[spec]:\n",
    "        ODEs[spec].remove('END')\n",
    "    while 'RETURN' in ODEs[spec]:\n",
    "        ODEs[spec].remove('RETURN')\n",
    "\n",
    "    \n",
    "    ## Translate\n",
    "\n",
    "    ## Conserved species\n",
    "    if cons == True:\n",
    "        ## replace brackets: ( --> [\n",
    "        print(len(ODEs[spec]))\n",
    "        for i in range(len(ODEs[spec])):\n",
    "            \n",
    "            ODEs[spec][i] = ODEs[spec][i].replace('(', '[')\n",
    "            ODEs[spec][i] = ODEs[spec][i].replace(')', ']')\n",
    "            ODEs[spec][i] = ODEs[spec][i].replace('[+', '(+')\n",
    "            ODEs[spec][i] = ODEs[spec][i].replace('[-', '(-')\n",
    "            ODEs[spec][i] = ODEs[spec][i].replace(']]', '])')\n",
    "\n",
    "            if ODEtype == 'torch':\n",
    "                line = list(ODEs[spec][i])\n",
    "                for j in range(len(line)):\n",
    "                    if line[j] == 'Y':\n",
    "                        line[j+1] = '[:,'\n",
    "\n",
    "                ODEs[spec][i] = \"\".join(line)\n",
    "\n",
    "    ## Non-conserved species\n",
    "    if cons == False:\n",
    "\n",
    "        ## prep for replacing final line --> sets the change in number density of species\n",
    "        last = ODEs[spec][-1]\n",
    "        while not last[0] == 'Y':       ## remove fortran statements & empty lines\n",
    "            ODEs[spec].remove(last)\n",
    "            last = ODEs[spec][-1]\n",
    "        ## get reaction number\n",
    "        for i in range(5,len(last)):    \n",
    "            if last[i] == ')':\n",
    "                stop = i\n",
    "                break\n",
    "        nb = last[5:stop]\n",
    "\n",
    "        nblen[spec] = len(nb)\n",
    "        nbs[spec] = int(nb)\n",
    "\n",
    "        ## replace brackets: ( --> [\n",
    "        for i in range(len(ODEs[spec])):\n",
    "            if ODEs[spec][i][0] == 'F' or ODEs[spec][i][0] == 'D':\n",
    "                ODEs[spec][i] = ODEs[spec][i].replace('(', '[')\n",
    "                ODEs[spec][i] = ODEs[spec][i].replace(')', ']')\n",
    "                ODEs[spec][i] = ODEs[spec][i].replace('F', 'F['+nb+']')\n",
    "                ODEs[spec][i] = ODEs[spec][i].replace('D', 'D['+nb+']')\n",
    "\n",
    "                if ODEtype == 'torch':\n",
    "                    line = list(ODEs[spec][i])\n",
    "                    for j in range(len(line)):\n",
    "                        if line[j] == 'Y' or line[j] == 'F' or line[j] == 'D':\n",
    "                            line[j+1] = '[:,'\n",
    "\n",
    "                    ODEs[spec][i] = \"\".join(line)\n",
    "\n",
    "        if spec == 'H':\n",
    "            ODEs[spec][-1] = 'F['+nb+']=F['+nb+']+HLOSS'\n",
    "        else:\n",
    "            ODEs[spec][-1] = ''\n",
    "             \n",
    "    return ODEs\n",
    "\n",
    "\n",
    "def parse_line(parsed, spec, nblen, cons = False, ODEtype = 'ivp'):\n",
    "    '''\n",
    "    Parse line to print is too long (Python doesn't like that much).\n",
    "    '''\n",
    "    \n",
    "    if cons == False:\n",
    "        line = parsed[-1]\n",
    "        \n",
    "        if len(line) > 1000:\n",
    "            if cons == False:\n",
    "                for i in range(1000):\n",
    "                    n=i+1000\n",
    "                    if n < len(line):\n",
    "                        if line[n] == '+' or line[n] == '-':\n",
    "                            break\n",
    "                    else:\n",
    "                        break\n",
    "                parsed[-1] = line[0:n]\n",
    "\n",
    "                length = nblen[spec]\n",
    "                \n",
    "                if ODEtype == 'ivp':     l =  3\n",
    "                if ODEtype == 'torch':   l =  5\n",
    "\n",
    "                parsed.append(line[0:length+l]+'='+line[0:length+l]+line[n:])            \n",
    "\n",
    "        else: \n",
    "            print('No parsing need')\n",
    "    \n",
    "    return \n",
    "\n",
    "def minus_one(match):\n",
    "    '''\n",
    "    Helper function that is used in the regular expression substitution.\n",
    "    '''\n",
    "    # Extract the mached string from the matched object\n",
    "    match = match.group(0)\n",
    "    # remove square brackets\n",
    "    match = match[1:-1]\n",
    "    # Cast to integer\n",
    "    match = int(match)    \n",
    "    # Minus one\n",
    "    match -= 1\n",
    "    # Add brackets again\n",
    "    match = f'[{match}]'\n",
    "\n",
    "    return match\n",
    "\n",
    "def minus_one_torch(match):\n",
    "    '''\n",
    "    Helper function that is used in the regular expression substitution.\n",
    "    '''\n",
    "    # Extract the mached string from the matched object\n",
    "    match = match.group(0)\n",
    "    # remove square brackets\n",
    "    match = match[3:-1]\n",
    "    # Cast to integer\n",
    "    match = int(match)    \n",
    "    # Minus one\n",
    "    match -= 1\n",
    "    # Add brackets again\n",
    "    match = f'[:,{match}]'\n",
    "\n",
    "    return match"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Read in file & prep"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "i= 10385\r"
     ]
    }
   ],
   "source": [
    "# file = 'src/ode/acodes.f'\n",
    "loc = 'src/ode/'\n",
    "filename = 'acodes'\n",
    "\n",
    "ODEtype = 'torch'\n",
    "\n",
    "with open(loc+filename+'.f', 'r') as f:\n",
    "    lines = f.readlines()\n",
    "\n",
    "ODEs = read_fortran_ODEs(loc+filename+'.f')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Translate"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Non-conserved species"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "H\n",
      "H+\n",
      "H-\n",
      "H2+\n",
      "H3+\n",
      "He\n",
      "He+\n",
      "HeH+\n",
      "C-\n",
      "C+\n",
      "C\n",
      "CH-\n",
      "CH+\n",
      "CH\n",
      "CH2+\n",
      "N+\n",
      "CH2\n",
      "N\n",
      "CH3\n",
      "NH\n",
      "NH+\n",
      "CH3+\n",
      "O+\n",
      "NH2+\n",
      "CH4+\n",
      "O-\n",
      "NH2\n",
      "O\n",
      "CH4\n",
      "OH-\n",
      "OH\n",
      "OH+\n",
      "CH5+\n",
      "NH3+\n",
      "NH3\n",
      "NH4+\n",
      "H2O\n",
      "H2O+\n",
      "F\n",
      "H3O+\n",
      "F+\n",
      "HF\n",
      "HF+\n",
      "H2F+\n",
      "Na\n",
      "Na+\n",
      "C2-\n",
      "C2+\n",
      "Mg\n",
      "Mg+\n",
      "C2\n",
      "C2H-\n",
      "C2H+\n",
      "C2H\n",
      "CN\n",
      "C2H2\n",
      "C2H2+\n",
      "CN-\n",
      "CN+\n",
      "HCN\n",
      "C2H3\n",
      "C2H3+\n",
      "HCN+\n",
      "HNC\n",
      "N2\n",
      "C2H4+\n",
      "Si+\n",
      "H2CN\n",
      "H2NC+\n",
      "N2+\n",
      "Si\n",
      "CO\n",
      "CO+\n",
      "C2H4\n",
      "HCNH+\n",
      "C2H5+\n",
      "HCO\n",
      "SiH+\n",
      "HCO+\n",
      "CH2NH\n",
      "SiH\n",
      "C2H5\n",
      "HOC+\n",
      "N2H+\n",
      "CH2NH2+\n",
      "SiH2\n",
      "CH3CH3+\n",
      "NO\n",
      "H2CO\n",
      "CH4N+\n",
      "H2CO+\n",
      "CH3CH3\n",
      "SiH2+\n",
      "NO+\n",
      "P+\n",
      "P\n",
      "H3CO+\n",
      "HNO+\n",
      "SiH3+\n",
      "HNO\n",
      "CF+\n",
      "C2H7+\n",
      "SiH3\n",
      "O2\n",
      "SiH4\n",
      "PH\n",
      "S\n",
      "H2NO+\n",
      "CH3OH\n",
      "SiH4+\n",
      "O2+\n",
      "S-\n",
      "O2-\n",
      "PH+\n",
      "CH3OH+\n",
      "S+\n",
      "CH3OH2+\n",
      "O2H\n",
      "PH2+\n",
      "HS\n",
      "HS+\n",
      "PH2\n",
      "O2H+\n",
      "SiH5+\n",
      "H2O2\n",
      "H2S+\n",
      "H2S\n",
      "PH3+\n",
      "H3S+\n",
      "Cl\n",
      "Cl+\n",
      "C3+\n",
      "HCl\n",
      "C3-\n",
      "C3\n",
      "HCl+\n",
      "H2Cl+\n",
      "C3H-\n",
      "C3H+\n",
      "C3H\n",
      "C2N\n",
      "C3H2+\n",
      "C2N+\n",
      "H2CCC\n",
      "CNC+\n",
      "C3H2\n",
      "C2NH+\n",
      "C3H3+\n",
      "CH2CCH\n",
      "CH2CCH+\n",
      "SiC\n",
      "CH2CCH2\n",
      "SiC+\n",
      "C2O\n",
      "C3H4+\n",
      "C2O+\n",
      "CH2CN\n",
      "CH3CCH\n",
      "CH2CN+\n",
      "C3H5+\n",
      "HCSi\n",
      "HCSi+\n",
      "CH3CN+\n",
      "HC2O+\n",
      "CH3CN\n",
      "C3H6+\n",
      "SiN\n",
      "CH3CNH+\n",
      "NH2CN\n",
      "CH3CHCH2\n",
      "SiN+\n",
      "SiCH2+\n",
      "SiCH2\n",
      "CNO\n",
      "CH2CO+\n",
      "OCN\n",
      "OCN+\n",
      "CH2CO\n",
      "CP+\n",
      "HNCO\n",
      "HCNO+\n",
      "HNCO+\n",
      "HONC\n",
      "SiCH3+\n",
      "HCNO\n",
      "CP\n",
      "C3H7+\n",
      "NH2CNH+\n",
      "HNSi+\n",
      "HOCN\n",
      "HNSi\n",
      "CH3CO+\n",
      "HOCN+\n",
      "SiCH3\n",
      "HONC+\n",
      "CS+\n",
      "HCNOH+\n",
      "N2O+\n",
      "CS\n",
      "SiCH4+\n",
      "CO2+\n",
      "H2OCN+\n",
      "HNCOH+\n",
      "SiNH2+\n",
      "HCP\n",
      "SiO+\n",
      "CH3CHO\n",
      "SiO\n",
      "H2CNO+\n",
      "CH3CHO+\n",
      "N2O\n",
      "H2NCO+\n",
      "HCP+\n",
      "CO2\n",
      "CH3CHOH+\n",
      "HCS+\n",
      "SiOH+\n",
      "PN+\n",
      "PN\n",
      "HCS\n",
      "PCH2+\n",
      "HN2O+\n",
      "HCO2+\n",
      "NS+\n",
      "H2SiO\n",
      "H2CS\n",
      "NS\n",
      "CH2PH\n",
      "HCOOH+\n",
      "C2H5OH+\n",
      "NO2\n",
      "CH3OCH3\n",
      "H2CS+\n",
      "HCOOH\n",
      "PCH3+\n",
      "H2SiO+\n",
      "CH3OCH3+\n",
      "C2H5OH\n",
      "NO2+\n",
      "HPN+\n",
      "CH3OCH4+\n",
      "SiF+\n",
      "PNH2+\n",
      "CCl\n",
      "H3CS+\n",
      "PO+\n",
      "CCl+\n",
      "HNS+\n",
      "PCH4+\n",
      "HCOOH2+\n",
      "PO\n",
      "C2H5OH2+\n",
      "H3SiO+\n",
      "C4-\n",
      "HPO+\n",
      "HPO\n",
      "C4\n",
      "SO\n",
      "C4+\n",
      "SO+\n",
      "PNH3+\n",
      "C4H+\n",
      "H2CCl+\n",
      "H2PO+\n",
      "C4H\n",
      "C4H-\n",
      "HSO+\n",
      "C4H2+\n",
      "C3N\n",
      "HC4H\n",
      "C3N+\n",
      "C3N-\n",
      "HC3N+\n",
      "C4H3\n",
      "HC3N\n",
      "ClO\n",
      "ClO+\n",
      "C4H3+\n",
      "HNC3\n",
      "C3O\n",
      "CH2CHCCH\n",
      "SiC2\n",
      "SiC2+\n",
      "NCCN\n",
      "HC3NH+\n",
      "C4H4+\n",
      "C2N2+\n",
      "C3O+\n",
      "HC3O+\n",
      "C4H5+\n",
      "SiC2H\n",
      "CH2CHCN\n",
      "CH2CHCN+\n",
      "SiC2H+\n",
      "NCCNH+\n",
      "SiC2H2+\n",
      "CH2CHCNH+\n",
      "CH2CHCHCH2\n",
      "C3H2O+\n",
      "SiC2H2\n",
      "SiNC+\n",
      "SiNC\n",
      "SiNCH+\n",
      "CCP\n",
      "C2H5CN\n",
      "SiC2H3+\n",
      "CCP+\n",
      "C4H7+\n",
      "H3C3O+\n",
      "HC2P\n",
      "HC2P+\n",
      "Fe+\n",
      "C2H5CNH+\n",
      "C2S\n",
      "C2S+\n",
      "Fe\n",
      "HC2S+\n",
      "PC2H2+\n",
      "CH3COCH3+\n",
      "CH3COCH3\n",
      "PC2H3+\n",
      "CH3COCH4+\n",
      "CH3CS+\n",
      "PC2H4+\n",
      "HCOOCH3\n",
      "C5\n",
      "OCS+\n",
      "SiO2\n",
      "C5-\n",
      "SiS+\n",
      "OCS\n",
      "SiS\n",
      "C5+\n",
      "COOCH4+\n",
      "HOCS+\n",
      "HSiS+\n",
      "C5H+\n",
      "HSiO2+\n",
      "C5H\n",
      "H5C2O2+\n",
      "C5H-\n",
      "C5H2+\n",
      "C4N+\n",
      "C5H2\n",
      "C4N\n",
      "HC4N+\n",
      "C5H3+\n",
      "CH3C4H\n",
      "SO2\n",
      "S2\n",
      "SiC3+\n",
      "CH3C4H+\n",
      "S2+\n",
      "SiC3\n",
      "H2C4N+\n",
      "SO2+\n",
      "HS2\n",
      "HSO2+\n",
      "SiC3H+\n",
      "CH3C3N+\n",
      "HS2+\n",
      "CH3C3N\n",
      "SiC3H\n",
      "C5H5+\n",
      "CH3C3NH+\n",
      "H2S2+\n",
      "SiC3H2+\n",
      "H2S2\n",
      "H3S2+\n",
      "NCCNCH3+\n",
      "C3P\n",
      "PC3H+\n",
      "C3S\n",
      "C3S+\n",
      "HC3S+\n",
      "C6-\n",
      "C6+\n",
      "C6\n",
      "C6H\n",
      "C6H+\n",
      "C6H-\n",
      "C6H2\n",
      "C6H2+\n",
      "C5N\n",
      "C5N-\n",
      "C5N+\n",
      "HC5N+\n",
      "HC5N\n",
      "C6H3+\n",
      "SiC4+\n",
      "C6H4+\n",
      "HC5NH+\n",
      "SiC4\n",
      "H3C5N+\n",
      "SiC4H+\n",
      "C6H5+\n",
      "C6H6\n",
      "C6H6+\n",
      "C4P\n",
      "C6H7+\n",
      "C4P+\n",
      "C4S\n",
      "PC4H+\n",
      "C4S+\n",
      "HC4S+\n",
      "C7+\n",
      "C7\n",
      "C7-\n",
      "C7H-\n",
      "C7H+\n",
      "C7H\n",
      "C7H2+\n",
      "C7H2\n",
      "C7H3+\n",
      "C7H4+\n",
      "CH3C6H\n",
      "CH3C5N\n",
      "C7H5+\n",
      "CH3C5NH+\n",
      "C8-\n",
      "C8\n",
      "C8+\n",
      "C8H\n",
      "C8H+\n",
      "C8H-\n",
      "C8H2\n",
      "C7N\n",
      "C7N+\n",
      "C8H2+\n",
      "HC7N\n",
      "HC7N+\n",
      "C8H3+\n",
      "C8H4+\n",
      "H2C7N+\n",
      "C8H5+\n",
      "H3C7N+\n",
      "C9\n",
      "C9+\n",
      "C9-\n",
      "C9H+\n",
      "C9H-\n",
      "C9H\n",
      "C9H2+\n",
      "C9H2\n",
      "C9H3+\n",
      "C9H4+\n",
      "CH3C7N\n",
      "C9H5+\n",
      "CH3C7NH+\n",
      "C10-\n",
      "C10\n",
      "C10+\n",
      "C10H+\n",
      "C10H-\n",
      "C10H\n",
      "C9N\n",
      "C9N+\n",
      "C10H2\n",
      "C10H2+\n",
      "HC9N+\n",
      "C10H3+\n",
      "HC9N\n",
      "H2C9N+\n",
      "H3C9N+\n",
      "C11\n",
      "C11+\n"
     ]
    }
   ],
   "source": [
    "nblen = dict()      ## dictionary to store the length of the reaction number per species\n",
    "nbs = dict()\n",
    "for spec in ODEs:\n",
    "    ODEs = fortran_to_python(ODEs, spec, nblen, nbs, ODEtype=ODEtype)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "H\n",
      "H+\n",
      "H-\n",
      "H2+\n",
      "H3+\n",
      "He\n",
      "He+\n",
      "HeH+\n",
      "C-\n",
      "C+\n",
      "C\n",
      "CH-\n",
      "CH+\n",
      "CH\n",
      "CH2+\n",
      "N+\n",
      "CH2\n",
      "N\n",
      "CH3\n",
      "NH\n",
      "NH+\n",
      "CH3+\n",
      "O+\n",
      "NH2+\n",
      "CH4+\n",
      "O-\n",
      "NH2\n",
      "O\n",
      "CH4\n",
      "OH-\n",
      "OH\n",
      "OH+\n",
      "CH5+\n",
      "NH3+\n",
      "NH3\n",
      "NH4+\n",
      "H2O\n",
      "H2O+\n",
      "F\n",
      "H3O+\n",
      "F+\n",
      "HF\n",
      "HF+\n",
      "H2F+\n",
      "Na\n",
      "Na+\n",
      "C2-\n",
      "C2+\n",
      "Mg\n",
      "Mg+\n",
      "C2\n",
      "C2H-\n",
      "C2H+\n",
      "C2H\n",
      "CN\n",
      "C2H2\n",
      "C2H2+\n",
      "CN-\n",
      "CN+\n",
      "HCN\n",
      "C2H3\n",
      "C2H3+\n",
      "HCN+\n",
      "HNC\n",
      "N2\n",
      "C2H4+\n",
      "Si+\n",
      "H2CN\n",
      "H2NC+\n",
      "N2+\n",
      "Si\n",
      "CO\n",
      "CO+\n",
      "C2H4\n",
      "HCNH+\n",
      "C2H5+\n",
      "HCO\n",
      "SiH+\n",
      "HCO+\n",
      "CH2NH\n",
      "SiH\n",
      "C2H5\n",
      "HOC+\n",
      "N2H+\n",
      "CH2NH2+\n",
      "SiH2\n",
      "CH3CH3+\n",
      "NO\n",
      "H2CO\n",
      "CH4N+\n",
      "H2CO+\n",
      "CH3CH3\n",
      "SiH2+\n",
      "NO+\n",
      "P+\n",
      "P\n",
      "H3CO+\n",
      "HNO+\n",
      "SiH3+\n",
      "HNO\n",
      "CF+\n",
      "C2H7+\n",
      "SiH3\n",
      "O2\n",
      "SiH4\n",
      "PH\n",
      "S\n",
      "H2NO+\n",
      "CH3OH\n",
      "SiH4+\n",
      "O2+\n",
      "S-\n",
      "O2-\n",
      "PH+\n",
      "CH3OH+\n",
      "S+\n",
      "CH3OH2+\n",
      "O2H\n",
      "PH2+\n",
      "HS\n",
      "HS+\n",
      "PH2\n",
      "O2H+\n",
      "SiH5+\n",
      "H2O2\n",
      "H2S+\n",
      "H2S\n",
      "PH3+\n",
      "H3S+\n",
      "Cl\n",
      "Cl+\n",
      "C3+\n",
      "HCl\n",
      "C3-\n",
      "C3\n",
      "HCl+\n",
      "H2Cl+\n",
      "C3H-\n",
      "C3H+\n",
      "C3H\n",
      "C2N\n",
      "C3H2+\n",
      "C2N+\n",
      "H2CCC\n",
      "CNC+\n",
      "C3H2\n",
      "C2NH+\n",
      "C3H3+\n",
      "CH2CCH\n",
      "CH2CCH+\n",
      "SiC\n",
      "CH2CCH2\n",
      "SiC+\n",
      "C2O\n",
      "C3H4+\n",
      "C2O+\n",
      "CH2CN\n",
      "CH3CCH\n",
      "CH2CN+\n",
      "C3H5+\n",
      "HCSi\n",
      "HCSi+\n",
      "CH3CN+\n",
      "HC2O+\n",
      "CH3CN\n",
      "C3H6+\n",
      "SiN\n",
      "CH3CNH+\n",
      "NH2CN\n",
      "CH3CHCH2\n",
      "SiN+\n",
      "SiCH2+\n",
      "SiCH2\n",
      "CNO\n",
      "CH2CO+\n",
      "OCN\n",
      "OCN+\n",
      "CH2CO\n",
      "CP+\n",
      "HNCO\n",
      "HCNO+\n",
      "HNCO+\n",
      "HONC\n",
      "SiCH3+\n",
      "HCNO\n",
      "CP\n",
      "C3H7+\n",
      "NH2CNH+\n",
      "HNSi+\n",
      "HOCN\n",
      "HNSi\n",
      "CH3CO+\n",
      "HOCN+\n",
      "SiCH3\n",
      "HONC+\n",
      "CS+\n",
      "HCNOH+\n",
      "N2O+\n",
      "CS\n",
      "SiCH4+\n",
      "CO2+\n",
      "H2OCN+\n",
      "HNCOH+\n",
      "SiNH2+\n",
      "HCP\n",
      "SiO+\n",
      "CH3CHO\n",
      "SiO\n",
      "H2CNO+\n",
      "CH3CHO+\n",
      "N2O\n",
      "H2NCO+\n",
      "HCP+\n",
      "CO2\n",
      "CH3CHOH+\n",
      "HCS+\n",
      "SiOH+\n",
      "PN+\n",
      "PN\n",
      "HCS\n",
      "PCH2+\n",
      "HN2O+\n",
      "HCO2+\n",
      "NS+\n",
      "H2SiO\n",
      "H2CS\n",
      "NS\n",
      "CH2PH\n",
      "HCOOH+\n",
      "C2H5OH+\n",
      "NO2\n",
      "CH3OCH3\n",
      "H2CS+\n",
      "HCOOH\n",
      "PCH3+\n",
      "H2SiO+\n",
      "CH3OCH3+\n",
      "C2H5OH\n",
      "NO2+\n",
      "HPN+\n",
      "CH3OCH4+\n",
      "SiF+\n",
      "PNH2+\n",
      "CCl\n",
      "H3CS+\n",
      "PO+\n",
      "CCl+\n",
      "HNS+\n",
      "PCH4+\n",
      "HCOOH2+\n",
      "PO\n",
      "C2H5OH2+\n",
      "H3SiO+\n",
      "C4-\n",
      "HPO+\n",
      "HPO\n",
      "C4\n",
      "SO\n",
      "C4+\n",
      "SO+\n",
      "PNH3+\n",
      "C4H+\n",
      "H2CCl+\n",
      "H2PO+\n",
      "C4H\n",
      "C4H-\n",
      "HSO+\n",
      "C4H2+\n",
      "C3N\n",
      "HC4H\n",
      "C3N+\n",
      "C3N-\n",
      "HC3N+\n",
      "C4H3\n",
      "HC3N\n",
      "ClO\n",
      "ClO+\n",
      "C4H3+\n",
      "HNC3\n",
      "C3O\n",
      "CH2CHCCH\n",
      "SiC2\n",
      "SiC2+\n",
      "NCCN\n",
      "HC3NH+\n",
      "C4H4+\n",
      "C2N2+\n",
      "C3O+\n",
      "HC3O+\n",
      "C4H5+\n",
      "SiC2H\n",
      "CH2CHCN\n",
      "CH2CHCN+\n",
      "SiC2H+\n",
      "NCCNH+\n",
      "SiC2H2+\n",
      "CH2CHCNH+\n",
      "CH2CHCHCH2\n",
      "C3H2O+\n",
      "SiC2H2\n",
      "SiNC+\n",
      "SiNC\n",
      "SiNCH+\n",
      "CCP\n",
      "C2H5CN\n",
      "SiC2H3+\n",
      "CCP+\n",
      "C4H7+\n",
      "H3C3O+\n",
      "HC2P\n",
      "HC2P+\n",
      "Fe+\n",
      "C2H5CNH+\n",
      "C2S\n",
      "C2S+\n",
      "Fe\n",
      "HC2S+\n",
      "PC2H2+\n",
      "CH3COCH3+\n",
      "CH3COCH3\n",
      "PC2H3+\n",
      "CH3COCH4+\n",
      "CH3CS+\n",
      "PC2H4+\n",
      "HCOOCH3\n",
      "C5\n",
      "OCS+\n",
      "SiO2\n",
      "C5-\n",
      "SiS+\n",
      "OCS\n",
      "SiS\n",
      "C5+\n",
      "COOCH4+\n",
      "HOCS+\n",
      "HSiS+\n",
      "C5H+\n",
      "HSiO2+\n",
      "C5H\n",
      "H5C2O2+\n",
      "C5H-\n",
      "C5H2+\n",
      "C4N+\n",
      "C5H2\n",
      "C4N\n",
      "HC4N+\n",
      "C5H3+\n",
      "CH3C4H\n",
      "SO2\n",
      "S2\n",
      "SiC3+\n",
      "CH3C4H+\n",
      "S2+\n",
      "SiC3\n",
      "H2C4N+\n",
      "SO2+\n",
      "HS2\n",
      "HSO2+\n",
      "SiC3H+\n",
      "CH3C3N+\n",
      "HS2+\n",
      "CH3C3N\n",
      "SiC3H\n",
      "C5H5+\n",
      "CH3C3NH+\n",
      "H2S2+\n",
      "SiC3H2+\n",
      "H2S2\n",
      "H3S2+\n",
      "NCCNCH3+\n",
      "C3P\n",
      "PC3H+\n",
      "C3S\n",
      "C3S+\n",
      "HC3S+\n",
      "C6-\n",
      "C6+\n",
      "C6\n",
      "C6H\n",
      "C6H+\n",
      "C6H-\n",
      "C6H2\n",
      "C6H2+\n",
      "C5N\n",
      "C5N-\n",
      "C5N+\n",
      "HC5N+\n",
      "HC5N\n",
      "C6H3+\n",
      "SiC4+\n",
      "C6H4+\n",
      "HC5NH+\n",
      "SiC4\n",
      "H3C5N+\n",
      "SiC4H+\n",
      "C6H5+\n",
      "C6H6\n",
      "C6H6+\n",
      "C4P\n",
      "C6H7+\n",
      "C4P+\n",
      "C4S\n",
      "PC4H+\n",
      "C4S+\n",
      "HC4S+\n",
      "C7+\n",
      "C7\n",
      "C7-\n",
      "C7H-\n",
      "C7H+\n",
      "C7H\n",
      "C7H2+\n",
      "C7H2\n",
      "C7H3+\n",
      "C7H4+\n",
      "CH3C6H\n",
      "CH3C5N\n",
      "C7H5+\n",
      "CH3C5NH+\n",
      "C8-\n",
      "C8\n",
      "C8+\n",
      "C8H\n",
      "C8H+\n",
      "C8H-\n",
      "C8H2\n",
      "C7N\n",
      "C7N+\n",
      "C8H2+\n",
      "HC7N\n",
      "HC7N+\n",
      "C8H3+\n",
      "C8H4+\n",
      "H2C7N+\n",
      "C8H5+\n",
      "H3C7N+\n",
      "C9\n",
      "C9+\n",
      "C9-\n",
      "C9H+\n",
      "C9H-\n",
      "C9H\n",
      "C9H2+\n",
      "C9H2\n",
      "C9H3+\n",
      "C9H4+\n",
      "CH3C7N\n",
      "C9H5+\n",
      "CH3C7NH+\n",
      "C10-\n",
      "C10\n",
      "C10+\n",
      "C10H+\n",
      "C10H-\n",
      "C10H\n",
      "C9N\n",
      "C9N+\n",
      "C10H2\n",
      "C10H2+\n",
      "HC9N+\n",
      "C10H3+\n",
      "HC9N\n",
      "H2C9N+\n",
      "H3C9N+\n",
      "C11\n",
      "C11+\n"
     ]
    }
   ],
   "source": [
    "## Parsing the ODEs of the non-conserved specs\n",
    "for spec in ODEs:\n",
    "    print(spec)\n",
    "\n",
    "    i = 0\n",
    "    while i < len(ODEs[spec]):\n",
    "\n",
    "        if len(ODEs[spec][i]) > 1000:\n",
    "            parsed = [ODEs[spec][i]]\n",
    "            while len(parsed[-1]) > 1000:\n",
    "                parse_line(parsed, spec, nblen, ODEtype = ODEtype)\n",
    "\n",
    "            ODEs[spec] = ODEs[spec][0:i] + parsed + ODEs[spec][i+1:]\n",
    "\n",
    "            i = i + len(parsed)\n",
    "        \n",
    "\n",
    "        else: i += 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Conserved species"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "i= 76\r"
     ]
    }
   ],
   "source": [
    "ODEs_cons = read_fortran_ODEs(loc+filename+'.f', True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "e-\n",
      "1\n",
      "H2\n",
      "1\n"
     ]
    }
   ],
   "source": [
    "for spec in ODEs_cons:\n",
    "    ODEs_cons = fortran_to_python(ODEs_cons, spec, nblen, nbs, cons = True, ODEtype=ODEtype)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Write out in as python function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "src/ode/acodes_torch.py\n"
     ]
    }
   ],
   "source": [
    "\n",
    "loc = 'src/ode/'\n",
    "\n",
    "if ODEtype == 'torch':\n",
    "    filename = filename+'_torch'\n",
    "\n",
    "print(loc+filename+'.py')\n",
    "\n",
    "## Y = number density of unconserved species\n",
    "## X = number density of conserved species\n",
    "## K = reaction rate\n",
    "## ACCR = loss of H due to dust interaction, kind of rate\n",
    "\n",
    "with open(loc+filename+'.py', 'w') as f:\n",
    "    ## imports\n",
    "    if ODEtype == 'ivp':\n",
    "        f.write('import numpy as np\\n')\n",
    "    if ODEtype == 'torch':\n",
    "        f.write('import torch\\n')\n",
    "\n",
    "    if filename == 'acodes' or filename =='acodes_torch' :\n",
    "        f.write(\"\\n'''\\nIMPORTANT NOTE:\\nThe content of this function is translated from Fortran 77. \\n \\nWritten by Silke Maes, April 2023.\\nOriginal Fortran code by rate10odecsT.pl on Sun, Jul 21 2019, \\nusing ratefile: rate16_IP_2330K_AP_6000K.rates, from http://udfa.ajmarkwick.net/index.php. \\n'''\\n\")\n",
    "    if filename == 'dcodes' or filename =='dcodes_torch':\n",
    "        f.write(\"\\n'''\\nIMPORTANT NOTE:\\nThe content of this function is translated from Fortran 77. \\n \\nWritten by Silke Maes, April 2023.\\nOriginal Fortran code by rate13dc.pl on Fri, May 31 2013, \\nusing ratefile: rate13.rates, from http://udfa.ajmarkwick.net/index.php. \\n'''\\n\")\n",
    "    f.write(\"\\n\")\n",
    "\n",
    "    ## start of the function\n",
    "    if ODEtype == 'ivp':\n",
    "        f.write('\\ndef ODE(t, Y, YDOT, X, TOTAL, K, HNR, ACCR):\\n')\n",
    "    if ODEtype == 'torch':\n",
    "        f.write('\\ndef torchODE(t, Y, args):\\n')\n",
    "    \n",
    "    ## documentation\n",
    "    f.write(\"\\t'''\")\n",
    "    f.write(\"\\n\\tCalculate the ODEs.\\n\")\n",
    "    f.write(\"\\n\\tY     = number density of the non-conserved species                      == n\")\n",
    "    f.write(\"\\n\\tYDOT  = change in number density of the conserved species                == ndot\")\n",
    "    f.write(\"\\n\\tX     = number density of the conserved species during calculation       == nconsvdot\")\n",
    "    f.write(\"\\n\\tK     = reaction rates, calculated in rates.py                           == k \")\n",
    "    f.write(\"\\n\\tACCR  = H accretion onto dust, calculated in main.py from input physics  == Haccr \")\n",
    "    f.write(\"\\n\\tTOTAL = initial number density of conserved species                      == nconsv\")\n",
    "    f.write(\"\\n\\tHNR   = input density                                                    == ρ\")\n",
    "    f.write(\"\\n\\t'''\\n\")\n",
    "\n",
    "    ## set cut off\n",
    "    f.write(\"\\n\\tcutoff = 1e-50\\n\")\n",
    "    f.write(\"\\n\\tY[Y < cutoff] = 0\\n\")\n",
    "\n",
    "    if ODEtype == 'torch':\n",
    "        f.write(\"\\n\\tYDOT, X, TOTAL, K, HNR, ACCR = args\")\n",
    "        f.write(\"\\n\\n\\tYDOT  = torch.from_numpy(YDOT)\")\n",
    "        f.write(  \"\\n\\tX     = torch.from_numpy(X)\")\n",
    "        f.write(\"\\n\\tTOTAL = torch.from_numpy(TOTAL)\")\n",
    "        f.write(\"\\n\\tK     = torch.from_numpy(K)\")\n",
    "\n",
    "    ## conserved species\n",
    "    f.write('\\n\\t## Conserved species\\n')\n",
    "    for key in ODEs_cons:\n",
    "        f.write('\\n\\t## '+str(key)+'\\n')\n",
    "        for line in ODEs_cons[key]:\n",
    "            f.write('\\t'+line+'\\n')\n",
    "\n",
    "    ## non-conserved species\n",
    "    f.write('\\n\\n\\t## Non-conserved species\\n')\n",
    "    if ODEtype == 'ivp':\n",
    "        f.write('\\tHLOSS=-ACCR*Y[1]\\n')\n",
    "        f.write('\\tF = np.zeros(len(Y))\\n')\n",
    "        f.write('\\tD = np.zeros(len(Y))\\n')\n",
    "    if ODEtype == 'torch':\n",
    "        f.write('\\tHLOSS=-ACCR*Y[:,1]\\n')\n",
    "        f.write('\\tF = torch.zeros(Y.shape[1]).view(1,-1)\\n')\n",
    "        f.write('\\tD = torch.zeros(Y.shape[1]).view(1,-1)\\n')\n",
    "    for key in ODEs:\n",
    "        f.write('\\n\\t## '+str(key)+'\\n')\n",
    "        for line in ODEs[key]:\n",
    "            f.write('\\t'+line+'\\n')\n",
    "\n",
    "    ## return\n",
    "    f.write('\\n\\tYDOT=F-(D*Y)\\n')\n",
    "    f.write('\\n\\treturn YDOT')\n",
    "\n",
    "# Read the original ODE file\n",
    "with open(loc+filename+'.py', 'r') as file:\n",
    "    content = file.read()\n",
    "\n",
    "# Substitute every number between square brackets with that number minus one\n",
    "content_new = re.sub('\\[\\d+\\]', minus_one, content)\n",
    "if ODEtype == 'torch':\n",
    "    content_new = re.sub('\\[:,\\d+\\]', minus_one_torch, content_new)\n",
    "\n",
    "# Write the new ODE file\n",
    "with open(loc+filename+'.py', 'w') as file:\n",
    "    file.write(content_new)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "'''\n",
    "acodes.py: ODEs for the conserved species, manually parsed --> plug in final file\n",
    "'''\n",
    "\n",
    "\n",
    "## e-\n",
    "    T = 0.+Y[1]-Y[2]+Y[3]+Y[4]+Y[6]+Y[7]-Y[8]+Y[9]-Y[11]+Y[12]+Y[14]+Y[15]+Y[20]+Y[21]+Y[22]+Y[23]+Y[24]-Y[25]-Y[29]+Y[31]+Y[32]+Y[33]+Y[35]+Y[37]+Y[39]+Y[40]+Y[42]+Y[43]+Y[45]-Y[46]+Y[47]+Y[49]-Y[51]+Y[52]+Y[56]-Y[57]+Y[58]+Y[61]+Y[62]+Y[65]+Y[66]+Y[68]+Y[69]+Y[72]+Y[74]+Y[75]+Y[77]+Y[78]+Y[82]+Y[83]+Y[84]+Y[86]+Y[89]+Y[90]+Y[92]+Y[93]+Y[94]+Y[96]+Y[97]+Y[98]+Y[100]+Y[101]+Y[107]+Y[109]+Y[110]-Y[111]-Y[112]+Y[113]+Y[114]+Y[115]+Y[116]+Y[118]+Y[120]+Y[122]+Y[123]+Y[125]+Y[127]+Y[128]+Y[130]+Y[131]-Y[133]+Y[135]+Y[136]-Y[137]+Y[138]+Y[141]+Y[142]+Y[144]+Y[146]+Y[147]+Y[149]+Y[152]+Y[154]+Y[155]+Y[158]+Y[159]+Y[161]+Y[162]+Y[163]+Y[165]+Y[167]+Y[170]+Y[171]+Y[174]+Y[176]+Y[178]+Y[180]+Y[181]+Y[183]+Y[186]+Y[187]+Y[188]+Y[191]+Y[192]+Y[194]+Y[195]+Y[196]+Y[197]+Y[199]+Y[200]+Y[201]+Y[202]+Y[203]+Y[205]+Y[208]+Y[209]+Y[211]+Y[212]+Y[214]+Y[215]+Y[216]+Y[217]+Y[220]+Y[221]+Y[222]+Y[223]+Y[228]+Y[229]+Y[232]+Y[234]+Y[235]+Y[236]+Y[238]+Y[239]+Y[240]+Y[241]+Y[242]+Y[244]+Y[245]+Y[246]+Y[247]+Y[248]+Y[249]+Y[251]+Y[252]-Y[253]\n",
    "    T = T+Y[254]+Y[258]+Y[259]+Y[260]+Y[261]+Y[262]+Y[263]-Y[265]+Y[266]+Y[267]+Y[270]-Y[271]+Y[272]+Y[276]+Y[277]+Y[282]+Y[284]+Y[285]+Y[286]+Y[287]+Y[288]+Y[289]+Y[292]+Y[293]+Y[294]+Y[295]+Y[296]+Y[298]+Y[300]+Y[302]+Y[305]+Y[306]+Y[307]+Y[308]+Y[310]+Y[311]+Y[312]+Y[314]+Y[316]+Y[317]+Y[318]+Y[320]+Y[321]+Y[322]+Y[323]+Y[326]-Y[328]+Y[329]+Y[332]+Y[333]+Y[334]+Y[335]+Y[336]+Y[337]+Y[339]-Y[340]+Y[341]+Y[342]+Y[345]+Y[346]+Y[350]+Y[351]+Y[352]+Y[354]+Y[355]+Y[357]+Y[358]+Y[359]+Y[360]+Y[363]+Y[364]+Y[365]+Y[366]+Y[368]+Y[369]+Y[371]+Y[373]+Y[374]-Y[375]+Y[376]+Y[379]-Y[380]+Y[382]-Y[384]+Y[385]+Y[386]+Y[388]+Y[389]+Y[390]+Y[391]+Y[393]+Y[394]+Y[395]+Y[397]+Y[399]+Y[400]+Y[402]+Y[403]+Y[404]+Y[405]-Y[407]-Y[408]+Y[409]+Y[411]+Y[413]+Y[414]+Y[417]+Y[418]-Y[419]+Y[421]+Y[423]-Y[424]+Y[427]+Y[428]+Y[430]+Y[431]+Y[432]+Y[433]+Y[434]+Y[435]+Y[437]-Y[438]+Y[439]-Y[440]+Y[442]+Y[444]+Y[445]+Y[447]+Y[448]-Y[449]+Y[451]+Y[452]-Y[453]+Y[456]+Y[458]+Y[459]+Y[460]+Y[462]+Y[463]+Y[465]\n",
    "    X[0]=0.+TOTAL[0]+T\n",
    "\n",
    "    ## H2\n",
    "    T = 0.+Y[0]+Y[1]+Y[2]+2*Y[3]+3*Y[4]+Y[7]+Y[11]+Y[12]+Y[13]+2*Y[14]+2*Y[16]+3*Y[18]+Y[19]+Y[20]+3*Y[21]+2*Y[23]+4*Y[24]+2*Y[26]+4*Y[28]+Y[29]+Y[30]+Y[31]+5*Y[32]+3*Y[33]+3*Y[34]+4*Y[35]+2*Y[36]+2*Y[37]+3*Y[39]+Y[41]+Y[42]+2*Y[43]+Y[51]+Y[52]+Y[53]+2*Y[55]+2*Y[56]+Y[59]+3*Y[60]+3*Y[61]+Y[62]+Y[63]+4*Y[65]+2*Y[67]+2*Y[68]+4*Y[73]+2*Y[74]+5*Y[75]+Y[76]+Y[77]+Y[78]+3*Y[79]+Y[80]+5*Y[81]+Y[82]+Y[83]+4*Y[84]+2*Y[85]+6*Y[86]+2*Y[88]+4*Y[89]+2*Y[90]+6*Y[91]+2*Y[92]+3*Y[96]+Y[97]+3*Y[98]+Y[99]+7*Y[101]+3*Y[102]+4*Y[104]+Y[105]+2*Y[107]+4*Y[108]+4*Y[109]+Y[113]+4*Y[114]+5*Y[116]+Y[117]+2*Y[118]+Y[119]+Y[120]+2*Y[121]+Y[122]+5*Y[123]+2*Y[124]+2*Y[125]+2*Y[126]+3*Y[127]+3*Y[128]+Y[132]+Y[135]+2*Y[136]+Y[137]+Y[138]+Y[139]+2*Y[141]+2*Y[143]+2*Y[145]+Y[146]+3*Y[147]+3*Y[148]+3*Y[149]+4*Y[151]+4*Y[154]+2*Y[156]+4*Y[157]+2*Y[158]+5*Y[159]+Y[160]+Y[161]+3*Y[162]+Y[163]+3*Y[164]+6*Y[165]+4*Y[167]+2*Y[168]+6*Y[169]+2*Y[171]+2*Y[172]+2*Y[174]+2*Y[177]+Y[179]+Y[180]+Y[181]+Y[182]+3*Y[183]+Y[184]+7*Y[186]+3*Y[187]+Y[188]+Y[189]+Y[190]\n",
    "    T = T+3*Y[191]+Y[192]+3*Y[193]+Y[194]+2*Y[196]+4*Y[199]+2*Y[201]+2*Y[202]+2*Y[203]+Y[204]+4*Y[206]+2*Y[208]+4*Y[209]+2*Y[211]+Y[212]+5*Y[214]+Y[215]+Y[216]+Y[219]+2*Y[220]+Y[221]+Y[222]+2*Y[224]+2*Y[225]+3*Y[227]+2*Y[228]+6*Y[229]+6*Y[231]+2*Y[232]+2*Y[233]+3*Y[234]+2*Y[235]+6*Y[236]+6*Y[237]+Y[239]+7*Y[240]+2*Y[242]+3*Y[244]+Y[247]+4*Y[248]+3*Y[249]+7*Y[251]+3*Y[252]+Y[254]+Y[255]+3*Y[260]+Y[261]+2*Y[262]+2*Y[263]+Y[264]+Y[265]+Y[266]+2*Y[267]+2*Y[269]+Y[272]+3*Y[273]+Y[274]+3*Y[277]+Y[278]+4*Y[280]+2*Y[284]+4*Y[285]+Y[288]+5*Y[289]+Y[290]+3*Y[291]+3*Y[292]+Y[293]+Y[294]+2*Y[295]+4*Y[296]+6*Y[297]+2*Y[298]+2*Y[299]+Y[302]+5*Y[304]+3*Y[305]+7*Y[307]+3*Y[308]+Y[309]+Y[310]+6*Y[312]+Y[316]+2*Y[317]+6*Y[318]+6*Y[319]+3*Y[320]+7*Y[321]+3*Y[322]+4*Y[323]+4*Y[324]+4*Y[333]+Y[334]+Y[335]+Y[336]+Y[337]+Y[338]+5*Y[339]+Y[340]+2*Y[341]+2*Y[343]+Y[345]+3*Y[346]+4*Y[347]+4*Y[351]+2*Y[354]+Y[356]+Y[357]+Y[358]+3*Y[359]+Y[360]+3*Y[361]+Y[362]+5*Y[363]+4*Y[364]+2*Y[365]+2*Y[366]+2*Y[367]+3*Y[368]+3*Y[369]+Y[371]+Y[374]+Y[378]\n",
    "    T = T+Y[379]+Y[380]+2*Y[381]+2*Y[382]+Y[386]+Y[387]+3*Y[388]+4*Y[390]+2*Y[391]+3*Y[393]+Y[394]+5*Y[395]+6*Y[396]+6*Y[397]+7*Y[399]+Y[402]+Y[404]+Y[408]+Y[409]+Y[410]+2*Y[411]+2*Y[412]+3*Y[413]+4*Y[414]+4*Y[415]+3*Y[416]+5*Y[417]+4*Y[418]+Y[422]+Y[423]+Y[424]+2*Y[425]+2*Y[428]+Y[429]+Y[430]+3*Y[431]+4*Y[432]+2*Y[433]+5*Y[434]+3*Y[435]+Y[439]+Y[440]+Y[441]+2*Y[442]+2*Y[443]+3*Y[444]+4*Y[445]+3*Y[446]+5*Y[447]+4*Y[448]+Y[452]+Y[453]+Y[454]+2*Y[457]+2*Y[458]+Y[459]+3*Y[460]+Y[461]+2*Y[462]+3*Y[463]\n",
    "    X[1]=0.+TOTAL[1]-0.5*T"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "''' \n",
    "dcodes.py: ODEs for the conserved species, manually parsed --> plug in final file (dc fortran code) \n",
    "'''\n",
    "\n",
    "## e-\n",
    "\tT=0.-Y[1]+Y[2]+Y[3]+Y[4]+Y[5]+Y[7]+Y[8]-Y[10]-Y[12]+Y[13]+Y[15]+Y[16]+Y[18]+Y[21]+Y[22]+Y[23]+Y[24]-Y[26]-Y[30]+Y[31]+Y[32]+Y[33]+Y[35]+Y[37]+Y[39]+Y[40]+Y[42]+Y[43]+Y[44]+Y[46]+Y[47]-Y[49]+Y[52]-Y[53]+Y[55]-Y[57]+Y[58]+Y[62]+Y[63]+Y[65]+Y[66]+Y[67]+Y[71]+Y[72]+Y[73]+Y[75]+Y[79]+Y[80]+Y[81]+Y[82]+Y[84]+Y[88]+Y[89]+Y[90]+Y[91]+Y[93]+Y[95]+Y[96]+Y[97]+Y[99]+Y[101]+Y[102]+Y[103]-Y[105]+Y[106]+Y[107]+Y[109]-Y[110]+Y[112]+Y[113]+Y[117]+Y[118]+Y[119]+Y[120]+Y[121]+Y[124]+Y[125]+Y[128]+Y[130]+Y[131]-Y[133]+Y[135]+Y[137]-Y[138]+Y[139]+Y[141]+Y[142]+Y[144]+Y[146]+Y[147]+Y[149]+Y[151]+Y[153]+Y[155]+Y[157]+Y[159]+Y[161]+Y[162]+Y[164]+Y[165]+Y[166]+Y[167]+Y[169]+Y[173]+Y[176]+Y[178]+Y[180]+Y[181]+Y[183]+Y[185]+Y[186]+Y[188]+Y[190]+Y[191]+Y[192]+Y[195]+Y[196]+Y[197]+Y[198]+Y[199]+Y[200]+Y[201]+Y[202]+Y[204]+Y[205]+Y[206]+Y[207]+Y[208]+Y[215]+Y[216]+Y[217]+Y[219]+Y[220]+Y[221]+Y[222]+Y[223]+Y[224]+Y[225]+Y[227]+Y[229]+Y[230]+Y[232]+Y[237]+Y[239]+Y[240]+Y[241]+Y[244]+Y[245]+Y[246]+Y[247]+Y[248]+Y[249]\n",
    "\tT=T+Y[250]+Y[251]+Y[252]+Y[253]+Y[254]+Y[257]+Y[258]-Y[259]+Y[261]+Y[262]+Y[264]-Y[265]+Y[266]+Y[267]-Y[268]+Y[271]+Y[274]+Y[275]+Y[278]+Y[279]+Y[281]+Y[284]+Y[285]+Y[286]+Y[288]+Y[289]+Y[291]+Y[292]+Y[293]+Y[295]+Y[297]+Y[298]+Y[299]+Y[302]+Y[304]+Y[305]+Y[306]+Y[307]+Y[310]+Y[312]+Y[313]+Y[315]+Y[316]+Y[317]+Y[318]+Y[319]+Y[321]+Y[322]+Y[323]+Y[324]+Y[326]+Y[327]+Y[329]-Y[333]+Y[334]-Y[336]+Y[337]+Y[338]+Y[339]+Y[340]+Y[342]+Y[343]+Y[345]+Y[346]+Y[348]+Y[350]+Y[351]+Y[352]+Y[355]+Y[357]+Y[358]+Y[360]+Y[361]+Y[363]+Y[364]+Y[365]+Y[366]+Y[369]+Y[370]+Y[371]+Y[372]+Y[374]-Y[375]+Y[376]+Y[379]-Y[380]+Y[381]-Y[383]+Y[385]+Y[387]+Y[388]+Y[389]+Y[390]+Y[391]+Y[393]+Y[394]+Y[395]+Y[396]+Y[398]+Y[400]+Y[401]+Y[402]+Y[404]+Y[405]-Y[407]+Y[408]-Y[409]+Y[412]+Y[413]+Y[415]+Y[417]+Y[418]+Y[420]-Y[421]+Y[422]-Y[423]+Y[425]+Y[428]+Y[429]+Y[430]+Y[432]+Y[433]+Y[434]+Y[435]+Y[436]-Y[437]-Y[439]+Y[440]+Y[443]+Y[444]+Y[445]+Y[446]+Y[448]+Y[449]-Y[450]-Y[452]+Y[454]+Y[457]+Y[458]+Y[459]+Y[460]+Y[462]+Y[463]+Y[465]\n",
    "\tX[0]=0.+TOTAL[0]+T\n",
    "\n",
    "\t## H2\n",
    "\tT=0.+Y[0]+Y[1]+Y[2]+2*Y[3]+3*Y[4]+Y[7]+Y[11]+Y[12]+Y[13]+2*Y[14]+2*Y[15]+Y[18]+Y[19]+3*Y[20]+3*Y[21]+2*Y[22]+4*Y[24]+2*Y[25]+4*Y[27]+3*Y[29]+Y[30]+Y[31]+3*Y[32]+5*Y[33]+Y[34]+4*Y[35]+2*Y[36]+2*Y[37]+3*Y[39]+Y[41]+Y[42]+2*Y[43]+Y[51]+Y[52]+Y[53]+2*Y[54]+2*Y[55]+Y[59]+Y[60]+3*Y[61]+3*Y[62]+Y[63]+2*Y[64]+4*Y[70]+2*Y[71]+2*Y[72]+4*Y[73]+Y[75]+5*Y[76]+Y[77]+Y[78]+Y[79]+Y[80]+5*Y[81]+Y[82]+3*Y[83]+2*Y[84]+2*Y[85]+6*Y[86]+2*Y[88]+4*Y[89]+4*Y[91]+2*Y[92]+6*Y[93]+3*Y[94]+3*Y[95]+3*Y[96]+7*Y[97]+Y[100]+Y[101]+4*Y[104]+4*Y[107]+4*Y[109]+Y[111]+2*Y[112]+Y[113]+4*Y[114]+Y[116]+2*Y[117]+Y[118]+5*Y[119]+5*Y[120]+Y[121]+2*Y[122]+Y[123]+2*Y[124]+3*Y[125]+2*Y[126]+2*Y[127]+3*Y[128]+Y[134]+Y[135]+Y[136]+2*Y[137]+Y[138]+Y[139]+2*Y[140]+2*Y[141]+2*Y[143]+3*Y[146]+Y[147]+3*Y[148]+3*Y[149]+4*Y[150]+2*Y[152]+2*Y[153]+4*Y[155]+4*Y[156]+5*Y[159]+3*Y[160]+3*Y[161]+Y[162]+Y[163]+Y[164]+2*Y[165]+2*Y[167]+2*Y[168]+6*Y[169]+6*Y[170]+2*Y[172]+2*Y[174]+4*Y[176]+7*Y[178]+3*Y[179]+Y[181]+Y[182]+Y[183]+Y[184]+3*Y[185]+3*Y[186]\n",
    "\tT=T+Y[187]+Y[188]+Y[190]+Y[191]+3*Y[192]+Y[193]+Y[194]+2*Y[197]+2*Y[199]+2*Y[200]+4*Y[202]+Y[204]+2*Y[205]+4*Y[206]+2*Y[207]+2*Y[208]+Y[210]+4*Y[211]+Y[214]+Y[215]+Y[216]+Y[217]+5*Y[220]+Y[221]+2*Y[222]+Y[223]+2*Y[224]+6*Y[225]+6*Y[226]+3*Y[227]+6*Y[228]+2*Y[230]+3*Y[231]+2*Y[233]+2*Y[234]+2*Y[237]+2*Y[238]+6*Y[239]+Y[240]+3*Y[241]+2*Y[245]+4*Y[246]+3*Y[249]+7*Y[250]+7*Y[251]+3*Y[252]+Y[253]+Y[256]+3*Y[257]+Y[261]+2*Y[262]+Y[263]+2*Y[264]+Y[265]+Y[266]+2*Y[270]+2*Y[271]+3*Y[272]+Y[274]+Y[276]+Y[277]+3*Y[278]+4*Y[280]+2*Y[281]+4*Y[284]+5*Y[288]+Y[289]+3*Y[290]+Y[291]+3*Y[292]+Y[293]+Y[294]+6*Y[296]+2*Y[297]+4*Y[298]+2*Y[299]+2*Y[301]+5*Y[303]+3*Y[304]+3*Y[305]+7*Y[306]+Y[307]+6*Y[312]+Y[314]+Y[315]+Y[316]+2*Y[317]+6*Y[318]+3*Y[319]+6*Y[320]+3*Y[321]+4*Y[322]+7*Y[323]+4*Y[324]+4*Y[328]+Y[334]+Y[335]+Y[336]+Y[337]+Y[338]+Y[339]+5*Y[340]+2*Y[342]+2*Y[344]+3*Y[345]+Y[346]+2*Y[348]+4*Y[351]+4*Y[353]+Y[356]+5*Y[357]+Y[358]+3*Y[359]+Y[360]+3*Y[361]+Y[362]+Y[363]+2*Y[364]+4*Y[365]+2*Y[366]+2*Y[367]\n",
    "\tT=T+3*Y[369]+3*Y[370]+Y[371]+Y[374]+Y[378]+Y[379]+Y[380]+2*Y[382]+2*Y[385]+Y[386]+Y[387]+3*Y[388]+2*Y[389]+4*Y[390]+Y[393]+5*Y[394]+3*Y[395]+6*Y[396]+6*Y[397]+7*Y[398]+Y[401]+Y[404]+Y[408]+Y[409]+Y[410]+2*Y[411]+2*Y[412]+3*Y[413]+4*Y[414]+4*Y[415]+3*Y[416]+5*Y[417]+4*Y[418]+Y[422]+Y[423]+Y[424]+2*Y[427]+2*Y[428]+Y[429]+3*Y[430]+Y[431]+4*Y[432]+2*Y[433]+3*Y[434]+5*Y[435]+Y[439]+Y[440]+Y[441]+2*Y[442]+2*Y[443]+3*Y[444]+4*Y[445]+5*Y[446]+3*Y[447]+4*Y[448]+Y[452]+Y[453]+Y[454]+2*Y[456]+2*Y[457]+Y[459]+3*Y[460]+Y[461]+2*Y[462]+3*Y[463]\n",
    "\tX[1]=0.+TOTAL[1]-0.5*T"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Translate to torch\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "## testing\n",
    "\n",
    "str1 = 'F[3]=0.+K[456]*X[1]*Y[40]*HNR+K[457]*X[1]*Y[6]*HNR+K[730]*X[1]+K[2564]*Y[1]*Y[76]*HNR+K[3064]*Y[0]*Y[7]*HNR+K[5921]*Y[4]+K[6092]*Y[1]*Y[0]*HNR+K[6396]*Y[4]+K[6732]*Y[4]'\n",
    "str2 = 'D[3]=0.+K[435]*Y[50]*HNR+K[436]*Y[55]*HNR+K[437]*Y[73]*HNR+K[438]*Y[53]*HNR+K[439]*Y[16]*HNR+K[440]*Y[91]*HNR+K[441]*Y[28]*HNR+K[442]*Y[13]*HNR+K[443]*Y[54]*HNR+K[444]*Y[71]*HNR+K[445]*Y[88]*HNR+K[446]*Y[36]*HNR+K[447]*Y[126]*HNR+K[448]*Y[59]*HNR+K[449]*Y[76]*HNR+K[450]*Y[26]*HNR+K[451]*Y[34]*HNR+K[452]*Y[19]*HNR+K[453]*Y[87]*HNR+K[454]*Y[103]*HNR+K[455]*Y[30]*HNR+K[488]*Y[0]*HNR+K[1238]*X[0]*HNR+K[2588]*Y[50]*HNR+K[2589]*Y[55]*HNR+K[2590]*Y[73]*HNR+K[2591]*Y[73]*HNR+K[2592]*Y[53]*HNR+K[2593]*Y[10]*HNR+K[2594]*Y[16]*HNR+K[2595]*Y[91]*HNR+K[2596]*Y[91]*HNR+K[2597]*Y[28]*HNR+K[2598]*Y[28]*HNR+K[2599]*Y[13]*HNR+K[2600]*Y[54]*HNR+K[2601]*Y[213]*HNR+K[2602]*Y[71]*HNR+K[2603]*X[1]*HNR+K[2604]*Y[88]*HNR+K[2605]*Y[36]*HNR+K[2606]*Y[126]*HNR+K[2607]*Y[126]*HNR+K[2608]*Y[76]*HNR+K[2609]*Y[5]*HNR+K[2610]*Y[64]*HNR+K[2611]*Y[17]*HNR+K[2612]*Y[19]*HNR+K[2613]*Y[87]*HNR+K[2614]*Y[103]*HNR+K[2615]*Y[27]*HNR+K[2616]*Y[30]*HNR+K[4921]*Y[2]*HNR+K[5902]+K[6377]+K[6713]'\n",
    "\t\n",
    "str1 = list(str1)\n",
    "\n",
    "for i in range(len(str1)):\n",
    "    if str1[i] == 'Y':\n",
    "        str1[i+1] = '[:,'\n",
    "    \n",
    "str1 = \"\".join(str1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'F[3]=0.+K[456]*X[1]*Y[:,40]*HNR+K[457]*X[1]*Y[:,6]*HNR+K[730]*X[1]+K[2564]*Y[:,1]*Y[:,76]*HNR+K[3064]*Y[:,0]*Y[:,7]*HNR+K[5921]*Y[:,4]+K[6092]*Y[:,1]*Y[:,0]*HNR+K[6396]*Y[:,4]+K[6732]*Y[:,4]'"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "str1"
   ]
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "4cfec3ae79f839ed35c55406a3ad8c3c1082e0129819951cbf0a51bf47981102"
  },
  "kernelspec": {
   "display_name": "Python 3.9.16 ('mace')",
   "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.13"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}