spectral_sim_ec_tinting.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bc5a0f03-1b1b-4bf8-b9e7-232ef0fb489f",
   "metadata": {},
   "outputs": [],
   "source": [
    "import datetime\n",
    "import frads as fr\n",
    "import matplotlib as mpl\n",
    "import numpy as np\n",
    "import os\n",
    "import pyradiance as pr\n",
    "import pywincalc\n",
    "import subprocess as sp\n",
    "from matplotlib.image import imread\n",
    "import matplotlib.pyplot as plt\n",
    "from typing import Tuple"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b22afce3",
   "metadata": {},
   "outputs": [],
   "source": [
    "total_frames = 180 #in 3fps for 1 min animation\n",
    "spectrum = np.arange(390, 780, 20)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3de9209f-4324-4d7d-9c52-69442ab19de0",
   "metadata": {},
   "outputs": [],
   "source": [
    "room = pr.Scene(\n",
    "    sid=\"room\",\n",
    "    surfaces=[\"./Groups/exterior.rad\", \"./Groups/room.rad\"],\n",
    "    materials=[\"./Materials/material.rad\", \"./Materials/ec60.rad\"]\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ae020d9b",
   "metadata": {},
   "outputs": [],
   "source": [
    "room_ec01 = pr.Scene(\n",
    "    sid=\"room\",\n",
    "    surfaces=[\"./Groups/exterior.rad\", \"./Groups/room.rad\"],\n",
    "    materials=[\"./Materials/material.rad\", \"./Materials/ec01.rad\"]\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6ec7f4c9",
   "metadata": {},
   "outputs": [],
   "source": [
    "# room with electric light\n",
    "# no sun\n",
    "room_elgt = pr.Scene(\n",
    "    sid=\"room_elgt\",\n",
    "    surfaces=[\"./Groups/exterior.rad\", \"./Groups/room.rad\", \"./Objects/elec_light.rad\"],\n",
    "    materials=[\"./Materials/material.rad\", \"./Materials/ec60.rad\"]\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0e868583",
   "metadata": {},
   "outputs": [],
   "source": [
    "ssky = pr.genssky(\n",
    "    dt=datetime.datetime(2024,12,21,12),\n",
    "    latitude=37,\n",
    "    longitude=122,\n",
    "    timezone=120,\n",
    "    year=2024,\n",
    ")\n",
    "for prim in pr.parse_primitive(ssky.decode()):\n",
    "    room.add_source(prim)\n",
    "\n",
    "room.add_source(pr.Primitive(\"skyfunc\", \"glow\", \"skyglow\", [], [1,1,1,0]))\n",
    "room.add_source(pr.Primitive(\"skyglow\", \"source\", \"skydome\", [], [0,0,1,180]))\n",
    "room.add_source(pr.Primitive(\"skyglow\", \"source\", \"groundplane\", [], [0,0,-1,180]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "544ebe3c",
   "metadata": {},
   "outputs": [],
   "source": [
    "for prim in pr.parse_primitive(ssky.decode()):\n",
    "    room_ec01.add_source(prim)\n",
    "\n",
    "room_ec01.add_source(pr.Primitive(\"skyfunc\", \"glow\", \"skyglow\", [], [1,1,1,0]))\n",
    "room_ec01.add_source(pr.Primitive(\"skyglow\", \"source\", \"skydome\", [], [0,0,1,180]))\n",
    "room_ec01.add_source(pr.Primitive(\"skyglow\", \"source\", \"groundplane\", [], [0,0,-1,180]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0e467c30",
   "metadata": {},
   "outputs": [],
   "source": [
    "view = pr.View(\n",
    "        vtype='a',\n",
    "        position=(70, 60, 49),\n",
    "        direction=(1, -1, 0),\n",
    "        horiz=180,\n",
    "        vert=180\n",
    "    )\n",
    "view_pt = (' '.join(str(item) for item in view.position + view.direction)).encode()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1d5700d3",
   "metadata": {},
   "outputs": [],
   "source": [
    "floor_primitives = fr.unpack_primitives(\"./Objects/floor.rad\")\n",
    "floor_prim = fr.parse_polygon(floor_primitives[0])\n",
    "grid = fr.gen_grid(floor_prim, 28.875, 12) #inches\n",
    "grid = ' '.join(' '.join(map(str, g)) for g in grid).encode()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "302195e7",
   "metadata": {},
   "outputs": [],
   "source": [
    "tr_params = pr.SamplingParameters(\n",
    "    ab = 2,\n",
    "    aa = 0,\n",
    "    ad = 2048,\n",
    "    lw = 1e-5,\n",
    "    dj = 0.8,\n",
    "    cs = 20,\n",
    "    co = True\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "132782c9",
   "metadata": {},
   "outputs": [],
   "source": [
    "render_params = pr.SamplingParameters(\n",
    "    ab = 2,\n",
    "    aa = 0,\n",
    "    ad = 16,\n",
    "    lw = 1e-5, \n",
    "    dj = 0.8,\n",
    "    cs = 20,\n",
    "    co = True\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "be4da0f1",
   "metadata": {},
   "outputs": [],
   "source": [
    "def interpolate_factor(range: Tuple[float, float], current_frame: int, frame_total: int) -> float:\n",
    "    '''\n",
    "    Calculate an interpolated factor using exponential interpolation.\n",
    "\n",
    "    Parameters:\n",
    "    range (tuple): A tuple of two numbers (start, end) representing the range to interpolate between.\n",
    "    current_frame (int): The current frame number for which to calculate the interpolated value.\n",
    "    frame_total (int): The total number of frames in the interpolation sequence. \n",
    "\n",
    "    Returns:\n",
    "    float: The interpolated factor for the current frame, normalized to a 0-1 range.\n",
    "    '''\n",
    "    x1, x2 = 1, frame_total\n",
    "    y1, y2 = range\n",
    "    b = (y2 / y1) ** (1 / (x2 - x1))\n",
    "    a = y1 / (b ** x1)\n",
    "    value = a * b ** current_frame\n",
    "    return value / y1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "78cb2b92",
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_ave_wpi(octree, sensors, params) -> float:\n",
    "    '''\n",
    "    Get average workplane illuminance\n",
    "    '''\n",
    "    wpi_dimming=pr.rtrace(sensors, octree, params=params, irradiance=True, header=False).decode()\n",
    "    wpi_dimming = [float(i) for i in wpi_dimming.split()]\n",
    "    ave_wpi_dimming = sum(wpi_dimming)/len(wpi_dimming)\n",
    "    return ave_wpi_dimming"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f77ab643",
   "metadata": {},
   "outputs": [],
   "source": [
    "spectral_transmittance_60= np.loadtxt(\"Materials/sageglass-sr20_7mm-lami-full-clear-60t.csv\",skiprows=38, delimiter=\",\")\n",
    "spectral_transmittance_60 = spectral_transmittance_60[np.isin(spectral_transmittance_60[:,0], spectrum/1000)][:,1]\n",
    "\n",
    "spectral_transmittance_01 = np.loadtxt(\"Materials/sageglass-sr20_7mm-lami-fully-tinted-1t.csv\",skiprows=38, delimiter=\",\")\n",
    "spectral_transmittance_01 = spectral_transmittance_01[np.isin(spectral_transmittance_01[:,0], spectrum/1000)][:,1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e4ce8a92",
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_tinted_factors(current_frame:int) -> np.array:\n",
    "    '''\n",
    "    Get tinted factors for the current frame \n",
    "\n",
    "    Parameters:\n",
    "    current_frame (int): The current frame number for which to calculate the tinted factor.\n",
    "\n",
    "    Returns:\n",
    "    float: The tinted factor for the given frame number.\n",
    "    '''\n",
    "    tinted_factor = []\n",
    "    for st60, st01 in zip(spectral_transmittance_60, spectral_transmittance_01):\n",
    "        tinted_factor.append(interpolate_factor((st60, st01), current_frame, total_frames))\n",
    "    return np.array(tinted_factor) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bd744356",
   "metadata": {},
   "outputs": [],
   "source": [
    "ird_tinted = pr.rtrace(view_pt, room.octree, params=tr_params.args(), irradiance=True, header=False).decode().split()\n",
    "ird_dimming = pr.rtrace(view_pt, room_elgt.octree, params=tr_params.args(), irradiance=True, header=False).decode().split()\n",
    "ird_tinted = list(map(float, ird_tinted)) # convert str to float\n",
    "ird_dimming = list(map(float, ird_dimming)) # convert str to float\n",
    "mlux_tinted = float(pr.rtrace(view_pt, room.octree, params=tr_params.args()+[\"-pM\"], irradiance=True, header=False).decode().split()[0])\n",
    "mlux_dimming = float(pr.rtrace(view_pt, room_elgt.octree, params=tr_params.args()+[\"-pM\"], irradiance=True, header=False).decode().split()[0])\n",
    "mlux_ec01 = float(pr.rtrace(view_pt, room_ec01.octree, params=tr_params.args()+[\"-pM\"], irradiance=True, header=False).decode().split()[0])\n",
    "ev_tinted = float(pr.rtrace(view_pt, room.octree, params=tr_params.args()+[\"-pY\"], irradiance=True, header=False).decode().split()[0])\n",
    "ev_dimming = float(pr.rtrace(view_pt, room_elgt.octree, params=tr_params.args()+[\"-pY\"], irradiance=True, header=False).decode().split()[0])\n",
    "ev_ec01 = float(pr.rtrace(view_pt, room_ec01.octree, params=tr_params.args()+[\"-pY\"], irradiance=True, header=False).decode().split()[0])\n",
    "wpi_tinted = get_ave_wpi(room.octree, grid, tr_params.args()+[\"-pY\"])\n",
    "wpi_dimming = get_ave_wpi(room_elgt.octree, grid, tr_params.args()+[\"-pY\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cc340c06",
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_dimming_factor(current_frame:int) -> float:\n",
    "    '''\n",
    "    Get dimming factor for the current frame \n",
    "\n",
    "    Parameters:\n",
    "    current_frame (int): The current frame number for which to calculate the dimming factor.\n",
    "\n",
    "    Returns:\n",
    "    float: The dimming factor for the given frame number.\n",
    "    '''\n",
    "    target_illum = 500\n",
    "    tvis60 = 0.677\n",
    "    tvis01 = 0.0108\n",
    "    tvis_i = interpolate_factor((tvis60, tvis01), current_frame, total_frames)\n",
    "    if wpi_tinted * tvis_i < target_illum:\n",
    "        dimming_factor = min((target_illum - wpi_tinted * tvis_i) / wpi_dimming, 1)\n",
    "    else:\n",
    "        dimming_factor = 0\n",
    "    return dimming_factor"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ed6b0be7",
   "metadata": {},
   "outputs": [],
   "source": [
    "os.makedirs(\"./Results\", exist_ok=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "24c723b1",
   "metadata": {},
   "outputs": [],
   "source": [
    "img_tinted = pr.rtpict(view, room.octree, nproc=8, xres=3200, yres=3200, params=render_params.args())\n",
    "with open(\"./Results/tinted.hsr\", \"wb\") as wtr:\n",
    "    wtr.write(img_tinted)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "003e38ee",
   "metadata": {},
   "outputs": [],
   "source": [
    "img_dimming = pr.rtpict(view, room_elgt.octree, nproc=8, xres=3200, yres=3200, params=render_params.args())\n",
    "with open(\"./Results/dimming.hsr\", \"wb\") as wtr:\n",
    "    wtr.write(img_dimming)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4dc72ad6",
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_dgp(hdr: bytes) -> float:\n",
    "    '''\n",
    "    Get DGP for the given hdr\n",
    "\n",
    "    Parameters:\n",
    "    hdr(bytes): The hdr to get the DGP for\n",
    "    \n",
    "    Returns:\n",
    "    float: The DGP for the given hdr\n",
    "    '''\n",
    "    evalglare_str = pr.evalglare(hdr, view=[\"-vta\", \"-vv\", \"180\", \"-vh\", \"180\"]).decode()\n",
    "\n",
    "    dgp = round(float(evalglare_str.split(\":\")[1].split(\" \")[1]),2)\n",
    "    return dgp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6428ff7f",
   "metadata": {},
   "outputs": [],
   "source": [
    "dgps = []\n",
    "timesteps = []\n",
    "for f in range(1, total_frames+1):\n",
    "    # get current time\n",
    "    if f == 1:\n",
    "        current_time = datetime.datetime(2024, 12, 21, 12)\n",
    "    else:\n",
    "        current_time += datetime.timedelta(seconds=1/3)\n",
    "    timesteps.append(current_time.strftime(\"%H:%M:%S.%f\")[:-5])\n",
    "    # combined tinted hdr and dimming hdr\n",
    "    tinted_factor = get_tinted_factors(f)[::-1] # reverse the order becauae channel split from red to blue\n",
    "    dimming_factor = get_dimming_factor(f)\n",
    "    tinted = pr.RcombInput(\"./Results/tinted.hsr\", scale=tinted_factor)\n",
    "    dimming = pr.RcombInput(\"./Results/dimming.hsr\", scale=[dimming_factor]*20)\n",
    "    combined_hdr = pr.rcomb([tinted, dimming], transform='RGB')\n",
    "    # downsample hdr\n",
    "    combined_ds_hdr = pr.pfilt(combined_hdr, xres='/3', yres='/3', gaussian_filter_radius=.6, one_pass=True)\n",
    "    with open(\"./Results/combined_ds.hdr\", \"wb\") as file:\n",
    "        file.write(combined_ds_hdr)\n",
    "\n",
    "    # falsecolor hdr\n",
    "    combined_ds_fc_hdr = pr.falsecolor(\"./Results/combined_ds.hdr\", contour=\"p\", scale='1e6', decades=6, palette=\"tbo\")\n",
    "    combined_ds_fc_ppm = pr.ra_ppm(combined_ds_fc_hdr)\n",
    "    with open(\"./Results/combined_ds_fc.ppm\", \"wb\") as file:\n",
    "        file.write(combined_ds_fc_ppm)\n",
    "    sp.run (f\"magick ./Results/combined_ds_fc.ppm ./Results/fc_{f:03d}.png\".split())\n",
    "    # _, yres, _, xres = pr.getinfo(combined_ds_hdr, dimension_only=True).decode().split()\n",
    "    # xres, yres = int(xres), int(yres)   \n",
    "    # img_data = pr.pvalue(combined_ds_hdr, header=False, outform='f', resstr=False)\n",
    "    # img_array = np.frombuffer(img_data, dtype=np.single).reshape(xres, yres, 3)\n",
    "    # luminance = img_array[:, :, 0] * 47.4 + img_array[:, :, 1] * 119.9 + img_array[:, :, 2] * 11.6\n",
    "    # plt.imshow(luminance, cmap='turbo', norm=mpl.colors.LogNorm(vmax=1e6, vmin=1))\n",
    "    # plt.colorbar()\n",
    "    # plt.axis('off')\n",
    "    # plt.savefig(f\"./Results/fc_{f:03d}.png\",bbox_inches='tight', dpi=500)\n",
    "    # plt.close()\n",
    "        \n",
    "    # get dgp\n",
    "    dgps.append(get_dgp(combined_ds_hdr))\n",
    "    # add human visual tone mapping\n",
    "    combined_ds_tm_hdr = pr.pcond(\"./Results/combined_ds.hdr\", human=True)\n",
    "    # convert to ppm\n",
    "    combined_ds_tm_ppm = pr.ra_ppm(combined_ds_tm_hdr)\n",
    "    with open(f\"./Results/combined_ds_tm.ppm\", \"wb\") as file:\n",
    "        file.write(combined_ds_tm_ppm)\n",
    "    # convert to png\n",
    "    sp.run(f\"magick ./Results/combined_ds_tm.ppm ./Results/{f:03d}.png\".split())\n",
    "    # remove temp files\n",
    "    sp.run(f\"rm ./Results/combined_ds_fc.ppm ./Results/combined_ds_tm.ppm ./Results/combined_ds.hdr\".split())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "31c49d61",
   "metadata": {},
   "outputs": [],
   "source": [
    "scene1 = pr.Scene(\"scene1\")\n",
    "sky_prims = pr.gendaylit(datetime.datetime(2024,12,21,12), 37, 122, 120, solar=True, dirnorm=800, diffhor=80).decode()\n",
    "for prim in pr.parse_primitive(sky_prims):\n",
    "    scene1.add_source(prim)\n",
    "scene1.add_source(pr.Primitive(\"skyfunc\", \"glow\", \"skyglow\", [], [1,1,1,0]))\n",
    "scene1.add_source(pr.Primitive(\"skyglow\", \"source\", \"skydome\", [], [0,0,1,180]))\n",
    "scene1.add_source(pr.Primitive(\"skyglow\", \"source\", \"groundplane\", [], [0,0,-1,180]))\n",
    "incident_irradiance = float(pr.rtrace(b\"0 -10 10 0 -1 0\", scene1.octree, header=False, params=[\"-ab\", \"1\", \"-I\"]).decode().split()[0])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f6adcbbc",
   "metadata": {},
   "outputs": [],
   "source": [
    "inside_environment = pywincalc.Environment(\n",
    "    air_temperature=294.15, \n",
    "    pressure=101325.0,\n",
    "    convection_coefficient = 0.0,\n",
    "    coefficient_model = pywincalc.BoundaryConditionsCoefficientModelType.CALCULATED_H,\n",
    "    radiation_temperature = 294.15,\n",
    "    emissivity = 1.0,\n",
    "    air_speed = 0.0,\n",
    "    air_direction = pywincalc.AirHorizontalDirection.NONE,\n",
    "    direct_solar_radiation = 0\n",
    ")\n",
    "\n",
    "outside_environment = pywincalc.Environment(\n",
    "    air_temperature=285.15, \n",
    "    pressure=101325.0,\n",
    "    convection_coefficient = 26.0,\n",
    "    coefficient_model = pywincalc.BoundaryConditionsCoefficientModelType.CALCULATED_H,\n",
    "    radiation_temperature = 255.15,\n",
    "    emissivity = 1.0,\n",
    "    air_speed = 5.0,\n",
    "    air_direction = pywincalc.AirHorizontalDirection.WINDWARD,\n",
    "    direct_solar_radiation = incident_irradiance\n",
    ")\n",
    "\n",
    "environmental_conditions = pywincalc.Environments(outside_environment, inside_environment)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "204a70b3",
   "metadata": {},
   "outputs": [],
   "source": [
    "ec01 = pywincalc.parse_json_file(\"./Materials/igsdb_product_7405.json\")\n",
    "solid_layers = [ec01]\n",
    "gs01 = pywincalc.GlazingSystem(solid_layers=solid_layers, environment=environmental_conditions)\n",
    "shgc_01 = gs01.shgc()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "05778082",
   "metadata": {},
   "outputs": [],
   "source": [
    "ec60 = pywincalc.parse_json_file(\"./Materials/igsdb_product_7406.json\")\n",
    "solid_layers = [ec60]\n",
    "gs_60 = pywincalc.GlazingSystem(solid_layers=solid_layers, environment=environmental_conditions)\n",
    "shgc_60 = gs_60.shgc() "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b9043f58",
   "metadata": {},
   "outputs": [],
   "source": [
    "l_ev = []\n",
    "l_mlux = []\n",
    "l_thml = []\n",
    "\n",
    "for f in np.arange(1, total_frames+1):\n",
    "    time = timesteps[f-1]\n",
    "    fig = plt.figure(figsize=(16, 9))\n",
    "\n",
    "    gs = fig.add_gridspec(nrows=4, ncols=3, \n",
    "                      width_ratios=[3, 2, 1], \n",
    "                      height_ratios=[1, 1, 1, 1],\n",
    "                    #   left=0.05, right=0.95, \n",
    "                    #   top=0.95, bottom=0.05, \n",
    "                    #   wspace=0, hspace=0\n",
    "                      )\n",
    "    \n",
    "    fig.subplots_adjust(left=0.15)\n",
    "\n",
    "    ax1 = fig.add_subplot(gs[:, 0])\n",
    "    ax1.imshow(imread(f\"./Results/{f:03d}.png\"))\n",
    "    ax1.set_anchor('N')\n",
    "    ax1.axis('off')\n",
    "\n",
    "    ax2 = fig.add_subplot(gs[:2, 1])\n",
    "    ax2.imshow(imread(f\"./Results/fc_{f:03d}.png\"))\n",
    "    ax2.axis('off')\n",
    "\n",
    "    ax3 = fig.add_subplot(gs[2:, 1])\n",
    "    ird_eye = []\n",
    "    tinted_factor = get_tinted_factors(f)\n",
    "    for it, tf, id in zip(ird_tinted, tinted_factor, ird_dimming):\n",
    "        ird_eye.append(it * tf + id * dimming_factor)\n",
    "    ax3.plot(spectrum, ird_eye, linewidth=3)\n",
    "    ax3.set(\n",
    "        xlabel=\"Wavelength (nm)\",\n",
    "        ylabel=r'Irradiance (W/m$^2$)',\n",
    "        ylim=(0, 225)\n",
    "    )\n",
    "    ax3.xaxis.set_major_locator(plt.MaxNLocator(5))\n",
    "    ax3.yaxis.set_major_locator(plt.MaxNLocator(5))\n",
    "    \n",
    "    ax4 = fig.add_subplot(gs[0, 2])\n",
    "    dimming_factor = get_dimming_factor(f)\n",
    "    ev_tf = interpolate_factor((ev_tinted, ev_ec01), f, total_frames)\n",
    "    ev = int(ev_tinted * ev_tf + ev_dimming * dimming_factor)\n",
    "    l_ev.append(ev)\n",
    "    ax4.plot(np.arange(1, f+1)/3, l_ev, linewidth=3)\n",
    "    ax4.set(xlim=(0, 60), ylim=(0, 37500), ylabel=\"Ev (lux)\")\n",
    "    ax4.yaxis.get_major_formatter().set_powerlimits((0, 1))\n",
    "\n",
    "    ax5 = fig.add_subplot(gs[1, 2])\n",
    "    mlux_tf = interpolate_factor((mlux_tinted, mlux_ec01), f, total_frames)\n",
    "    mlux = int(mlux_tinted * mlux_tf + mlux_dimming * dimming_factor)\n",
    "    l_mlux.append(mlux)\n",
    "    ax5.plot(np.arange(1, f+1)/3, l_mlux, linewidth=3)\n",
    "    ax5.set(xlim=(0, 60), ylim=(0, 37500), ylabel=\"Mlux (EML)\")\n",
    "    #scientific notation\n",
    "    ax5.yaxis.get_major_formatter().set_powerlimits((0, 1))\n",
    "\n",
    "    ax6 = fig.add_subplot(gs[2, 2])\n",
    "    ax6.plot(np.arange(1, f+1)/3, dgps[:f], linewidth=3)\n",
    "    ax6.set(xlim=(0, 60), ylim=(0, 1.2), ylabel=\"DGP\")\n",
    "\n",
    "    ax7 = fig.add_subplot(gs[3, 2])\n",
    "    sghc = interpolate_factor((shgc_60, shgc_01), f, total_frames)\n",
    "    thml = int(sghc * incident_irradiance)\n",
    "    l_thml.append(thml)\n",
    "    ax7.plot(np.arange(1, f+1)/3, l_thml, linewidth=3)\n",
    "    ax7.set(xlim=(0, 60), ylim=(0, 1000), xlabel=\"Time (s)\", ylabel=\"Thermal load (w/m2)\")\n",
    "    ax7.yaxis.get_major_formatter().set_powerlimits((0, 1))\n",
    "\n",
    "    text = f\"Time: {time}\\nPhotopic illuminance: {ev} lux\\nMelanopic illuminance: {mlux} EML\\nDGP: {dgps[f-1]}\\nThermal load: {thml} w/m2\"\n",
    "    ax1.text(0, -0.03, f\"{text}\", fontsize=20, transform=ax1.transAxes, ha='left', va='top', bbox=dict(facecolor='white', alpha=0.8, edgecolor='none'))\n",
    "\n",
    "    fig.savefig(f\"./Results/compile_{f:03d}.png\", dpi=500)\n",
    "\n",
    "    plt.close(fig)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "33ae042b",
   "metadata": {},
   "outputs": [],
   "source": [
    "!ffmpeg -framerate 3 -i Results/compile_%03d.png -c:v libx264 -preset slow -crf 18 -vf \"scale=800:800,format=yuv420p\" final.mp4"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1506e04c",
   "metadata": {},
   "source": [
    "# extra data visualization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "692cacb9",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "for i in np.arange(180):\n",
    "    tf = get_tinted_factors(i)\n",
    "    ax.plot(spectrum, tf * spectral_transmittance_60, color = \"0.5\")\n",
    "    ax.set_ylabel(\"Visible Transmittance\")\n",
    "    ax.set_xlabel(\"Wavelength [nm]\")\n",
    "    ax.set_title(\"Spectral Transmittance Distribution Per Electrochromic Tinted State\")\n",
    "fig.savefig(\"spectral_transmittance.png\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e87e2724",
   "metadata": {},
   "outputs": [],
   "source": [
    "tfs = []\n",
    "for f in np.arange(1, total_frames+1):\n",
    "    tint_factors = get_tinted_factors(f)\n",
    "    tfs.append(tint_factors)\n",
    "fig, ax = plt.subplots()\n",
    "for n, sp in enumerate(spectrum):\n",
    "    if sp == 490:\n",
    "        color = \"blue\"\n",
    "        label = \"480-500 nm\"\n",
    "    else:\n",
    "        color = \"lightgray\"\n",
    "        label = \"370 - 780 nm\"\n",
    "    ax.plot(np.arange(180)/3, [tf[n] for tf in tfs], color=color, label=label)\n",
    "\n",
    "label, legend = ax.get_legend_handles_labels()\n",
    "ax.legend(label[5:7], legend[5:7])\n",
    "\n",
    "ax.set_xlabel(\"Time [second]\")\n",
    "ax.set_ylabel(\"Normalized Tinted Factors\")\n",
    "ax.set_title(\"Electrochromic Tinted Factors Per Wavelength\")\n",
    "\n",
    "fig.savefig(\"tinted_factors.png\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e3515683",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "ax.plot(spectrum, ird_tinted, label=\"tvis: 60%\")\n",
    "for f in np.arange(2, total_frames+1):\n",
    "    if f == 180:\n",
    "        color = \"C1\"\n",
    "        label = \"tvis: 1%\"\n",
    "    else:\n",
    "        color = \"lightgray\"\n",
    "        label = \"interpolated\"\n",
    "    ax.plot(spectrum, ird_tinted*get_tinted_factors(f), color=color, label=label)\n",
    "label, legend = ax.get_legend_handles_labels()\n",
    "ax.legend(label[:1]+label[-2:], legend[:1]+legend[-2:], frameon=False)\n",
    "ax.set(ylabel=r\"Irradiance [cd/m$^2$]\", xlabel=\"Wavelength [nm]\")\n",
    "\n",
    "fig.savefig(\"irradiance.png\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "75062232",
   "metadata": {},
   "outputs": [],
   "source": [
    "ev = []\n",
    "mlux = []\n",
    "fig, ax = plt.subplots()\n",
    "for f in np.arange(1, total_frames+1):\n",
    "    dimming_factor = get_dimming_factor(f)\n",
    "    ev_tf = interpolate_factor((ev_tinted, ev_ec01), f, total_frames)\n",
    "    mlux_tf = interpolate_factor((mlux_tinted, mlux_ec01), f, total_frames)\n",
    "    ev.append(int(ev_tinted * ev_tf + ev_dimming * dimming_factor))\n",
    "    mlux.append(int(mlux_tinted * mlux_tf + mlux_dimming * dimming_factor))\n",
    "ax.plot(np.arange(1, total_frames+1)/3, ev, label=\"ev\")\n",
    "ax1 = ax.twinx()\n",
    "ax1.plot(np.arange(1, total_frames+1)/3, mlux, label=\"mlux\", color=\"C1\")\n",
    "# ax.plot(np.arange(1, total_frames+1), [500]*total_frames, linestyle=\"dashed\", color=\"lightgray\", label = r\"target wpi = 500 lux$\")\n",
    "\n",
    "ax.set_xlabel(xlabel=\"Time [second]\")\n",
    "ax.set_ylabel(ylabel=\"Vertical Illuminance [lux]\", color=\"C0\")\n",
    "ax1.set_ylabel(ylabel=\"Melanopic lux [EML]\", color=\"C1\")\n",
    "ax.tick_params(axis='y', colors='C0')\n",
    "ax1.tick_params(axis='y', colors='C1')\n",
    "\n",
    "fig.tight_layout()\n",
    "fig.savefig(\"ev_mlux.png\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "96987400",
   "metadata": {},
   "outputs": [],
   "source": [
    "# wpi_wo_lgt = []\n",
    "# wpi_lgt = []\n",
    "# dimming_factors= []\n",
    "# target_illum = 500\n",
    "# tvis60 = 0.677\n",
    "# tvis01 = 0.0108\n",
    "\n",
    "# fig, ax = plt.subplots()\n",
    "# for f in np.arange(1, total_frames+1):\n",
    "#     dimming_factor = get_dimming_factor(f)\n",
    "#     dimming_factors.append(dimming_factor)\n",
    "#     tvis_i = interpolate_factor((tvis60, tvis01), f, total_frames)\n",
    "#     wpi_wo_lgt.append(wpi_tinted * tvis_i)\n",
    "#     wpi_lgt.append(wpi_tinted * tvis_i + dimming_factor * wpi_dimming)\n",
    "# ax.plot([500]*60, linestyle=\"dashed\", color=\"lightgray\", label = r\"target wpi = 500 lux$\")\n",
    "# ax.plot(np.arange(1, total_frames+1)/3, wpi_wo_lgt, label=\"without electric light\")\n",
    "# ax.plot(np.arange(1, total_frames+1)/3, wpi_lgt, label=\"with electric light\", linestyle=\"dashed\")\n",
    "\n",
    "# ax1 = ax.twinx()\n",
    "# ax1.plot(np.arange(1, total_frames+1)/3, dimming_factors, label=\"dimming factor\", color=\"C1\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c5a4ce80",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": ".venv",
   "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.11.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}