c5-simulation-driver.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Simulation command generation from simulation sheet\n",
    ">'It's turtles all the way down!'\n",
    "\n",
    "Plan: \n",
    "- Have array of sample settings, maybe optimized to give comparable $\\sigma t$ to give comparable signal strength."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from definitions import *\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import os \n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Rs = [2e-6, 300e-9, 50e-9]\n",
    "\n",
    "thickness_map = {\n",
    "    (3.112, 2e-6):  3e-3,\n",
    "    (3.112, 300e-9): 10e-3,\n",
    "    (3.112, 50e-9): 20e-3,\n",
    "    (9.9177, 2e-6):  0.5e-3,\n",
    "    (9.9177, 300e-9): 3e-3,\n",
    "    (9.9177, 50e-9): 10e-3,\n",
    "    (4.321, 2e-6):  1e-3,\n",
    "    (8, 2e-6):  1e-3,\n",
    "    (4.321, 300e-9): 10e-3,\n",
    "    (8, 300e-9): 5e-3,\n",
    "    (4.321, 50e-9): 10e-3,\n",
    "    (8, 50e-9): 10e-3,\n",
    "}\n",
    "for (L0, R) in thickness_map:\n",
    "    t = thickness_map[(L0,R)]\n",
    "    # print(L0,R,t)\n",
    "    print(F\"R = {R * 1e9}nm, t = {round(t * 1e3,2)}mm, lambda = {L0} Å:\")\n",
    "    st = s_t(R,t, L0 * 1e-10, phi, delta_rho)\n",
    "    permissible = st < 0.8 and st >= 0.1\n",
    "    print(f\"\\ts*t: {round(st, 4)}\\n\\t0.1 <= s*t <= 0.8: {permissible}\\n\\texp(-s*t): {round(np.exp(-st), 4)}\")\n",
    "    sigma = s_t(R,t, L0 * 1e-10) / t\n",
    "    t_optimum = 0.38 / sigma"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from util import *\n",
    "from instrument import Instrument\n",
    "\n",
    "instrs = load_instruments('instruments.csv')\n",
    "\n",
    "N_foil = 10000000\n",
    "N_iso_wsp = 100000000\n",
    "\n",
    "N_steps = 5\n",
    "\n",
    "shift_Ls = False\n",
    "\n",
    "with open('simulate.sh', 'w') as f:\n",
    "    f.write('#!/bin/bash\\n# Simulation script automatically generated by simulation-driver.ipynb, use this to create variants of it\\n')\n",
    "    f.write('rm -rf data\\n')\n",
    "    for instr in instrs[3:-6]:\n",
    "        prec_type = instr.prec_type\n",
    "        if instr.prec_type == 'wsp':\n",
    "            prec_type = 'iwsp'\n",
    "\n",
    "        # if instr.prec_type == 'foil':\n",
    "        #     continue\n",
    "        # By = 0.01\n",
    "        \n",
    "        d_min, d_max = instr.delta_range()\n",
    "        if shift_Ls:\n",
    "            L_s_scaling = 20e-9 / d_min\n",
    "            L_s_new = instr.L_s * L_s_scaling\n",
    "\n",
    "            # print(L_s_new)\n",
    "            theta_a_max = 15e-3\n",
    "            L_min = round(detector_size / (2 * np.tan(theta_a_max)),4)\n",
    "            L_s_new_corrected = max(L_s_new, L_min)\n",
    "            print(L_s_new_corrected, L_min, L_s_new)\n",
    "            instr.L_s = L_s_new_corrected\n",
    "        d_min, d_max = instr.delta_range()\n",
    "        print(\"-------------------------\\n\"+str(instr))\n",
    "        for R in Rs:\n",
    "            t = thickness_map[(instr.L0 * 1e10, R)]\n",
    "            mode = 'GPU'\n",
    "            if instr.prec_type == 'foil':\n",
    "                mode = 'CPU'\n",
    "            # print(f\"delta_y range: {round(d_min * 1e9,2)} - {round(d_max * 1e9, 2)}nm\")\n",
    "            # print(f\"\\tRange of interest for R = {R * 1e9}nm: 0 - {3 * R * 1e9}\")\n",
    "            print(f\"\\tR = {R * 1e9}nm:\")\n",
    "            print(f\"\\t\\tRange of interest:  {round(0.1 * R * 1e9)} - {round(3 * R * 1e9)}nm\")\n",
    "\n",
    "            \n",
    "\n",
    "            overlap = find_overlap((d_min, d_max), (0.0, 3 * R))\n",
    "            if overlap == None:\n",
    "                print(\"No overlap!\")\n",
    "                f.write(f'echo \"Skipping {instr.id}_{int(R * 1e10)} due to non-overlapping ranges!\\n')\n",
    "\n",
    "                raise Exception(\"No overlap!\")\n",
    "            else:\n",
    "                (a,b) = overlap\n",
    "                fraction = (b - a) / (3 * R) * 100\n",
    "                log_fraction = log_overlap_percentage((d_min, d_max), (0.1 * R, 3 * R))\n",
    "                print(f\"\\t\\tOverlap: {round(a * 1e9,2)} - {round(b * 1e9, 2)}nm ({round(fraction,1)}% of linear range, {round(log_fraction,1)}% of log range)\")\n",
    "                d_min_B_field, d_max_B_field = instr.delta_range_B_field()\n",
    "                max_ratio = b / d_max_B_field\n",
    "                min_ratio = a / d_min_B_field\n",
    "                By_min = instr.By_min * min_ratio\n",
    "                # The simulation parameter By range corresponds to the setting of B1\n",
    "                # The maximal setting is when B1 = Bmax * L2 / L1 so that\n",
    "                # B2 = B1 * L1/L2 = Bmax, having the greater field\n",
    "                # Take the ratio from this depending on the overlap\n",
    "                By_max = instr.By_max * instr.L_2 / instr.L_1 * max_ratio\n",
    "                print(f\"\\t\\t=>Simulating B field range {By_min} - {By_max}\")\n",
    "                # print(By_max)\n",
    "                if instr.prec_type == 'foil':\n",
    "                    N = N_foil\n",
    "                else:\n",
    "                    N = N_iso_wsp\n",
    "                sim_descr = f\"# Name: {instr.name}, R = {round(R * 1e9)} nm\\n# Simulating delta range {round(a * 1e9, 4)} - {round(b * 1e9, 4)} nm using B field range {round(By_min * 1e3,4)} - {round(By_max * 1e3,4)} mT (limits: {round(instr.By_min * 1e3,4)} - {round(instr.By_max * 1e3,4)} mT)\"\n",
    "\n",
    "                # print(sim_cmd)\n",
    "                f.write(f'{sim_descr}\\n')\n",
    "                sim_cmd = f\"./full-simulation.sh {N} 4 {N_steps} {By_min},{By_max} {instr.L0 * 1e10} {instr.DL* 1e10} {R * 1e10} {t} {instr.theta_0} {instr.L_s} {prec_type} {prec_type}_empty {mode}; rm -rf data_{instr.id}_{int(R * 1e10)}; mv data data_{instr.id}_{int(R * 1e10)}\"\n",
    "                # print(sim_cmd)\n",
    "                f.write(f'{sim_cmd}\\n')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Modulated plot for cover"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from util import *\n",
    "from instrument import Instrument\n",
    "\n",
    "instrs = load_instruments('instruments_2.csv')\n",
    "\n",
    "N_foil = 10000000\n",
    "N_iso_wsp = 100000000\n",
    "\n",
    "N_steps = 30\n",
    "\n",
    "shift_Ls = False\n",
    "R = 3000e-10\n",
    "\n",
    "with open('simulate.sh', 'w') as f:\n",
    "    f.write('#!/bin/bash\\n# Simulation script automatically generated by simulation-driver.ipynb, use this to create variants of it\\n')\n",
    "    f.write('rm -rf data\\n')\n",
    "    for instr in instrs:\n",
    "        prec_type = instr.prec_type\n",
    "        if instr.prec_type == 'wsp':\n",
    "            prec_type = 'iwsp'\n",
    "\n",
    "        # if instr.prec_type == 'foil':\n",
    "        #     continue\n",
    "        # By = 0.01\n",
    "        \n",
    "        d_min, d_max = instr.delta_range()\n",
    "        if shift_Ls:\n",
    "            L_s_scaling = 20e-9 / d_min\n",
    "            L_s_new = instr.L_s * L_s_scaling\n",
    "\n",
    "            # print(L_s_new)\n",
    "            theta_a_max = 15e-3\n",
    "            L_min = round(detector_size / (2 * np.tan(theta_a_max)),4)\n",
    "            L_s_new_corrected = max(L_s_new, L_min)\n",
    "            print(L_s_new_corrected, L_min, L_s_new)\n",
    "            instr.L_s = L_s_new_corrected\n",
    "        d_min, d_max = instr.delta_range()\n",
    "        print(\"-------------------------\\n\"+str(instr))\n",
    "        \n",
    "        t = 0.001\n",
    "        mode = 'GPU'\n",
    "        if instr.prec_type == 'foil':\n",
    "            mode = 'CPU'\n",
    "        # print(f\"delta_y range: {round(d_min * 1e9,2)} - {round(d_max * 1e9, 2)}nm\")\n",
    "        # print(f\"\\tRange of interest for R = {R * 1e9}nm: 0 - {3 * R * 1e9}\")\n",
    "        print(f\"\\tR = {R * 1e9}nm:\")\n",
    "        print(f\"\\t\\tRange of interest:  {round(0.1 * R * 1e9)} - {round(3 * R * 1e9)}nm\")\n",
    "\n",
    "        \n",
    "\n",
    "        overlap = find_overlap((d_min, d_max), (0.0, 3 * R))\n",
    "        if overlap == None:\n",
    "            print(\"No overlap!\")\n",
    "            f.write(f'echo \"Skipping {instr.id}_{int(R * 1e10)} due to non-overlapping ranges!\\n')\n",
    "\n",
    "            raise Exception(\"No overlap!\")\n",
    "        else:\n",
    "            (a,b) = overlap\n",
    "            fraction = (b - a) / (3 * R) * 100\n",
    "            log_fraction = log_overlap_percentage((d_min, d_max), (0.1 * R, 3 * R))\n",
    "            print(f\"\\t\\tOverlap: {round(a * 1e9,2)} - {round(b * 1e9, 2)}nm ({round(fraction,1)}% of linear range, {round(log_fraction,1)}% of log range)\")\n",
    "            d_min_B_field, d_max_B_field = instr.delta_range_B_field()\n",
    "            max_ratio = b / d_max_B_field\n",
    "            min_ratio = a / d_min_B_field\n",
    "            By_min = instr.By_min * min_ratio\n",
    "            # The simulation parameter By range corresponds to the setting of B1\n",
    "            # The maximal setting is when B1 = Bmax * L2 / L1 so that\n",
    "            # B2 = B1 * L1/L2 = Bmax, having the greater field\n",
    "            # Take the ratio from this depending on the overlap\n",
    "            By_max = instr.By_max * instr.L_2 / instr.L_1 * max_ratio\n",
    "            print(f\"\\t\\t=>Simulating B field range {By_min} - {By_max}\")\n",
    "            # print(By_max)\n",
    "            if instr.prec_type == 'foil':\n",
    "                N = N_foil\n",
    "            else:\n",
    "                N = N_iso_wsp\n",
    "            sim_descr = f\"# Name: {instr.name}, R = {round(R * 1e9)} nm\\n# Simulating delta range {round(a * 1e9, 4)} - {round(b * 1e9, 4)} nm using B field range {round(By_min * 1e3,4)} - {round(By_max * 1e3,4)} mT (limits: {round(instr.By_min * 1e3,4)} - {round(instr.By_max * 1e3,4)} mT)\"\n",
    "\n",
    "            # print(sim_cmd)\n",
    "            f.write(f'{sim_descr}\\n')\n",
    "            sim_cmd = f\"./full-simulation.sh {N} 4 {N_steps} {By_min},{By_max} {instr.L0 * 1e10} {instr.DL* 1e10} {R * 1e10} {t} {instr.theta_0} {instr.L_s} {prec_type} {prec_type}_empty {mode}; rm -rf data_{instr.id}_{int(R * 1e10)}; mv data data_{instr.id}_{int(R * 1e10)}\"\n",
    "            # print(sim_cmd)\n",
    "            f.write(f'{sim_cmd}\\n')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Optimized instrument simulation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from util import *\n",
    "from instrument import Instrument\n",
    "import pickle"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with open('objects/optimized_instruments_2.pkl', 'rb') as file: \n",
    "    instrs = pickle.load(file) \n",
    "\n",
    "N_foil = 10000000\n",
    "N_iso_wsp = 1000000\n",
    "# N_steps = 10\n",
    "# N_steps = 30\n",
    "N_steps = 50\n",
    "\n",
    "shift_Ls = False\n",
    "\n",
    "with open('simulate.sh', 'w') as f:\n",
    "    f.write('#!/bin/bash\\n# Simulation script automatically generated by simulation-driver.ipynb, use this to create variants of it\\n')\n",
    "    f.write('rm -rf data\\n')\n",
    "    for (instr,_) in instrs:\n",
    "        prec_type = instr.prec_type\n",
    "        if instr.prec_type == 'wsp':\n",
    "            prec_type = 'iwsp'\n",
    "\n",
    "        if instr.prec_type != 'foil2':\n",
    "            continue\n",
    "        instr.id = instr.name\n",
    "        instr.id=instr.id.replace(' ','_')\n",
    "\n",
    "        instr.prec_type='foil'\n",
    "\n",
    "        # By = 0.01\n",
    "        \n",
    "        d_min, d_max = instr.delta_range()\n",
    "        if shift_Ls:\n",
    "            L_s_scaling = 20e-9 / d_min\n",
    "            L_s_new = instr.L_s * L_s_scaling\n",
    "\n",
    "            # print(L_s_new)\n",
    "            theta_a_max = 15e-3\n",
    "            L_min = round(detector_size / (2 * np.tan(theta_a_max)),4)\n",
    "            L_s_new_corrected = max(L_s_new, L_min)\n",
    "            print(L_s_new_corrected, L_min, L_s_new)\n",
    "            instr.L_s = L_s_new_corrected\n",
    "        d_min, d_max = instr.delta_range()\n",
    "        print(\"-------------------------\\n\"+str(instr))\n",
    "        for R in Rs:\n",
    "            t = thickness_map[(instr.L0 * 1e10, R)]\n",
    "            mode = 'GPU'\n",
    "            # if instr.prec_type == 'foil':\n",
    "            #     mode = 'CPU'\n",
    "            # print(f\"delta_y range: {round(d_min * 1e9,2)} - {round(d_max * 1e9, 2)}nm\")\n",
    "            # print(f\"\\tRange of interest for R = {R * 1e9}nm: 0 - {3 * R * 1e9}\")\n",
    "            print(f\"\\tR = {R * 1e9}nm:\")\n",
    "            print(f\"\\t\\tRange of interest:  {round(0.1 * R * 1e9)} - {round(3 * R * 1e9)}nm\")\n",
    "\n",
    "            \n",
    "\n",
    "            overlap = find_overlap((d_min, d_max), (0.0, 3 * R))\n",
    "            if overlap == None:\n",
    "                print(\"No overlap!\")\n",
    "                f.write(f'echo \"Skipping {instr.id}_{int(R * 1e10)} due to non-overlapping ranges!\\n')\n",
    "\n",
    "                raise Exception(\"No overlap!\")\n",
    "            else:\n",
    "                (a,b) = overlap\n",
    "                fraction = (b - a) / (3 * R) * 100\n",
    "                log_fraction = log_overlap_percentage((d_min, d_max), (0.1 * R, 3 * R))\n",
    "                print(f\"\\t\\tOverlap: {round(a * 1e9,2)} - {round(b * 1e9, 2)}nm ({round(fraction,1)}% of linear range, {round(log_fraction,1)}% of log range)\")\n",
    "                d_min_B_field, d_max_B_field = instr.delta_range_B_field()\n",
    "                max_ratio = b / d_max_B_field\n",
    "                min_ratio = a / d_min_B_field\n",
    "                By_min = instr.By_min * min_ratio\n",
    "                # The simulation parameter By range corresponds to the setting of B1\n",
    "                # The maximal setting is when B1 = Bmax * L2 / L1 so that\n",
    "                # B2 = B1 * L1/L2 = Bmax, having the greater field\n",
    "                # Take the ratio from this depending on the overlap\n",
    "                By_max = instr.By_max * instr.L_2 / instr.L_1 * max_ratio\n",
    "                print(f\"\\t\\t=>Simulating B field range {By_min} - {By_max}\")\n",
    "                # print(By_max)\n",
    "                if instr.prec_type == 'foil':\n",
    "                    N = N_foil\n",
    "                else:\n",
    "                    N = N_iso_wsp\n",
    "                sim_descr = f\"# Name: {instr.name}, R = {round(R * 1e9)} nm\\n# Simulating delta range {round(a * 1e9, 4)} - {round(b * 1e9, 4)} nm using B field range {round(By_min * 1e3,4)} - {round(By_max * 1e3,4)} mT (limits: {round(instr.By_min * 1e3,4)} - {round(instr.By_max * 1e3,4)} mT)\"\n",
    "\n",
    "                # print(sim_cmd)\n",
    "                f.write(f'{sim_descr}\\n')\n",
    "                sim_cmd = f\"./full-simulation-optimized.sh {N} 4 {N_steps} {By_min},{By_max} {instr.L0 * 1e10} {instr.DL* 1e10} {R * 1e10} {t} {instr.theta_0} {instr.L_1} {instr.L_2} {instr.L_s} {prec_type} {prec_type}_empty {mode}; rm -rf data_{instr.id}_{int(R * 1e10)}; mv data data_{instr.id}_{int(R * 1e10)}\"\n",
    "                print(sim_cmd)\n",
    "                f.write(f'{sim_cmd}\\n')\n",
    "        break"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "L_1 = instr.L_1\n",
    "L_2 = instr.L_2\n",
    "length_boost = 4.75 / L_1\n",
    "length_boost\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Solve the equation\n",
    "$$R_{pixel} = \\frac{\\sin(\\pi p f)}{\\pi p f} = 0.75$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sympy import sin, nsolve, Symbol, pi\n",
    "x = Symbol('x')\n",
    "\n",
    "V_reduction = 0.7\n",
    "solution = float(nsolve(sin(x) / x  - V_reduction, 0.4))\n",
    "x_range = np.linspace(-4 * np.pi, 4 * np.pi, 10000)\n",
    "y = np.sin(x_range) / x_range\n",
    "plt.plot(x_range, y)\n",
    "plt.axvline(solution, linestyle='--')\n",
    "plt.axhline(V_reduction, linestyle='--')\n",
    "plt.grid()\n",
    "# solution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f_min = solution / (np.pi * detector_pixel_size)\n",
    "f_min, f_min * 2e-10 * 1.6"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "bep",
   "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.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}