moon_brightness_temp_sim_v3.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "9cdb7fc2-9efa-4148-9827-a7bdda9287ec",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "# Average Lunar Flux Brightness Simulation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "bd44c3d2-8f58-415e-8752-187a71c8fd80",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "# import various packages\n",
    "import numpy as np\n",
    "# import cupy as cp"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "732da333-75fe-4344-a383-2cff0c1152ae",
   "metadata": {},
   "source": [
    "## Inititalize Facet Regions"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fe0c7a23-dc02-421b-8526-36c5987e7735",
   "metadata": {},
   "source": [
    "Let's split up the moon into facets along lines of latitude and longitude "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3f5d2885-149e-463c-9d3d-7a647dcf3798",
   "metadata": {},
   "source": [
    "Only necessary to structure the latitude from 0 to 90 degrees because of the symmetry of the northern/southern hemispheres, will need to multiply the average flux density by 2 at the end of the simulation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "5a468dd8-782d-4945-b0e8-340e1b192f08",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define ranges for phase, longitude, latitude, and declination\n",
    "phase_array = np.arange(0, 360, dtype = np.int16) # angular positions of the sun overhead w.r.t latitude = 0, longitude = 0 position\n",
    "longitude_array= np.arange(-180, 181, dtype = np.int16)  # setting longitude range in degrees\n",
    "latitude_array = np.arange(0, 91, dtype = np.int16)  # setting latitude range in degrees, (0,90 exclusive) I'm considering the centre point of each facet\n",
    "declination = 0  # The angular position of the sun at solar noon in degrees"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "299d443e-cd59-4887-a582-00255a4e3daa",
   "metadata": {},
   "source": [
    "## Measure Facets Areas"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "ac5eb207-b21d-4d19-9b82-ef794e983109",
   "metadata": {},
   "outputs": [],
   "source": [
    "def lunar_facets_areas(latitude_array, longitude_array, moon_radius):\n",
    "    \"\"\"\n",
    "    Generation of facet areas over the surface of the Moon using latitude and longitude arrays.\n",
    "    \n",
    "    - Parameters:\n",
    "        - latitude_array: Array of latitudes in degrees.\n",
    "        - longitude_array: Array of longitudes in degrees.\n",
    "        - moon_radius: Radius of the Moon in meters.\n",
    "    \n",
    "    - Returns: facet_areas: 2D array of facet areas (m²) corresponding to each latitude and longitude grid cell.\n",
    "    \"\"\"\n",
    "    \n",
    "    # Convert latitude and longitude arrays from degrees to radians\n",
    "    latitudes = np.radians(latitude_array)\n",
    "    longitudes = np.radians(longitude_array)\n",
    " \n",
    "    # Compute the areas of the facets\n",
    "    facet_areas = (moon_radius**2) * np.abs(-1* np.sin(latitudes[np.newaxis, 1:]) + np.sin(latitudes[np.newaxis, :-1])) * (longitudes[1:, np.newaxis] - longitudes[:-1, np.newaxis])\n",
    "    \n",
    "    return facet_areas"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fe99758a-9fb1-4141-9719-8dfddb8a4ef4",
   "metadata": {},
   "source": [
    "### Test of generate_moon_facets_areas\n",
    "If the radius is 1 the sum of all facet areas should be 2pi for a hemisphere"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "c09d5b4c-4663-47b6-b1f8-85bc38a62f72",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "shape =  (360, 90) area of surface =  6.283185 , Difference: 2pi - area_of_facets =  3.0199159795074593e-07\n"
     ]
    }
   ],
   "source": [
    "moon_radius_test = 1  # Radius of the Moon in meters\n",
    "area_of_facets_test = lunar_facets_areas(latitude_array, longitude_array, moon_radius_test)\n",
    "print('shape = ', area_of_facets_test.shape, 'area of surface = ' , np.sum(area_of_facets_test), ', Difference: 2pi - area_of_facets = ', 2 * np.pi - np.sum(area_of_facets_test))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3ae12729-64d7-458e-bdcf-beba5904a032",
   "metadata": {},
   "source": [
    "### Implement the funciton with the radius of the Moon in meters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "5d890a18-cba2-4bf6-aba3-8ba6d8b94f4d",
   "metadata": {},
   "outputs": [],
   "source": [
    "moon_radius = 1.737e6  # Radius of the Moon in meters\n",
    "area_of_facets = lunar_facets_areas(latitude_array, longitude_array, moon_radius)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "46473022-790c-489c-9131-e2293624564e",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(66338.2734375, 257.5621739260251)"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.abs(2*np.pi*moon_radius**2 - np.sum(area_of_facets)), np.sqrt(np.abs(2*np.pi*moon_radius**2 - np.sum(area_of_facets)))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6ac264ad-6b77-4897-af4f-29bbccc567cb",
   "metadata": {},
   "source": [
    "The area is off by a bit, how much? a grid approximately 257 m * 257 m in area, I think thats acceptable."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ebad241b-35b9-4706-84e6-611cc6bc33fb",
   "metadata": {},
   "source": [
    "## Solar Zenith Angle for Phase cycle"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "20bd09bd-3101-4d39-82c2-a4d17c88c653",
   "metadata": {},
   "outputs": [],
   "source": [
    "# def calculate_zenith_angle(phase, latitude, longitude, declination):\n",
    "#     \"\"\"\n",
    "#     Vectorized calculation of the zenith angle on each facet of lat, lon.\n",
    "#     - Parameters: \n",
    "#         - Latitude: a numpy array of latitudes\n",
    "#         - Longutide: a numpy array of longitudes\n",
    "#         - Phase: a numpy array of phases\n",
    "#     - Returns: Outputs a 3D array of zenith angles with shape (360, 360, 91)\n",
    "#     \"\"\"\n",
    "#     # Convert angles from degrees to radians for numpy trigonometric functions\n",
    "#     phase_rad = np.radians(phase)[:, np.newaxis, np.newaxis]  # shape (360, 1, 1)\n",
    "#     longitude_rad = np.radians(longitude)[np.newaxis, :, np.newaxis]  # shape (1, 360, 1)\n",
    "#     latitude_rad = np.radians(latitude)[np.newaxis, np.newaxis, :]  # shape (1, 1, 91)\n",
    "#     declination_rad = np.radians(declination)\n",
    "\n",
    "#     # Calculate the zenith angle for all combinations of phase, lon, and lat\n",
    "#     zenith_angle_array = np.arccos(\n",
    "#         np.sin(declination_rad) * np.sin(latitude_rad)\n",
    "#         + np.cos(declination_rad) * np.cos(latitude_rad)\n",
    "#         * np.cos(longitude_rad - phase_rad)\n",
    "#     )\n",
    "\n",
    "#     return zenith_angle_array"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "2d0e393f-5deb-4c01-ae3f-24aac0fd85e1",
   "metadata": {},
   "outputs": [],
   "source": [
    "def calculate_zenith_angle(phase, latitude, longitude, declination):\n",
    "    \"\"\"\n",
    "    Vectorized calculation of the zenith angle on each facet of lat, lon.\n",
    "    - Parameters: \n",
    "        - Latitude: a numpy array of latitudes\n",
    "        - Longutide: a numpy array of longitudes\n",
    "        - Phase: a numpy array of phases\n",
    "    - Returns: Outputs a 3D array of zenith angles with shape (360, 360, 91)\n",
    "    \"\"\"\n",
    "    # Convert angles from degrees to radians for numpy trigonometric functions\n",
    "    # phase_rad = np.radians(phase)[:, np.newaxis, np.newaxis]  # shape (360, 1, 1)\n",
    "    longitude_rad = np.radians(longitude)[:, np.newaxis]  # shape (1, 360, 1)\n",
    "    latitude_rad = np.radians(latitude)[np.newaxis, :]  # shape (1, 1, 91)\n",
    "    declination_rad = np.radians(declination)\n",
    "\n",
    "    # Calculate the zenith angle for all combinations of phase, lon, and lat\n",
    "    zenith_angle_array = np.arccos(\n",
    "        np.sin(declination_rad) * np.sin(latitude_rad)\n",
    "        + np.cos(declination_rad) * np.cos(latitude_rad)\n",
    "        * np.cos(longitude_rad)\n",
    "    )\n",
    "\n",
    "    return zenith_angle_array"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "0d325566-c62c-48fc-a913-5c325e454b32",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(361, 91)\n"
     ]
    }
   ],
   "source": [
    "# Call the function\n",
    "zenith_angle_array = calculate_zenith_angle(phase_array,\n",
    "                                            latitude_array,\n",
    "                                            longitude_array,\n",
    "                                            declination)\n",
    "\n",
    "# check the shape\n",
    "print(zenith_angle_array.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "96904bf8-c655-4c29-8434-51b26a0e7350",
   "metadata": {
    "jp-MarkdownHeadingCollapsed": true,
    "tags": []
   },
   "source": [
    "#### zenith_angle sanity checks"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "67560c9b-7d3c-4ad3-8ce7-17c83b8b6413",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# np.degrees(zenith_angle_array[0,180,:]) \n",
    "# sanity check for what the latitude angles would look like for a sub solar position at theta = 0, longitude, = 0, and phase = 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "5505d00d-207a-417c-b8a5-1b7b31417529",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# np.degrees(zenith_angle_array[0,225,:]) \n",
    "# sanity check for what the latitude angles would look like for longitude = 45, and phase = 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "f0f5cc0f-deff-429e-b2c4-5432a9cad2dc",
   "metadata": {},
   "outputs": [],
   "source": [
    "# np.degrees(zenith_angle_array[45,225,:]) \n",
    "# sanity check for what the latitude angles would look like for longitude = 45, and phase = 45"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "09204aa3-0f8d-484e-a5c1-58231080629e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# np.degrees(zenith_angle_array[0,270,:]) \n",
    "# sanity check for what the latitude angles would look like for longitude = 90, and phase = 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "fffa9f4d-3750-4fb5-b363-884911e93a7f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# np.degrees(zenith_angle_array[0,0,:]) \n",
    "# sanity check for what the latitude angles would look like for longitude = 0, and phase = 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "7d1f05fa-569f-46e4-bd47-b8fc6e5d4750",
   "metadata": {},
   "outputs": [],
   "source": [
    "# vector_zenith_angle_array = vector_calculate_zenith_angle(phase_array,\n",
    "#                                                           latitude_array,\n",
    "#                                                           longitude_array,\n",
    "#                                                           declination)\n",
    "# zenith_angle_array.shape # checking shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "b4bee4d3-e3a1-4711-be60-64a603cb2293",
   "metadata": {},
   "outputs": [],
   "source": [
    "# np.sum(vector_zenith_angle_array) == np.sum(zenith_angle_array) \n",
    "# because they are the essentially the same function, they should result in the same sum of angles"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "525a9635-a73f-4134-adb9-aeb6c7839391",
   "metadata": {},
   "outputs": [],
   "source": [
    "# zenith_angle_array[1,:,:].shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e3334eae-03bb-4034-add0-63ea82e4ad0d",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "5d328692-cbdc-4558-8610-f2ceb5be9c47",
   "metadata": {},
   "source": [
    "## Calculate Heat Input"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "611f489d-fa5b-4248-b1f7-2b53eadb014c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# def flux_input_func(W_0, zenith_angle_array, albedo, phase_step):\n",
    "#     \"\"\"\n",
    "#     - Parameters:\n",
    "#         - W_0: energy flux W / m**2\n",
    "#         - zenith_angle_array: array of incident vector angle w.r.t sun\n",
    "#         - albedo: Lunar Albedo\n",
    "#         - phase_step: phase ranges from [0,360) degrees\n",
    "    \n",
    "#     -Return: the flux input from the Sun. The time is indexed in degrees of phase, which allows isolating the values\n",
    "#     of latitude and longitude at that specific phase.\n",
    "#     \"\"\"\n",
    "#     # Calculate the heat input based on the zenith angle\n",
    "#     W_t = W_0 * np.clip(np.cos(zenith_angle_array[phase_step, :, :]), 0, None) * (1 - albedo)\n",
    "#     # the clip returns a 0 for any values less than 0 resulting from an angle greater than 90 degrees\n",
    "#     # this way no negative heat is introduced.\n",
    "    \n",
    "#     return W_t\n",
    "               \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "b9a0ad1c-f4dd-4337-a369-206e7de7790a",
   "metadata": {},
   "outputs": [],
   "source": [
    "def flux_input_func(W_0, zenith_angle_array, albedo):\n",
    "    \"\"\"\n",
    "    - Parameters:\n",
    "        - W_0: energy flux W / m**2\n",
    "        - zenith_angle_array: array of incident vector angle w.r.t sun\n",
    "        - albedo: Lunar Albedo\n",
    "        - phase_step: phase ranges from [0,360) degrees\n",
    "    \n",
    "    -Return: the flux input from the Sun. The time is indexed in degrees of phase, which allows isolating the values\n",
    "    of latitude and longitude at that specific phase.\n",
    "    \"\"\"\n",
    "    # Calculate the heat input based on the zenith angle\n",
    "    W_t = W_0 * np.clip(np.cos(zenith_angle_array), 0, None) * (1 - albedo)\n",
    "    # the clip returns a 0 for any values less than 0 resulting from an angle greater than 90 degrees\n",
    "    # this way no negative heat is introduced.\n",
    "    \n",
    "    return W_t\n",
    "               \n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "125f5744-acd3-4f3b-83c5-128b3bae43bc",
   "metadata": {},
   "source": [
    "#### Solar Irradiance Test"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "1907d62e-a426-49d4-a106-1eb0b8713f9e",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "W_0 = 1365.0 # Solar irradiance in W.m-2, the energy flux arriving from the Sun\n",
    "albedo = 0.15 # the lunar albedo\n",
    "phase_step = 0 # One degree in phase corresponds to 6700 seconds, phase ranges from [0,360) degrees\n",
    "                # corresponds to [0, 2,405,300 seconds), [0, 668 hours), [0,27.8 days) phase_step inputs must be integers\n",
    "\n",
    "flux_input_array = flux_input_func(W_0, \n",
    "                                   zenith_angle_array, \n",
    "                                   albedo)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "6f8dc3b0-6de9-4ec6-ad32-6158633c0394",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(361, 91)"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# shape check\n",
    "flux_input_array.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "3630add7-39f5-4cf8-94ad-2ef279cba5ca",
   "metadata": {
    "editable": true,
    "scrolled": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "function output 1160.25 manual calculation 1160.25\n"
     ]
    }
   ],
   "source": [
    "# verify that its matching the expected value\n",
    "print('function output', flux_input_array[180,:][0], 'manual calculation', (1-albedo) * 1365)\n",
    "# sanity check for when longitude = 0, latitude = 0\n",
    "# if time = 0 then subsolar position is theta = 0, longitude = 0"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6580d84a-3c6b-452a-a210-60dd2a7c6a12",
   "metadata": {},
   "source": [
    "### Test of Lunar Flux"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "98b32026-896b-4da2-9a98-de66f3c13c1e",
   "metadata": {},
   "source": [
    "Combining the functions from flux input and facet areas"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "3dbff62e-809c-4b6d-8129-61dd273e00f1",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "((360, 90), (361, 91))"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "area_of_facets.shape, flux_input_array.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dcf2363a-0b2d-486a-a408-533116d57980",
   "metadata": {},
   "source": [
    "If we are instead interested in the centres of the facets, the flux input needs its ends trimed off the last values in the longitude and latitude values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "1ba30143-e380-4ffe-b634-58bc6955003b",
   "metadata": {},
   "outputs": [],
   "source": [
    "lunar_total_flux = area_of_facets * flux_input_array[:-1,:-1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "3053c161-2944-4edf-8530-b3e4037fc301",
   "metadata": {},
   "outputs": [],
   "source": [
    "# # test of the 0 longitude positions\n",
    "# lunar_total_flux[180]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "89682299-a40e-46ee-afc2-3f88da2ba8e1",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5528970500000000.0"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.sum(lunar_total_flux)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "91d9d173-1493-4759-b641-f116561a29ca",
   "metadata": {},
   "source": [
    "## Define Thermal Energy Tranfer Functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "3f318d55-7e0f-4dd0-a1c6-aebcfc513c3f",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def rate_of_change_surface_thermal_energy_func(flux_input_array, area_of_facets, T_0, T_1, sigma, BasCond, tau):\n",
    "    \"\"\"\n",
    "    A function to calculate the rate of change of thermal energy at the surface \n",
    "    - Parameter: \n",
    "        - heat_input_array: a multi-dimensional array\n",
    "        - T_0: temperatue at surface position\n",
    "        - T_1: temperature at depth below surface\n",
    "        - sigma: Stefan's constant 5.67e-8 (W/m²K⁴)\n",
    "        - BasCond: 1500, Basalt conductivity in W.m-3.K-1\n",
    "        - tau: 0.01, is the thickness of the slice [m]\n",
    "        \n",
    "    -Return: A scalar value of the rate of change of thermal energy of the surface layer\n",
    "    \"\"\"\n",
    "    # previous version\n",
    "    # dQ0_dt = heat_input_array - sigma * T_0**4 + (T_1 - T_0) / R \n",
    "\n",
    "    # Ken's Version\n",
    "    # LTK[j] =  LTK[j] + (A*(I - sigma*LTK[j]**4) + (LTK[j+1] - LTK[j])*beta)*delta_t/gamma        \n",
    "\n",
    "    # New Version\n",
    "    dQ0_dt = area_of_facets * (flux_input_array - sigma * T_0**4) + ((T_1 - T_0) * BasCond * area_of_facets/tau)\n",
    "\n",
    "    return dQ0_dt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "id": "26e04752-20b0-4c2b-96a5-5274835977b1",
   "metadata": {},
   "outputs": [],
   "source": [
    "# beta = 0.15         #BasCond * A/delta_z         #conductivity parameter\n",
    "# gamma = 10050 #BasCap * m                 #thermal capacity parameter\n",
    "# BasCond = 1.50e3                 #Basalt conductivity in W.m-3.K-1\n",
    "# BasCap  = 670.0                  #Basalt thermal capacity in J Kg-1.K-1 \n",
    "# m = A * delta_z * rho            #mass of element\n",
    "# rho     = 1500.0                 #Lunar regolith density in Kg.m-3\n",
    "# delta_z = 0.01                   #thickness of element (metres)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "id": "9244822e-d06b-4400-8f5a-30d3731ee39f",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def rate_of_change_intermediate_thermal_energy_func(T_np1, T_n, T_nm1, BasCond, area_of_facets, tau):\n",
    "    \"\"\"\n",
    "    A function to calculate the rate of change of thermal energy in the intermediate layers\n",
    "    - Parameter: \n",
    "        - T_nm1: Previous Layer (T_n-1)\n",
    "        - T_np1: Next Layer (T_n+1)\n",
    "        - T_n: Layer being evaluated\n",
    "        - beta:\n",
    "        \n",
    "    -Return: A scalar value of the rate of change of thermal energy of the surface layer applied to each element for which there \n",
    "    \"\"\"\n",
    "    # previous version\n",
    "    # dQn_dt = ((T_np1 - T_n) / R) + ((T_nm1 - T_n) / R)\n",
    "\n",
    "    # Ken's Version from moon_temp_model\n",
    "    # LTK[j] = LTK[j] + ((LTK[j+1]-LTK[j])*beta + (LTK[j-1] - LTK[j])*beta)*delta_t/gamma \n",
    "\n",
    "    dQn_dt = ((T_np1 - T_n) + (T_nm1 - T_n)) * (BasCond * area_of_facets/tau)\n",
    "\n",
    "    \n",
    "    return dQn_dt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "id": "b56fffa3-4e63-4305-ac1d-07eb4a6b5de5",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def rate_of_change_final_thermal_energy_func(T_const, T_f, T_nm1, BasCond, area_of_facets, tau):\n",
    "    \"\"\"\n",
    "    The rate of change of thermal energy at the final depth \n",
    "    - Parameter: \n",
    "        - T_nm1: Previous Layer (T_n-1)\n",
    "        - T_const: The temperature below the layers being evaluated, a constant value\n",
    "        - T_f: The final layer being evaluated that between the previous layer and the T_constant layer\n",
    "        - R: Thermal Resistance\n",
    "    -Return: A scalar value of the rate of change of thermal energy of the surface layer applied to each element for which there \n",
    "\n",
    "    \"\"\"\n",
    "    dQf_dt = ((T_const - T_f) + (T_nm1 - T_f)) * (BasCond * area_of_facets/tau)\n",
    "    \n",
    "    return dQf_dt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "id": "64c5237b-7c0c-4eea-a908-e1c3ba7d3a75",
   "metadata": {},
   "outputs": [],
   "source": [
    "# old version\n",
    "# def thermal_capacitance(beta, area, tau):\n",
    "#     \"\"\"\n",
    "#     Calculation of the thermal capacitance\n",
    "#     - Parameters:\n",
    "#         - beta: the thermal capacity per unit volume\n",
    "#         - area: is the area of the slice of material\n",
    "#         - tau: is the thickness of the slice\n",
    "#     - Return: Q, thermal capacitance [J/K]\n",
    "#     \"\"\"\n",
    "#     Q = beta * area * tau\n",
    "    \n",
    "#     return Q"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "id": "beb368aa-64b1-48b6-8a6a-8e2df85cf763",
   "metadata": {},
   "outputs": [],
   "source": [
    "# New Version\n",
    "def thermal_capacitance(BasCap, area, tau, density):\n",
    "    \"\"\"\n",
    "    Calculation of the thermal capacitance\n",
    "    - Parameters:\n",
    "        - BasCap: Basalt thermal capacity in J Kg-1.K-1 \n",
    "        - area: is the area of the slice of material, passed as an array of area's\n",
    "        - tau: 0.01, is the thickness of the slice [m]\n",
    "        - density: 1500.0, Lunar regolith density in Kg.m-3\n",
    "    - Return: gamma, thermal capacitance [J/K]\n",
    "    \"\"\"\n",
    "    gamma = BasCap * area * tau * density\n",
    "    \n",
    "    return gamma"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "id": "de8755a0-8c69-4ba3-b0b3-69174055b43a",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# not actually used in the function, simplified to output temperature\n",
    "\n",
    "# def thermal_charge(Q_i, dQi_dt, dt):\n",
    "#     \"\"\"\n",
    "#     The thermal charge, Q in one of the layers is updated where,\n",
    "#     Q_i is the presently evaluated thermal charge,\n",
    "#     dQi_dt is the rate of change of thermal energy,\n",
    "#     dt is the time interval\n",
    "#     \"\"\"\n",
    "#     Q_i += dQi_dt * dt"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "16e32d30-294d-4b03-9b77-b865008c0d49",
   "metadata": {},
   "source": [
    "## Develop the Heat Transfer Model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "id": "436553e8-f9f4-4efc-8fb1-e919a2030f1f",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# # Heat transfer model function\n",
    "# def heat_transfer_model_func(heat_transfer_array, W_0, albedo, T_const, BasCap, area_of_facets, tau, rho, delta_t, model_run_time):\n",
    "#     \"\"\"\n",
    "#     Heat transfer model to update temperature in each layer over time.\n",
    "\n",
    "#     Parameters:\n",
    "#     - heat_transfer_array: 3D NumPy array representing solar flux with dependance on incident angle\n",
    "#     - W_0: Energy flux arriving from the Sun\n",
    "#     - albedo: Albedo (reflectivity) of the surface\n",
    "#     - T_const: Constant temperature for the final boundary condition\n",
    "#     - R: Thermal resistance between layers\n",
    "#     - beta\n",
    "#     - area\n",
    "#     - tau\n",
    "#     - delta_t: Time step for updating temperatures\n",
    "#     - model_run_time: Total simulation time in seconds\n",
    "\n",
    "#     Returns:\n",
    "#     - Updated heat_transfer_array with new temperatures for each layer.\n",
    "#     \"\"\"\n",
    "#     phase_step = 0\n",
    "    \n",
    "#     for second in range(0, model_run_time, delta_t): # range(start,stop,step)\n",
    "        \n",
    "#         # Step 1: Initialize Thermal Capacitance\n",
    "#         gamma = thermal_capacitance(BasCap, area_of_facets, tau, rho)\n",
    "\n",
    "        \n",
    "#         # Step 2: Update the surface layer (T_0)\n",
    "#         T_0 = heat_transfer_array[phase_step,:,:,0]\n",
    "#         T_1 = heat_transfer_array[phase_step,:,:,1]\n",
    "#         W_t = flux_input_func(W_0, zenith_angle_array, albedo, phase_step)\n",
    "#         dQ0_dt = rate_of_change_surface_thermal_energy_func(W_t[:-1,:-1], area_of_facets, T_0, T_1, sigma, BasCond, tau)\n",
    "#         T_0 += dQ0_dt * delta_t / gamma\n",
    "\n",
    "\n",
    "#         # Step 3: Update intermediate layers\n",
    "#         T_n = heat_transfer_array[phase_step,:,:,1:-1]  # All intermediate layers\n",
    "#         T_np1 = heat_transfer_array[phase_step,:,:,2:]  # Next layer in each direction\n",
    "#         T_nm1 = heat_transfer_array[phase_step,:,:,:-2]  # Previous layer in each direction\n",
    "#         dQn_dt = rate_of_change_intermediate_thermal_energy_func(T_np1, T_n, T_nm1, BasCond, area_of_facets[:,:,np.newaxis], tau)\n",
    "#         T_n += dQn_dt * delta_t / gamma[:,:,np.newaxis]\n",
    "\n",
    "#         # Step 4: Update the final layer (T_f)\n",
    "#         T_f = heat_transfer_array[phase_step,:,:,-1]\n",
    "#         T_fm1 = heat_transfer_array[phase_step,:,:,-2]\n",
    "#         T_nm1 = T_nm1[:,:,-1]\n",
    "#         dQf_dt = rate_of_change_final_thermal_energy_func(T_const, T_f, T_nm1, BasCond, area_of_facets, tau)\n",
    "#         T_f += dQf_dt * delta_t / gamma\n",
    "\n",
    "#         # Step 5: Update the heat transfer array with new temperatures\n",
    "#         heat_transfer_array[phase_step,:,:,0] = T_0\n",
    "#         heat_transfer_array[phase_step,:,:,1:-1] = T_n\n",
    "#         heat_transfer_array[phase_step,:,:,-1] = T_f\n",
    "\n",
    "#         # uncomment for diagnosing\n",
    "#         # print(f\"At second {second}: T_0 = {T_0}, T_1 = {T_1}, W_t = {W_t}, dQ0_dt = {dQ0_dt}\")\n",
    "#             # Independent phase_step increment based on the condition\n",
    "#         # if second % 6700 == 0:\n",
    "#         #     phase_step += 1\n",
    "#     return heat_transfer_array\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "id": "1d7b6041-925c-4430-bc09-1723d962d467",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Heat transfer model function\n",
    "def heat_transfer_model_func(heat_transfer_array, W_0, albedo, T_const, BasCap, area_of_facets, tau, rho, delta_t, model_run_time):\n",
    "    \"\"\"\n",
    "    Heat transfer model to update temperature in each layer over time.\n",
    "\n",
    "    Parameters:\n",
    "    - heat_transfer_array: 3D NumPy array representing solar flux with dependance on incident angle\n",
    "    - W_0: Energy flux arriving from the Sun\n",
    "    - albedo: Albedo (reflectivity) of the surface\n",
    "    - T_const: Constant temperature for the final boundary condition\n",
    "    - R: Thermal resistance between layers\n",
    "    - beta\n",
    "    - area\n",
    "    - tau\n",
    "    - delta_t: Time step for updating temperatures\n",
    "    - model_run_time: Total simulation time in seconds\n",
    "\n",
    "    Returns:\n",
    "    - Updated heat_transfer_array with new temperatures for each layer.\n",
    "    \"\"\"\n",
    "    \n",
    "    for second in range(0, model_run_time, delta_t): # range(start,stop,step)\n",
    "        \n",
    "        # Step 1: Initialize Thermal Capacitance\n",
    "        gamma = thermal_capacitance(BasCap, area_of_facets, tau, rho)\n",
    "\n",
    "        \n",
    "        # Step 2: Update the surface layer (T_0)\n",
    "        T_0 = heat_transfer_array[:,:,0]\n",
    "        T_1 = heat_transfer_array[:,:,1]\n",
    "        W_t = flux_input_func(W_0, zenith_angle_array, albedo)\n",
    "        dQ0_dt = rate_of_change_surface_thermal_energy_func(W_t[:-1,:-1], area_of_facets, T_0, T_1, sigma, BasCond, tau)\n",
    "        T_0 += dQ0_dt * delta_t / gamma\n",
    "\n",
    "\n",
    "        # Step 3: Update intermediate layers\n",
    "        T_n = heat_transfer_array[:,:,1:-1]  # All intermediate layers\n",
    "        T_np1 = heat_transfer_array[:,:,2:]  # Next layer in each direction\n",
    "        T_nm1 = heat_transfer_array[:,:,:-2]  # Previous layer in each direction\n",
    "        dQn_dt = rate_of_change_intermediate_thermal_energy_func(T_np1, T_n, T_nm1, BasCond, area_of_facets[:,:,np.newaxis], tau)\n",
    "        T_n += dQn_dt * delta_t / gamma[:,:,np.newaxis]\n",
    "\n",
    "        # Step 4: Update the final layer (T_f)\n",
    "        T_f = heat_transfer_array[:,:,-1]\n",
    "        T_fm1 = heat_transfer_array[:,:,-2]\n",
    "        T_nm1 = T_nm1[:,:,-1]\n",
    "        dQf_dt = rate_of_change_final_thermal_energy_func(T_const, T_f, T_nm1, BasCond, area_of_facets, tau)\n",
    "        T_f += dQf_dt * delta_t / gamma\n",
    "\n",
    "        # Step 5: Update the heat transfer array with new temperatures\n",
    "        heat_transfer_array[:,:,0] = T_0\n",
    "        heat_transfer_array[:,:,1:-1] = T_n\n",
    "        heat_transfer_array[:,:,-1] = T_f\n",
    "\n",
    "        # uncomment for diagnosing\n",
    "        # print(f\"At second {second}: T_0 = {T_0}, T_1 = {T_1}, W_t = {W_t}, dQ0_dt = {dQ0_dt}\")\n",
    "            # Independent phase_step increment based on the condition\n",
    "        # if second % 6700 == 0:\n",
    "        #     phase_step += 1\n",
    "    return heat_transfer_array\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fa667af7-f0cc-4459-b521-273c0e25ab84",
   "metadata": {
    "tags": []
   },
   "source": [
    "## Create the heat transfer array and run function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "id": "d19fa635-9682-409d-92c5-e8fe247b3b95",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# we know the heat_input_array is of the shape (360,91) for which there are 360 phase degree steps.\n",
    "# The full phase model will be of shape (360,360,91) and now we want to add a depth component of 100 layers\n",
    "# This will bring the shape to (360,360,91,100) of which its components are (phase, longitude, latitude, depth)\n",
    "# of which there is one depth for component for each latitude, longitude pair per phase.\n",
    "# This turns out to be 1.1 billion elements, lets instead start with a per phase_step"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "id": "93120c42-61b8-45be-95e8-78916c44d19d",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# T_const = 225# The initial temperature to which all values are set to initially\n",
    "\n",
    "# # creating an array of (longitude, latitude, depth) filled with T_const\n",
    "# # Ensure the heat_transfer_array is float64 for floating-point operations\n",
    "# # must be run each time to reset the heat_transfer_array \n",
    "# heat_transfer_array = np.full((360, 360, 90, 100), T_const).astype(np.float16)\n",
    "\n",
    "# # setting various varibles\n",
    "# sigma = 5.67e-8  # Stefan-Boltzmann constant (W/m²K⁴)\n",
    "# W_0 = 1365.0   # the energy flux arriving from the Sun\n",
    "# BasCond = 1.50e3                 #Basalt conductivity in W.m-3.K-1\n",
    "# albedo = 0.15 # lunar albedo\n",
    "# BasCap = 670.0 #Basalt thermal capacity in J Kg-1.K-1 \n",
    "# rho = 1500 # Lunar regolith density in Kg.m-3\n",
    "# tau = 0.01 # tau is the thickness of the slice (meters)\n",
    "# delta_t = 500 # time over which each loop calculation is preformed\n",
    "# model_run_time = 500 # simulation run time in seconds\n",
    "\n",
    "# # run the function \n",
    "# heat_transfer_model_array = heat_transfer_model_func(heat_transfer_array, \n",
    "#                                                      W_0, \n",
    "#                                                      albedo, \n",
    "#                                                      T_const, \n",
    "#                                                      BasCap, \n",
    "#                                                      area_of_facets, \n",
    "#                                                      tau, \n",
    "#                                                      rho, \n",
    "#                                                      delta_t, \n",
    "#                                                      model_run_time)\n",
    "\n",
    "# heat_transfer_model_array.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "id": "3c03de60-2700-43e0-a28c-5fe740ae3750",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\seane\\AppData\\Local\\Temp\\ipykernel_26956\\2017439074.py:21: RuntimeWarning: overflow encountered in power\n",
      "  dQ0_dt = area_of_facets * (flux_input_array - sigma * T_0**4) + ((T_1 - T_0) * BasCond * area_of_facets/tau)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(360, 90, 100)"
      ]
     },
     "execution_count": 63,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "T_const = 225# The initial temperature to which all values are set to initially\n",
    "\n",
    "# creating an array of (longitude, latitude, depth) filled with T_const\n",
    "# Ensure the heat_transfer_array is float64 for floating-point operations\n",
    "# must be run each time to reset the heat_transfer_array \n",
    "heat_transfer_array = np.full((360, 90, 100), T_const).astype(np.float16)\n",
    "\n",
    "# setting various varibles\n",
    "sigma = 5.67e-8  # Stefan-Boltzmann constant (W/m²K⁴)\n",
    "W_0 = 1365.0   # the energy flux arriving from the Sun\n",
    "BasCond = 1.50e3                 #Basalt conductivity in W.m-3.K-1\n",
    "albedo = 0.15 # lunar albedo\n",
    "BasCap = 670.0 #Basalt thermal capacity in J Kg-1.K-1 \n",
    "rho = 1500 # Lunar regolith density in Kg.m-3\n",
    "tau = 0.01 # tau is the thickness of the slice (meters)\n",
    "delta_t = 500 # time over which each loop calculation is preformed\n",
    "model_run_time = 500 # simulation run time in seconds\n",
    "\n",
    "# run the function \n",
    "heat_transfer_model_array = heat_transfer_model_func(heat_transfer_array, \n",
    "                                                     W_0, \n",
    "                                                     albedo, \n",
    "                                                     T_const, \n",
    "                                                     BasCap, \n",
    "                                                     area_of_facets, \n",
    "                                                     tau, \n",
    "                                                     rho, \n",
    "                                                     delta_t, \n",
    "                                                     model_run_time)\n",
    "\n",
    "heat_transfer_model_array.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ac433ee0-5c54-4401-9903-ef38245697e0",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "# heat_transfer_model_array[0,180,0,:] # Sanity check, values where longitude = 0, theta = 0, and the full depth of 100 layers"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fb4b76c2-fcf7-4d0f-82bc-7be18307ac2f",
   "metadata": {
    "jp-MarkdownHeadingCollapsed": true,
    "tags": []
   },
   "source": [
    "## The Radio Brightness Temperature"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4c4ba02a-260d-44a8-a5e3-7aa84bd299dd",
   "metadata": {},
   "outputs": [],
   "source": [
    "def radio_brightness_temp_func(alpha, T_const, heat_transfer_model_array):\n",
    "    \"\"\"\n",
    "    Converting the temperature at each layer to a Temperature Brightness [Kelvin]\n",
    "    \n",
    "    Parameters:\n",
    "    - T_const: The heat at the base layer, does not change\n",
    "    - alpha: \n",
    "    \"\"\"\n",
    "    # act upon the bottom layer\n",
    "    Chi_bottom_layer = alpha * T_const + (1 - alpha) * heat_transfer_model_array[:,:,-1]\n",
    "    \n",
    "    # act upon the intermediate layers\n",
    "    Chi_middle_layers = alpha * heat_transfer_model_array"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "82e43396-f33a-4a04-becd-b2cbf7713246",
   "metadata": {},
   "outputs": [],
   "source": [
    "def radio_brightness_temp_func(alpha, T_const, heat_transfer_model_array):\n",
    "    \"\"\"\n",
    "    Converts the temperature at each layer to brightness temperature [Kelvin].\n",
    "\n",
    "    Parameters:\n",
    "    - alpha: Weighting factor for the brightness temperature formula.\n",
    "    - T_const: Constant temperature at the base layer (T*).\n",
    "    - heat_transfer_model_array: 3D array representing temperatures across layers, \n",
    "      with the last axis representing the depth (from bottom to top).\n",
    "\n",
    "    Returns:\n",
    "    - chi_0: The brightness temperature at the top layer radiated into space.\n",
    "    - chi_layers: The brightness temperature for all layers.\n",
    "    \"\"\"\n",
    "    \n",
    "    # Step 1: Get the number of layers (depth is the last axis)\n",
    "    num_layers = heat_transfer_model_array.shape[2]\n",
    "    \n",
    "    # Step 2: Initialize the chi array to hold brightness temperatures for each layer\n",
    "    chi_layers = np.zeros_like(heat_transfer_model_array)\n",
    "    \n",
    "    # Step 3: Calculate the brightness temperature for the bottom layer (last layer in the stack)\n",
    "    chi_layers[:, :, -1] = alpha * T_const + (1 - alpha) * heat_transfer_model_array[:, :, -1]\n",
    "    \n",
    "    # Step 4: From the bottom to the top of the stack\n",
    "    for n in range(num_layers - 2, -1, -1):\n",
    "        chi_layers[:, :, n] = alpha * chi_layers[:, :, n+1] + (1 - alpha) * heat_transfer_model_array[:, :, n]\n",
    "    \n",
    "    # Step 5: Return the brightness temperature at the top layer (chi_0) and all layers\n",
    "    chi_0 = chi_layers[:, :, 0]\n",
    "    return chi_0 # , chi_layers\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d0a55bcb-8c41-45d4-b09b-932cd8e0ef21",
   "metadata": {
    "jp-MarkdownHeadingCollapsed": true,
    "tags": []
   },
   "source": [
    "## The Radio Flux Density"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5c9924b0-cba6-4e1f-8b1d-b700a454d31c",
   "metadata": {},
   "source": [
    "If Ωj is the solid angle subtended by the jth facet as viewed from the Earth, the contribution to the radio flux density is:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0e581192-0393-426d-a249-9dace76dbde5",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Constants\n",
    "k_B = 1.38e-23  # Boltzmann constant in J/K\n",
    "moon_radius = 1.737e6  # Radius of the Moon in meters\n",
    "distance_to_earth = 3.844e8  # Average distance from the Moon to Earth in meters\n",
    "\n",
    "def radio_flux_density(chi_0, solid_angles, wavelength):\n",
    "    \"\"\"\n",
    "    Calculate the radio flux density contribution from each facet.\n",
    "\n",
    "    Parameters:\n",
    "    - chi_0: 2D array of brightness temperatures at the top layer for each facet (K).\n",
    "    - solid_angles: 2D array of solid angles subtended by each facet (sr).\n",
    "    - wavelength: Wavelength of the radio emission (in meters).\n",
    "\n",
    "    Returns:\n",
    "    - total_flux_density: The total radio flux density (W/m²/Hz).\n",
    "    - delta_S: 2D array of individual contributions to the flux density from each facet.\n",
    "    \"\"\"\n",
    "    \n",
    "    # Calculate the flux density contribution for each facet using the formula\n",
    "    delta_S = (2 * k_B * solid_angles * chi_0) / wavelength**2\n",
    "    \n",
    "    # Sum over all facets to get the total flux density\n",
    "    total_flux_density = np.sum(delta_S)\n",
    "    \n",
    "    return total_flux_density #, delta_S"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fe7ce764-3b5e-4578-b52f-d05471ced9b2",
   "metadata": {},
   "outputs": [],
   "source": [
    "def average_radio_brightness_temperature(total_flux_density, wavelength, solid_angle_moon):\n",
    "    \"\"\"\n",
    "    Calculate the average radio brightness temperature of the lunar disc.\n",
    "\n",
    "    Parameters:\n",
    "    - S: Total radio flux density (W/m²/Hz).\n",
    "    - wavelength: Wavelength of the radio emission (in meters).\n",
    "    - solid_angle_moon: Solid angle subtended by the Moon (in steradians).\n",
    "\n",
    "    Returns:\n",
    "    - T_b: The average radio brightness temperature of the lunar disc (K).\n",
    "    \"\"\"\n",
    "    \n",
    "    # Calculate the average brightness temperature\n",
    "    T_b = (total_flux_density * wavelength**2) / (2 * k_B * solid_angle_moon)\n",
    "    \n",
    "    return T_b"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7f297154-a0d8-4d4d-b824-f3ad7b44fca4",
   "metadata": {
    "tags": []
   },
   "source": [
    "## Solid Angle Calculations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "df12bac5-a7b9-4235-909e-302a943d94e5",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def generate_moon_facets_areas_vectorized(latitude_array, longitude_array, moon_radius):\n",
    "    \"\"\"\n",
    "    Generation of facet areas over the surface of the Moon using latitude and longitude arrays.\n",
    "    \n",
    "    facet_area = R2×∣sin(latitudei+1)−sin(latitudei)∣×dlongitude\n",
    "    \n",
    "    https://math.stackexchange.com/questions/4102850/area-of-surface-between-two-lines-of-latitude\n",
    "    \n",
    "    Parameters:\n",
    "    - latitude_array: Array of latitudes in degrees.\n",
    "    - longitude_array: Array of longitudes in degrees.\n",
    "    - moon_radius: Radius of the Moon in meters.\n",
    "    \n",
    "    Returns:\n",
    "    - facet_areas: 2D array of facet areas (m²) corresponding to each latitude and longitude grid cell.\n",
    "    \"\"\"\n",
    "    \n",
    "    # Convert latitude and longitude arrays from degrees to radians\n",
    "    latitudes = np.radians(latitude_array)\n",
    "    longitudes = np.radians(longitude_array)\n",
    "    \n",
    "    # Calculate differences in latitude and longitude\n",
    "    dlat = np.diff(latitudes)[:, np.newaxis]  # Latitude differences (shape: (len(latitudes)-1, 1))\n",
    "    dlon = np.diff(longitudes)[np.newaxis, :]  # Longitude differences (shape: (1, len(longitudes)-1))\n",
    "    \n",
    "    # Compute the areas of the facets\n",
    "    facet_areas = (moon_radius**2) * np.abs(np.sin(latitudes[1:, np.newaxis]) - np.sin(latitudes[:-1, np.newaxis])) * dlon\n",
    "    \n",
    "    return facet_areas\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b9337d52-be74-4553-b3d8-bbd86b257a3b",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def calculate_solid_angles_for_moon_with_projection(moon_radius, distance_to_earth, latitude_array, longitude_array):\n",
    "    \"\"\"\n",
    "    Vectorized calculation of the solid angles subtended by all facets on the Moon using latitude and longitude arrays,\n",
    "    accounting for projection effects.\n",
    "    \n",
    "    Parameters:\n",
    "    - moon_radius: Radius of the Moon (m).\n",
    "    - distance_to_earth: Distance from the Moon to Earth (m).\n",
    "    - latitude_array: Array of latitudes in degrees.\n",
    "    - longitude_array: Array of longitudes in degrees.\n",
    "\n",
    "    Returns:\n",
    "    - solid_angles: 2D array of solid angles for each facet (steradians).\n",
    "    \"\"\"\n",
    "    \n",
    "    # Generate facet areas using the vectorized function\n",
    "    facet_areas = generate_moon_facets_areas_vectorized(latitude_array, longitude_array, moon_radius)\n",
    "    \n",
    "    # Convert latitude and longitude arrays from degrees to radians\n",
    "    latitudes = np.radians(latitude_array)\n",
    "    longitudes = np.radians(longitude_array)\n",
    "    \n",
    "    # # Get the central latitude and longitude (assuming observer is directly above equator)\n",
    "    # observer_lat = 0.0  # Observer on Earth looking at the Moon's equator\n",
    "    # observer_lon = 0.0  # Assume the sub-Earth point is at longitude 0\n",
    "    \n",
    "    # Broadcast latitudes and longitudes into 2D grids\n",
    "    lat_grid, lon_grid = np.meshgrid(latitudes, longitudes, indexing='ij')\n",
    "    \n",
    "    # Calculate the unit normal vectors for each facet\n",
    "    x_facet = np.cos(lat_grid) * np.cos(lon_grid)\n",
    "    y_facet = np.cos(lat_grid) * np.sin(lon_grid)\n",
    "    z_facet = np.sin(lat_grid)\n",
    "    \n",
    "    # Facet normal vectors (shape: [len(latitudes), len(longitudes), 3])\n",
    "    facet_normals = np.stack([x_facet, y_facet, z_facet], axis=-1)\n",
    "    \n",
    "    # Observer direction vector (assumed to be [0, 0, 1] in this case)\n",
    "    observer_vector = np.array([0, 0, 1])\n",
    "    \n",
    "    # Calculate the cosine of the angle between facet normals and observer vector\n",
    "    cos_theta = np.dot(facet_normals, observer_vector)  # Dot product gives cos(theta)\n",
    "    \n",
    "    # Ensure no negative cosines (only facets facing Earth contribute)\n",
    "    cos_theta = np.clip(cos_theta, 0, 1)  # Set cos(theta) to 0 for facets not facing Earth\n",
    "    \n",
    "    # Calculate solid angles, including the projection factor\n",
    "    solid_angles = (facet_areas * cos_theta[:-1, :-1]) / distance_to_earth**2\n",
    "    # solid_angles = (facet_areas * cos_theta) / distance_to_earth**2\n",
    "\n",
    "    \n",
    "    return solid_angles\n",
    "\n",
    "\n",
    "# Latitude and longitude arrays \n",
    "latitude_array = np.arange(0, 91)  # degrees\n",
    "longitude_array = np.arange(-180, 180)  # degrees\n",
    "\n",
    "# Calculate solid angles\n",
    "solid_angles_with_projection = calculate_solid_angles_for_moon_with_projection(moon_radius, distance_to_earth, latitude_array, longitude_array)\n",
    "\n",
    "solid_angles_with_projection.shape  # Display the vectorized solid angles"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fd8bf8c0-8787-4a48-b5d5-11168d6a0373",
   "metadata": {},
   "outputs": [],
   "source": [
    "solid_angles_with_projection"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dbf95223-c664-4a9b-acad-08182bd22f0e",
   "metadata": {},
   "source": [
    "### Kens moon_temp_model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "aad2b61b-6206-4d15-aec8-c3993ba47b98",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "# Model Constants\n",
    "DELTA_PHASE = 0.05                # Phase increment in degrees\n",
    "LUNATION_COUNT = 16               # Number of lunations\n",
    "LUNATION_SECONDS = ((29 * 24 + 12) * 60 + 44) * 60\n",
    "\n",
    "# File setup for output\n",
    "fname = \"moontemp_phaseloop.csv\"\n",
    "\n",
    "# Physical Constants\n",
    "PI = np.pi\n",
    "SIGMA = 5.67e-8                # Stefan-Boltzmann Constant (W/m^2K^4)\n",
    "BAS_COND = 1.50e3              # Basalt conductivity (W.m-3.K-1)\n",
    "BAS_CAP = 670.0                # Basalt thermal capacity (J Kg-1.K-1)  \n",
    "RHO = 1500.0                   # Lunar regolith density (Kg.m-3)\n",
    "ALBEDO = 0.15                  # Lunar albedo (dimensionless)\n",
    "SOLAR_CONSTANT = 1365.0        # Solar irradiance (W/m^2)\n",
    "T_MEAN = 225.0                 # Mean lunar soil temperature in Kelvins\n",
    "RF_ABSORPTION_COEFF = 1        # RF Transmission Coefficient\n",
    "\n",
    "# Model Parameters\n",
    "A = 1.0                        # Surface area of element (square meters)\n",
    "DELTA_Z = 0.01                 # Thickness of element (meters)\n",
    "MASS_ELEMENT = A * DELTA_Z * RHO  # Mass of element (kg)\n",
    "\n",
    "# Derived constants\n",
    "I0 = (1 - ALBEDO) * SOLAR_CONSTANT  # Solar energy absorbed by lunar surface\n",
    "BETA = BAS_COND * A / DELTA_Z       # Conductivity parameter\n",
    "GAMMA = BAS_CAP * MASS_ELEMENT      # Thermal capacity parameter\n",
    "ALPHA = np.exp(-RF_ABSORPTION_COEFF * DELTA_Z)  # RF Transmission Coefficient\n",
    "\n",
    "def calculate_zenith_angle(phase):\n",
    "    \"\"\"Calculate the zenith angle of the sun based on the phase.\"\"\"\n",
    "    zenith_angle = abs(phase - 180)\n",
    "    return zenith_angle\n",
    "\n",
    "def calculate_solar_irradiance(phase):\n",
    "    \"\"\"Calculate solar irradiance at the surface based on phase.\"\"\"\n",
    "    zenith_angle = calculate_zenith_angle(phase)\n",
    "    if zenith_angle > 90:\n",
    "        return 0\n",
    "    irradiance = I0 * np.cos(np.radians(zenith_angle))\n",
    "    return irradiance\n",
    "\n",
    "def update_temperature_intermediate_layers(LTK, irradiance, delta_t):\n",
    "    \"\"\"Update the temperature of intermediate soil layers.\"\"\"\n",
    "    for j in range(29):\n",
    "        if j == 0:\n",
    "            # Surface layer with solar irradiance and radiation loss\n",
    "            LTK[j] += (A * (irradiance - SIGMA * LTK[j] ** 4) + (LTK[j + 1] - LTK[j]) * BETA) * delta_t / GAMMA\n",
    "        else:\n",
    "            # Intermediate layers with only conduction\n",
    "            LTK[j] += ((LTK[j + 1] - LTK[j]) * BETA + (LTK[j - 1] - LTK[j]) * BETA) * delta_t / GAMMA\n",
    "    return LTK\n",
    "\n",
    "def calculate_radio_brightness_temperature(LTK):\n",
    "    \"\"\"Calculate the radio brightness temperature based on the soil temperature profile.\"\"\"\n",
    "    Tb = T_MEAN  # Initial brightness temperature\n",
    "    for k in range(29, -1, -1):\n",
    "        Tb = ALPHA * Tb + (1 - ALPHA) * LTK[k]\n",
    "    LTK[30] = Tb\n",
    "    return LTK\n",
    "\n",
    "def soilstack(LTK, phase, delta_phase):\n",
    "    \"\"\"Simulate the lunar soil temperature and update brightness temperature.\"\"\"\n",
    "    delta_t = delta_phase * LUNATION_SECONDS / 360\n",
    "    irradiance = calculate_solar_irradiance(phase)\n",
    "    LTK = update_temperature_intermediate_layers(LTK, irradiance, delta_t)\n",
    "    LTK = calculate_radio_brightness_temperature(LTK)\n",
    "    return LTK\n",
    "\n",
    "# Initialize soil temperatures to T_mean (30 soil temperatures, LTK[0] to LTK[29])\n",
    "LTK = [T_MEAN] * 31\n",
    "\n",
    "# Run the lunation and phase loop and write results to a CSV file\n",
    "with open(fname, 'w') as opfile:\n",
    "    for lunation in range(LUNATION_COUNT):\n",
    "        phase = 0.0\n",
    "        while phase <= 360:\n",
    "            # Run temperature update for each phase step\n",
    "            print(f\"Running soilstack for lunation {lunation}, phase {phase} degrees\")\n",
    "            LTK = soilstack(LTK, phase, DELTA_PHASE)\n",
    "            phase += DELTA_PHASE\n",
    "\n",
    "            # Format temperature values for output\n",
    "            temp_values = [f\"{0.1 * int(10 * temp):.1f}\" for temp in LTK]\n",
    "            opstring = f\"{lunation}, {phase:.1f}, \" + \", \".join(temp_values)\n",
    "\n",
    "            # Only output if lunation > 14\n",
    "            if lunation > 14:\n",
    "                print(opstring)\n",
    "                opfile.write(opstring + \"\\n\")\n",
    "\n",
    "print(\"Simulation complete\")\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}