build_Moon_scatter.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "19197a16",
   "metadata": {},
   "source": [
    "# Build a Scattering Model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3a725114",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "e04c0d7f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# The code for setting up the 3D models is based on this code: \n",
    "# https://github.com/williameaton/axisem3d-shapes/src\n",
    "# I made some modifications to the code (modified files have _latlon in the filename.)\n",
    "# Note that I have not tested the code for other geometries\n",
    "\n",
    "import numpy as np\n",
    "import sys\n",
    "import os\n",
    "import netCDF4 as nc4\n",
    "\n",
    "sys.path.append('./shapes_latlon')\n",
    "from gen_scripts_latlon import latlon_to_cartesian\n",
    "\n",
    "MOON_RADIUS_IN_KM = 1737.1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "ad4ad1b8",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "initiating model_latlon\n",
      "self.min_wavelength 4000.0\n",
      "x, y, z  5 300 1350\n",
      "x, y, z  5 301 1351\n",
      "1717100.0 1737100.0 20000.0\n",
      "[0.0, 0.125, 0.25, 0.375, 0.5]\n",
      "depth list  [20000.0, 15000.0, 10000.0, 5000.0, 0.0]\n",
      "grid radius  [1717100. 1722100. 1727100. 1732100. 1737100.]\n",
      "\n",
      "rho\n",
      "5 301 1351\n",
      "Spacing  4000.0 0.132890365448505 0.13323464100666174\n",
      "2033255\n",
      "rho - top layer min max\n",
      "-0.49999778298821107 0.4999949685397729\n",
      "rho - bottom layer min max\n",
      "0.0 0.0\n",
      "shape  (5, 301, 1351)\n",
      "initiated\n",
      "Writing spherical model...\n",
      "PARAVIEW FLAG = FALSE: MODEL OKAY FOR SIMULATION\n",
      "Data written to file  3D_models_output/linear_20.0.nc\n",
      "# ___________________________________________________________________________________\n",
      "    - 3D_models_output/linear_20.0:\n",
      "        activated: true\n",
      "        class_name: StructuredGridV3D\n",
      "        nc_data_file: 3D_models_output/linear_20.0.nc\n",
      "        coordinates:\n",
      "            horizontal: LATITUDE_LONGITUDE\n",
      "            vertical: RADIUS\n",
      "            ellipticity: FILL THIS IN - true/false\n",
      "            depth_below_solid_surface: FILL THIS IN - true/false\n",
      "            nc_variables: [radius, lon, lat]\n",
      "            data_rank: [2, 1, 0]\n",
      "            length_unit: m\n",
      "            angle_unit: degree\n",
      "            undulated_geometry: false\n",
      "            whole_element_inplane: false\n",
      "        properties:\n",
      "            - VP:\n",
      "                nc_var: vp\n",
      "                factor: 1\n",
      "                reference_kind: FILL THIS IN - ABS/REF1D/REF3D\n",
      "            - VS:\n",
      "                nc_var: vs\n",
      "                factor: 1\n",
      "                reference_kind: FILL THIS IN - ABS/REF1D/REF3D\n",
      "            - RHO:\n",
      "                nc_var: rho\n",
      "                factor: 1\n",
      "                reference_kind: FILL THIS IN - ABS/REF1D/REF3D\n",
      "        store_grid_only_on_leaders: true\n",
      "# ___________________________________________________________________________________\n"
     ]
    }
   ],
   "source": [
    "# Build a model with perturbation at each grid point\n",
    "# It's a ring running along longitude = 0.\n",
    "\n",
    "linear_inhomogeneity_depth = 20.\n",
    "\n",
    "# We dont need to make a model that spans the whole domain, just a segment near longitude = 0 \n",
    "radius = MOON_RADIUS_IN_KM*1000\n",
    "lat_lim = [-90, 90]\n",
    "long_lim = [-20, 20]\n",
    "radius_lim = [(MOON_RADIUS_IN_KM-linear_inhomogeneity_depth)*1000, MOON_RADIUS_IN_KM*1000]\n",
    "\n",
    "# If type=Spherical, the convention is that x = radius (out of the page), y = longitude, z = latitude, 0,0,0 at the centre of the body\n",
    " \n",
    "\n",
    "# Create our global model, including inhomogeneity\n",
    "glob_m = Model(\"spherical\", x_lim=radius_lim, y_lim=long_lim, z_lim=lat_lim, elements_per_wavelength=1, \n",
    "               dominant_freq=0.5, min_velocity=2000, oversaturation=1, a=radius, \n",
    "               linear_inhomogeneity_depth=linear_inhomogeneity_depth)\n",
    "#         If type=Spherical, the convention is that x = radius (out of the page), y = longitude, z = latitude, 0,0,0 at the centre of the body\n",
    "\n",
    "#\n",
    "os.makedirs('3D_models_output', exist_ok=True)\n",
    "# Write to NetCDF file \n",
    "glob_m.writeNetCDF(\"3D_models_output/linear_{}.nc\".format(linear_inhomogeneity_depth))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6b329e2b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# # Build a model with perturbation at each grid point\n",
    "# # It's a ring running along longitude = 0.\n",
    "\n",
    "# linear_inhomogeneity_depth = 50.\n",
    "\n",
    "# # We dont need to make a model that spans the whole domain, just a segment near longitude = 0 =L\n",
    "# radius = MOON_RADIUS_IN_KM*1000\n",
    "# lat_lim = [-90, 90]\n",
    "# long_lim = [-20, 20]\n",
    "# radius_lim = [(MOON_RADIUS_IN_KM-linear_inhomogeneity_depth)*1000, MOON_RADIUS_IN_KM*1000]\n",
    "\n",
    "# # If type=Spherical, the convention is that x = radius (out of the page), y = longitude, z = latitude, 0,0,0 at the centre of the body\n",
    "\n",
    "\n",
    "# # Create our global model, including inhomogeneity\n",
    "# glob_m = Model(\"spherical\", x_lim=radius_lim, y_lim=long_lim, z_lim=lat_lim, elements_per_wavelength=1, \n",
    "#                dominant_freq=0.5, min_velocity=2200, oversaturation=1, a=radius, \n",
    "#                linear_inhomogeneity_depth=linear_inhomogeneity_depth)\n",
    "# #         If type=Spherical, the convention is that x = radius (out of the page), y = longitude, z = latitude, 0,0,0 at the centre of the body\n",
    "\n",
    "# # Write to NetCDF file \n",
    "# glob_m.writeNetCDF(\"3D_models_output/linear_{}.nc\".format(linear_inhomogeneity_depth))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4c6f12ec",
   "metadata": {},
   "outputs": [],
   "source": [
    "# # Build a model with perturbation at each grid point\n",
    "# # It's a ring running along longitude = 0.\n",
    "\n",
    "# linear_inhomogeneity_depth = 80.\n",
    "\n",
    "# # We dont need to make a model that spans the whole domain, just a segment near longitude = 0 \n",
    "# radius = MOON_RADIUS_IN_KM*1000\n",
    "# lat_lim = [-90, 90]\n",
    "# long_lim = [-20, 20]\n",
    "# radius_lim = [(MOON_RADIUS_IN_KM-linear_inhomogeneity_depth)*1000, MOON_RADIUS_IN_KM*1000]\n",
    "\n",
    "# # If type=Spherical, the convention is that x = radius (out of the page), y = longitude, z = latitude, 0,0,0 at the centre of the body\n",
    "\n",
    "\n",
    "# # Create our global model, including inhomogeneity\n",
    "# glob_m = Model(\"spherical\", x_lim=radius_lim, y_lim=long_lim, z_lim=lat_lim, elements_per_wavelength=1, \n",
    "#                dominant_freq=0.5, min_velocity=2200, oversaturation=1, a=radius, \n",
    "#                linear_inhomogeneity_depth=linear_inhomogeneity_depth)\n",
    "# #         If type=Spherical, the convention is that x = radius (out of the page), y = longitude, z = latitude, 0,0,0 at the centre of the body\n",
    "\n",
    "# # Write to NetCDF file \n",
    "# glob_m.writeNetCDF(\"3D_models_output/linear_{}.nc\".format(linear_inhomogeneity_depth))\n",
    "\n",
    "# # approx 4.4 km spacing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d13545a7",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "74283b90",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Build a model with perturbation at each grid point\n",
    "# Test a whole Moon model\n",
    "\n",
    "linear_inhomogeneity_depth = 20.\n",
    "\n",
    "# Span the whole domain\n",
    "radius = MOON_RADIUS_IN_KM*1000\n",
    "lat_lim = [-90, 90]\n",
    "long_lim = [-180, 180]\n",
    "radius_lim = [(MOON_RADIUS_IN_KM-linear_inhomogeneity_depth)*1000, MOON_RADIUS_IN_KM*1000]\n",
    "\n",
    "# If type=Spherical, the convention is that x = radius (out of the page), y = longitude, z = latitude, 0,0,0 at the centre of the body\n",
    "\n",
    "\n",
    "# Create our global model, including inhomogeneity\n",
    "glob_m = Model(\"spherical\", x_lim=radius_lim, y_lim=long_lim, z_lim=lat_lim, elements_per_wavelength=1, \n",
    "               dominant_freq=0.5, min_velocity=2000, oversaturation=1, a=radius, \n",
    "               linear_inhomogeneity_depth=linear_inhomogeneity_depth)\n",
    "#         If type=Spherical, the convention is that x = radius (out of the page), y = longitude, z = latitude, 0,0,0 at the centre of the body\n",
    "\n",
    "# Write to NetCDF file \n",
    "# This makes a large file, so is currently commented out\n",
    "# glob_m.writeNetCDF(\"3D_models_output/linear_{}_whole.nc\".format(linear_inhomogeneity_depth))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ff048e18",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Build a model with perturbation at each grid point\n",
    "# Test a whole Moon model\n",
    "\n",
    "linear_inhomogeneity_depth = 50.\n",
    "\n",
    "# Span the whole domain\n",
    "radius = MOON_RADIUS_IN_KM*1000\n",
    "lat_lim = [-90, 90]\n",
    "long_lim = [-180, 180]\n",
    "radius_lim = [(MOON_RADIUS_IN_KM-linear_inhomogeneity_depth)*1000, MOON_RADIUS_IN_KM*1000]\n",
    "\n",
    "# If type=Spherical, the convention is that x = radius (out of the page), y = longitude, z = latitude, 0,0,0 at the centre of the body\n",
    "\n",
    "\n",
    "# Create our global model, including inhomogeneity\n",
    "glob_m = Model(\"spherical\", x_lim=radius_lim, y_lim=long_lim, z_lim=lat_lim, elements_per_wavelength=1, \n",
    "               dominant_freq=0.5, min_velocity=2000, oversaturation=1, a=radius, \n",
    "               linear_inhomogeneity_depth=linear_inhomogeneity_depth)\n",
    "#         If type=Spherical, the convention is that x = radius (out of the page), y = longitude, z = latitude, 0,0,0 at the centre of the body\n",
    "\n",
    "# Write to NetCDF file \n",
    "# This makes a large file, so is currently commented out\n",
    "# glob_m.writeNetCDF(\"3D_models_output/linear_{}_whole.nc\".format(linear_inhomogeneity_depth))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "81d2b3b8-1578-4235-86de-25c05f8a96b6",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "initiating model_latlon\n",
      "self.min_wavelength 4000.0\n",
      "x, y, z  20 2700 1350\n",
      "x, y, z  21 2701 1351\n",
      "1657100.0 1737100.0 80000.0\n",
      "[0.0, 0.025000000000000022, 0.04999999999999999, 0.07500000000000001, 0.09999999999999998, 0.125, 0.15000000000000002, 0.175, 0.2, 0.22499999999999998, 0.25, 0.275, 0.3, 0.325, 0.35, 0.375, 0.4, 0.425, 0.45, 0.475, 0.5]\n",
      "depth list  [80000.0, 76000.0, 72000.0, 68000.0, 64000.0, 60000.0, 56000.0, 52000.0, 48000.0, 44000.0, 40000.0, 36000.0, 32000.0, 28000.0, 24000.0, 20000.0, 16000.0, 12000.0, 8000.0, 4000.0, 0.0]\n",
      "grid radius  [1657100. 1661100. 1665100. 1669100. 1673100. 1677100. 1681100. 1685100.\n",
      " 1689100. 1693100. 1697100. 1701100. 1705100. 1709100. 1713100. 1717100.\n",
      " 1721100. 1725100. 1729100. 1733100. 1737100.]\n",
      "\n",
      "rho\n",
      "21 2701 1351\n",
      "Spacing  3809.5238095238096 0.13328396890040725 0.13323464100666174\n",
      "76630071\n",
      "rho - top layer min max\n",
      "-0.49999996605542407 0.4999999914689215\n",
      "rho - bottom layer min max\n",
      "-0.0 -0.0\n",
      "shape  (21, 2701, 1351)\n",
      "initiated\n",
      "Writing spherical model...\n",
      "PARAVIEW FLAG = FALSE: MODEL OKAY FOR SIMULATION\n",
      "Data written to file  3D_models_output/linear_80.0_whole.nc\n",
      "# ___________________________________________________________________________________\n",
      "    - 3D_models_output/linear_80.0_whole:\n",
      "        activated: true\n",
      "        class_name: StructuredGridV3D\n",
      "        nc_data_file: 3D_models_output/linear_80.0_whole.nc\n",
      "        coordinates:\n",
      "            horizontal: LATITUDE_LONGITUDE\n",
      "            vertical: RADIUS\n",
      "            ellipticity: FILL THIS IN - true/false\n",
      "            depth_below_solid_surface: FILL THIS IN - true/false\n",
      "            nc_variables: [radius, lon, lat]\n",
      "            data_rank: [2, 1, 0]\n",
      "            length_unit: m\n",
      "            angle_unit: degree\n",
      "            undulated_geometry: false\n",
      "            whole_element_inplane: false\n",
      "        properties:\n",
      "            - VP:\n",
      "                nc_var: vp\n",
      "                factor: 1\n",
      "                reference_kind: FILL THIS IN - ABS/REF1D/REF3D\n",
      "            - VS:\n",
      "                nc_var: vs\n",
      "                factor: 1\n",
      "                reference_kind: FILL THIS IN - ABS/REF1D/REF3D\n",
      "            - RHO:\n",
      "                nc_var: rho\n",
      "                factor: 1\n",
      "                reference_kind: FILL THIS IN - ABS/REF1D/REF3D\n",
      "        store_grid_only_on_leaders: true\n",
      "# ___________________________________________________________________________________\n"
     ]
    }
   ],
   "source": [
    "# Build a model with perturbation at each grid point\n",
    "# Test a whole Moon model\n",
    "\n",
    "linear_inhomogeneity_depth = 80.\n",
    "\n",
    "# Span the whole domain\n",
    "radius = MOON_RADIUS_IN_KM*1000\n",
    "lat_lim = [-90, 90]\n",
    "long_lim = [-180, 180]\n",
    "radius_lim = [(MOON_RADIUS_IN_KM-linear_inhomogeneity_depth)*1000, MOON_RADIUS_IN_KM*1000]\n",
    "\n",
    "# If type=Spherical, the convention is that x = radius (out of the page), y = longitude, z = latitude, 0,0,0 at the centre of the body\n",
    "\n",
    "\n",
    "# Create our global model, including inhomogeneity\n",
    "glob_m = Model(\"spherical\", x_lim=radius_lim, y_lim=long_lim, z_lim=lat_lim, elements_per_wavelength=1, \n",
    "               dominant_freq=0.5, min_velocity=2000, oversaturation=1, a=radius, \n",
    "               linear_inhomogeneity_depth=linear_inhomogeneity_depth)\n",
    "#         If type=Spherical, the convention is that x = radius (out of the page), y = longitude, z = latitude, 0,0,0 at the centre of the body\n",
    "\n",
    "# Write to NetCDF file \n",
    "# This makes a large file, so is currently commented out\n",
    "# glob_m.writeNetCDF(\"3D_models_output/linear_{}_whole.nc\".format(linear_inhomogeneity_depth))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "71d8f98c-569d-4bd9-9757-0512191007f6",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "lat  -90.0 90.0\n",
      "lon  -180.0 180.0\n",
      "radius  1657100.0 1737100.0\n",
      "36494248 / 76630071 47.62392559965134 %\n",
      "rho shape  (21, 2701, 1351)\n",
      "rho - top layer min max\n",
      "-0.49999997 0.5\n",
      "rho - bottom layer min max\n",
      "-0.0 -0.0\n"
     ]
    }
   ],
   "source": [
    "# # view the file just made - currently commented out\n",
    "# ds = nc4.Dataset('3D_models_output/linear_80.0_whole.nc')\n",
    "# lat = ds.variables['lat'][:]\n",
    "# print('lat ',lat.min(),lat.max())\n",
    "# lon = ds.variables['lon'][:]\n",
    "# print('lon ',lon.min(),lon.max())\n",
    "# radius = ds.variables['radius'][:]\n",
    "# print('radius ',radius.min(),radius.max())\n",
    "\n",
    "# rho = ds.variables['rho'][:]\n",
    "# rho_set = np.count_nonzero((np.where(rho < 0, True, False)))\n",
    "# print(rho_set, '/', rho.size, (rho_set/rho.size*100), '%')\n",
    "\n",
    "# print('rho shape ', rho.shape)\n",
    "# print('rho - top layer min max')\n",
    "# print(rho[-1,:,:].min(),rho[-1,:,:].max())\n",
    "# print('rho - bottom layer min max')\n",
    "# print(rho[0,:,:].min(),rho[0,:,:].max())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5a78f2ac-247f-4331-9803-7bc6bffa3b7f",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.8.20"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}