SCHA.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<p style='text-align: justify;'>After the $n_k(m)$ solved, the potential field can be differentiated to obtain the field on different components. This code takes an input set of X, Y and Z magnetic field components and from them calculates a spherical cap harmonic model of the magnetic field within the region of a half degree radius cap around a point specified by the user. The declination, inclination and intensity of the model field are then plotted on a map. The spherical cap harmonic calculation is much more complicated than spherical harmonic, with each legendre polynomial being a summation, and n being non-integer.</p>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 120,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pickle\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import math\n",
    "import os\n",
    "import os.path"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<p style='text-align: justify;'>Define function to calculate distance on a sphere.</p>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 121,
   "metadata": {},
   "outputs": [],
   "source": [
    "def distance_on_unit_sphere(lat1, long1, lat2, long2):\n",
    "    # Convert latitude and longitude to spherical coordinates in radians.\n",
    "    degrees_to_radians = math.pi/180.0\n",
    "    \n",
    "    # phi = 90 - latitude\n",
    "    phi1 = (90.0 - lat1)*degrees_to_radians\n",
    "    phi2 = (90.0 - lat2)*degrees_to_radians\n",
    "    \n",
    "    # theta = longitude\n",
    "    theta1 = long1*degrees_to_radians\n",
    "    theta2 = long2*degrees_to_radians\n",
    "    \n",
    "    # Compute spherical distance from spherical coordinates.\n",
    "    cos = (math.sin(phi1)*math.sin(phi2)*math.cos(theta1 - theta2) + math.cos(phi1)*math.cos(phi2))\n",
    "    arc = math.acos( cos )\n",
    "    \n",
    "    return arc"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<p style='text-align: justify;'>The geomagnetic measurements data positions and components on the Earth's surface are usually given in a geographic coordinate system. Therefore, it is necessary to transform both coordinate locations and the vector-pointing directions to the spherical cap coordinate system to simplify calculations. From the spherical cap coordinate system, the transformed longitude and colatitude for any data site can be found by considering a triangle on the Earth’s surface with the centre of the spherical cap, the North geographic pole, and the data points as shown in following figure.</p>\n",
    "\n",
    "<img src=\"img/transformation.png\" width=\"400\"/>\n",
    "\n",
    "<p style='text-align: justify;'>Consider an arbitrary point $P$ in the geographic coordinate system which has coordinate $(\\varphi_P,\\theta_P)$ and the spherical cap pole with coordinate $(\\varphi_C,\\theta_C)$. The meridian of these points $(\\theta_C)$ and $(\\theta_P)$ pass through the centre of the cap and data sites, and simply become the sides of the spherical triangle with longitude difference between these two sides become the angle $(\\varphi_C-\\varphi_P)$. Using the spherical version of the cosine rule, the transformed colatitude $(\\theta_T)$ of the data point, i.e. relative to the centre of the cap, can be determined as follows \\begin{align}\\theta_T=\\cos^{-1⁡}(\\cos\\,⁡\\theta_C\\;\\cos\\,⁡\\theta_P+\\sin⁡\\,\\theta_C\\;\\sin\\,\\theta_P\\;\\cos\\,⁡(\\varphi_C-\\varphi_P))\\end{align}</p>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 122,
   "metadata": {},
   "outputs": [],
   "source": [
    "def sph_cosine_side(a, b, C):\n",
    "    # Input and output in radian\n",
    "    c = np.arccos(np.cos(a)*np.cos(b)+np.sin(a)*np.sin(b)*np.cos(C))\n",
    "    return c"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<p style='text-align: justify;'>While the transformed longitude $(\\varphi_T)$ can be determined using the spherical version of the cosine rule as follows.</p>\n",
    "\\begin{align}\\varphi_T=\\cos^{-1} \\left (\\frac{\\cos\\,\\theta_P-\\cos⁡\\,\\theta_T\\;\\cos⁡\\,\\theta_C}{\\sin⁡\\,\\theta_T\\;\\sin\\,\\theta_C} \\right)\\end{align}\n",
    "<p></p>\n",
    "<p style='text-align: justify;'>Similarly, the geomagnetic data at point $P$ also need to be transformed into the cap coordinate system. The rotation angle $\\gamma$ needs to be determined for every data point using the spherical version of the cosine rule.</p>\n",
    "\\begin{align}\\gamma=\\cos^{-1} \\left (\\frac{\\cos\\,\\theta_C-\\cos⁡\\,\\theta_T\\;\\cos⁡\\,\\theta_P}{\\sin⁡\\,\\theta_T\\;\\sin\\,\\theta_P} \\right)\\end{align}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 123,
   "metadata": {},
   "outputs": [],
   "source": [
    "def sph_cosine_angle(a, b, c):\n",
    "    C = np.arccos((np.cos(c)-np.cos(a)*np.cos(b))/(np.sin(a)*np.sin(b)))\n",
    "    return C"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<p style='text-align: justify;'>The Schmidt normalised Associated Legendre Polynomials in the solution above can be determined using</p>\n",
    "\n",
    "\\begin{align}P_n^m (\\cos\\,\\theta)=\\sum_{k=0}^K A_k (m,n) \\; {\\sin^{2k} \\, (\\theta/2)} \\end{align}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 124,
   "metadata": {},
   "outputs": [],
   "source": [
    "def assoc_Leg_func(n, A, thetaT):\n",
    "    P = np.zeros(len(n))\n",
    "    for i in range(len(n)):\n",
    "        for j in range(np.size(A,1)):\n",
    "            P[i] = P[i] + A[i,j]*(np.sin(thetaT/2)**(2*j))\n",
    "    return P"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<p style='text-align: justify;'>while the derivation of Schmidt normalised Associated Legendre Polynomials is</p>\n",
    "\n",
    "\\begin{align}\\frac{\\partial P_n^m (\\cos\\,\\theta)}{\\partial \\theta} = \\frac{\\sin\\,\\theta}{2} \\sum_{k=1}^K kA_k (m,n)\\;\\sin^{2(k-1)} ⁡(\\theta/2) + \\cos\\,\\theta\\; \\left [\\frac{m}{\\sin\\,\\theta} \\; P_n^m (\\cos⁡\\,\\theta) \\right]\\end{align}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 125,
   "metadata": {},
   "outputs": [],
   "source": [
    "def deriv_assoc_Leg_func(n, m, A, thetaT, P):\n",
    "    sigma = np.zeros(len(n))\n",
    "    dP = np.zeros(len(n))\n",
    "    for i in range(len(n)):\n",
    "        for j in range(1,np.size(A,1)):\n",
    "            sigma[i] = sigma[i] + j*A[i,j]*(np.sin(thetaT/2)**(2*(j-1)))\n",
    "        dP[i] = np.sin(thetaT)/2*sigma[i] + np.cos(thetaT)*(m[i]/np.sin(thetaT)*P[i])\n",
    "    return dP"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<p style='text-align: justify;'>The eastern component of geomagnetic field can be found as follows.</p>\n",
    "\n",
    "\\begin{align}B_x=\\frac{1}{r} \\frac{\\partial U(r,\\theta,\\varphi)} {\\partial \\theta}\\end{align}\n",
    "\n",
    "\\begin{align}=\\sum_{m=0}^\\infty \\, \\sum_{k=m}^\\infty \\left (\\frac{a}{r} \\right)^{n_k(m)+2} \\cdot \\frac{\\partial P_{n_k(m)}^m \\; (\\cos⁡ \\,\\theta_T)}{\\partial \\theta_T} \\; \\{g_k^m  \\cos⁡ \\,m\\varphi_T + h_k^m  \\sin⁡ \\, m\\varphi_T\\}\\end{align}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 126,
   "metadata": {},
   "outputs": [],
   "source": [
    "def Bx(m, dP, phiT, thetaT):\n",
    "    X_coeff=[]\n",
    "    for r in range(len(m)):\n",
    "        if m[r]==0: # m=0, we only compute g\n",
    "            X_coeff = X_coeff + [(np.sin(thetaT)**(m[r]))*dP[r]*np.cos(m[r]*phiT)] # g\n",
    "        else: # Otherwise, we need to compute g and h\n",
    "            X_coeff = X_coeff + [(np.sin(thetaT)**(m[r]))*dP[r]*np.cos(m[r]*phiT)] # g\n",
    "            X_coeff = X_coeff + [(np.sin(thetaT)**(m[r]))*dP[r]*np.sin(m[r]*phiT)] # h\n",
    "    return X_coeff"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<p style='text-align: justify;'>The northern component of geomagnetic field can be found as follows.</p>\n",
    "\n",
    "\\begin{align}B_y = - \\frac{1}{r \\sin\\,\\theta} \\; \\frac{\\partial U(r,\\theta,\\varphi)} {\\partial \\varphi}\\end{align}\n",
    "\n",
    "\\begin{align}=\\sum_{m=0}^\\infty \\sum_{k=m}^\\infty m \\left(\\frac{a}{r} \\right)^{n_k (m)+2} \\cdot \\frac{P_{n_k(m)}^m \\; (\\cos⁡\\,\\theta_T)} {\\sin⁡\\,\\theta_T} \\; \\{g_k^m \\sin⁡ \\,m\\varphi_T - h_k^m \\cos⁡ \\,m\\varphi_T\\}\\end{align}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 127,
   "metadata": {},
   "outputs": [],
   "source": [
    "def By(n, m, P, phiT, thetaT):\n",
    "    Y_coeff=[]\n",
    "    for r in range(len(m)):\n",
    "        if m[r]==0: # m=0, we only compute g\n",
    "            Y_coeff = Y_coeff + [(np.sin(thetaT)**(m[r]))*m[r]/np.sin(thetaT)*P[r]*np.sin(m[r]*phiT)] # g\n",
    "        else: # Otherwise, we need to compute g and h\n",
    "            Y_coeff = Y_coeff + [(np.sin(thetaT)**(m[r]))*m[r]/np.sin(thetaT)*P[r]*np.sin(m[r]*phiT)] # g\n",
    "            Y_coeff = Y_coeff + [(np.sin(thetaT)**(m[r]))*(-m[r])/np.sin(thetaT)*P[r]*np.cos(m[r]*phiT)] # g\n",
    "    return Y_coeff"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<p style='text-align: justify;'>The vertical component of geomagnetic field can be found as follows.</p>\n",
    "\n",
    "\\begin{align}B_z=\\frac{\\partial U(r,\\theta,\\varphi)}{\\partial r}\\end{align}\n",
    "\n",
    "\\begin{align}=-\\sum_{m=0}^\\infty \\sum_{k=m}^\\infty \\; (n_k(m)+1) \\; \\left (\\frac{a}{r} \\right)^{n_k(m)+2} \\cdot P_{n_k(m)}^m \\; (\\cos⁡ \\,\\theta_T) \\; \\{g_k^m  \\cos⁡\\,m\\varphi_T + h_k^m \\sin⁡\\,m\\varphi_T\\}\\end{align}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 128,
   "metadata": {},
   "outputs": [],
   "source": [
    "def Bz(n, m, P, phiT, thetaT):\n",
    "    Z_coeff = []\n",
    "    for r in range(len(m)):\n",
    "        if m[r]==0: # m=0, we only compute g\n",
    "            Z_coeff = Z_coeff + [(np.sin(thetaT)**m[r])*(-(n[r]+1)*P[r]*np.cos(m[r]*phiT))] # g\n",
    "        else: # Otherwise, we need to compute g and h\n",
    "            Z_coeff = Z_coeff + [(np.sin(thetaT)**m[r])*(-(n[r]+1)*P[r]*np.cos(m[r]*phiT))] # g\n",
    "            Z_coeff = Z_coeff + [(np.sin(thetaT)**m[r])*(-(n[r]+1)*P[r]*np.sin(m[r]*phiT))] # h     \n",
    "    return Z_coeff"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<p style='text-align: justify;'>The model parameter can be determined here</p>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 129,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Order of the model\n",
    "mmax = 7\n",
    "# Spherical cap parameter\n",
    "Latcentre = -3\n",
    "Longcentre = 122\n",
    "theta0 = 30\n",
    "year = 2015"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<p style='text-align: justify;'>The data used in this study are the X, Y, and Z components compiled from 68 geomagnetic repeat stations in Indonesia, covering the period 1985 - 2015 from BMKG (Badan Meteorologi Klimatologi dan Geofisika). However, two locations are not included in the calculation and will be used later to validate the model. These validation observatories are Bandar Lampung (105.175, -5.240) and Tual (132.735, -5.662). The repeat station data are combined with definitive data from five BMKG’s geomagnetic observatories from 1985 - 2019. To complete the dataset, the definitive data from a similar period from 13 INTERMAGNET observatories are used. In spherical cap harmonic modelling it is essential to have data around the edges of the chosen spherical cap. Because of the irregularities in the data distribution in the research area, the synthetic cartesian X, Y, and Z components at sea level at 17 fixed points are calculated from IGRF-13 at five-year intervals from 1985 – 2015. The IGRF-13 is also used to remove the main field before applying the techniques to calculate the regional field.</p>\n",
    "\n",
    "<img src=\"img/data_points.png\" width=\"400\"/>\n",
    "\n",
    "<p style='text-align: justify;'>In this research, the coordinate of the spherical cap pole is at 122$^{\\circ}$E and 3$^{\\circ}$S, and the spherical cap half-angle is 30$^{\\circ}$. This covers the Indonesia region, as shown in the following figure. According to Torta, Gaya-Piqué, & De Santis (2006), the use of a large spherical cap is necessary to obtain the main field and its secular variation; otherwise, the results are unrealistic. Research by Qamili et al. (2010) used a half-angle of 16$^{\\circ}$, although their data mostly concentrated in the central 4$^{\\circ}$ half-angle. Torta, García, Curto, & De Santis (1992) and Torta & De Santis (1996) explain that the RMS residual of the SCHA model decreases as the cap half-angle increase. However, Düzgit & Malin (2000) obtained different results and reported that the quality of the fit of the SCHA model is relatively insensitive to the cap radius. Nevertheless, the cap half-angle chosen in this research is the most optimal as it is cover the Indonesia region, and there are several geomagnetic observatories located near the boundary of the cap.</p>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 130,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load sample data\n",
    "with open(os.path.join(os.path.expanduser('~'), 'data/2015.pkl'), 'rb') as f:\n",
    "    df = pickle.load(f)\n",
    "with open(os.path.join(os.path.expanduser('~'), 'data/Grid_IGRF.pkl'), 'rb') as f:\n",
    "    data_grid_IGRF = pickle.load(f)\n",
    "# Load A and n from nk_solver.py\n",
    "with open(os.path.join(os.path.expanduser('~'), 'data/output_'+str(mmax)+'_'+str(theta0)+'.pkl'), 'rb') as f:\n",
    "    n,A,Kmn = pickle.load(f)\n",
    "\n",
    "# Gives names to different data columns (assuming order as above)\n",
    "Lat = df['Lat']\n",
    "Long = df['Lon']\n",
    "Xobs = df['X Obs']\n",
    "Yobs = df['Y Obs']\n",
    "Zobs = df['Z Obs']\n",
    "XIGRF = df['X IGRF']\n",
    "YIGRF = df['Y IGRF']\n",
    "ZIGRF = df['Z IGRF']\n",
    "\n",
    "# Remove main field\n",
    "Xres = Xobs-XIGRF\n",
    "Yres = Yobs-YIGRF\n",
    "Zres = Zobs-ZIGRF\n",
    "\n",
    "# Number of data\n",
    "num_data = len(Lat)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<p style='text-align: justify;'>Convert coordinate into radians</p>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 131,
   "metadata": {},
   "outputs": [],
   "source": [
    "m = []\n",
    "for a in range(mmax + 1):     # We need to add 1 because python list start at 0\n",
    "    b = list(range(a + 1))\n",
    "    m = m + b\n",
    "\n",
    "# Colatitude of the center of the cap\n",
    "colatradc = (90-Latcentre)*np.pi/180\n",
    "# Longitude of the center of the cap\n",
    "longradc = Longcentre*np.pi/180\n",
    "\n",
    "# spherical harmonic analysis works with colatitude, and in radians, so convert original data\n",
    "colat = 90-Lat\n",
    "# Colatitude of the data\n",
    "colatrad = colat*np.pi/180\n",
    "# Longitude of the data\n",
    "longrad = Long*np.pi/180\n",
    "# number of data sites\n",
    "num = len(colatrad)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<p style='text-align: justify;'>Firstly need to rotate given X, Y and Z values so that they are in the coordinate system of the cap</p>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 132,
   "metadata": {},
   "outputs": [],
   "source": [
    "colatcap = sph_cosine_side(colatradc, colatrad, longradc-longrad)\n",
    "longcap = np.zeros(num)\n",
    "gamma = np.zeros(num)\n",
    "\n",
    "# longitude in cap\n",
    "for i in range(num):\n",
    "    if Long[i]==Longcentre and Lat[i]>=Latcentre: # all points above centre on meridian\n",
    "        longcap[i] = 0\n",
    "    elif Long[i]==Longcentre and Lat[i]<Latcentre: # all points below centre on meridian\n",
    "        longcap[i] = np.pi\n",
    "    else:\n",
    "        longcap[i] = sph_cosine_angle(colatradc, colatcap[i], colatrad[i])\n",
    "    \n",
    "    # also need to rotate the X and Y values by angle gamma, such that:\n",
    "    if Long[i]==Longcentre and Lat[i]<=Latcentre:        # points below centre, on meridian\n",
    "        gamma[i] = 0.0\n",
    "    elif Long[i]==Longcentre and Lat[i]>Latcentre: # points above centre, on meridian\n",
    "        gamma[i] = np.pi\n",
    "    else:\n",
    "        gamma[i] = sph_cosine_angle(colatrad[i], colatcap[i], colatradc)\n",
    "        \n",
    "    # place the longitudes and gammas in right quadrant\n",
    "    if longrad[i]>longradc:\n",
    "        longcap[i] = 2*np.pi-longcap[i]\n",
    "        gamma[i] = 2*np.pi-gamma[i]\n",
    "        "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<p style='text-align: justify;'>Rotate the geomagnetic components. The transformed geomagnetic data for the X and Y components can be determined as<\\p>\n",
    "\n",
    "\\begin{align} B_X^\\prime=B_X \\cos\\,\\gamma + B_Y \\sin\\,\\gamma \\end{align}\n",
    "    \n",
    "\\begin{align} B_Y^\\prime=-B_X \\sin\\,\\gamma + B_Y \\cos\\,\\gamma \\end{align}\n",
    "    \n",
    "\\begin{align} B_Z^\\prime=B_Z \\end{align}\n",
    "    \n",
    "<p style='text-align: justify;'>The downward component ($B_Z$) is unaffected by the coordinate transformation</p>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 133,
   "metadata": {},
   "outputs": [],
   "source": [
    "Xcap = Xres*np.cos(gamma) + Yres*np.sin(gamma)\n",
    "Ycap = -Xres*np.sin(gamma) + Yres*np.cos(gamma)\n",
    "Zcap = Zres"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<p style='text-align: justify;'>Calculates the Gauss coefficients</p>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 134,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Each coefficient consists of two parts, each with a summation over the A(m,n) values as given in the matrix above.\n",
    "field = np.zeros(3*num)\n",
    "coeff = list(np.zeros((3*num)))\n",
    "\n",
    "for i in range(num):\n",
    "    # Associate Legendre function\n",
    "    P = assoc_Leg_func(n, A, colatcap[i])\n",
    "    # Derived associate Legendre function\n",
    "    dP = deriv_assoc_Leg_func(n, m, A, colatcap[i], P)\n",
    "    \n",
    "    # X coefficients are 1/r dU/d(colatcap)\n",
    "    coeff[i] = Bx(m, dP, longcap[i], colatcap[i])\n",
    "    field[i] = Xcap[i]\n",
    "    \n",
    "    # Y coefficients are -1/rsin(colatcap) dU/d(longcap)\n",
    "    coeff[num+i] = By(n, m, P, longcap[i], colatcap[i])\n",
    "    field[num+i] = Ycap[i]\n",
    "    \n",
    "    # Z coefficients are du/dr\n",
    "    coeff[2*num+i] = Bz(n, m, P, longcap[i], colatcap[i])\n",
    "    field[2*num+i] = Zcap[i]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<p style='text-align: justify;'>The gauss coefficients are found by a least squares method for the over determined set of equations.</p>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 135,
   "metadata": {},
   "outputs": [],
   "source": [
    "gauss = np.linalg.lstsq(coeff,field,rcond=None)[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<p style='text-align: justify;'>Calculates the statistics of the model (Mean Absolute Error (MAE) and Root Mean Square Deviation (RMSD))</p>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 136,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate the Mean Absolute Error (MAE) or mean difference\n",
    "# Need to calculate the modeled field first (field=model*gauss)\n",
    "Xmodcap = np.dot(coeff[0:num_data],gauss)              # X model field in cap coordinate\n",
    "Ymodcap = np.dot(coeff[num:num+num_data],gauss)      # Y model field in cap coordinate\n",
    "Zmodcap = np.dot(coeff[2*num:2*num+num_data],gauss)  # Z model field in cap coordinate\n",
    "\n",
    "Xmod = Xmodcap*np.cos(gamma[0:num_data])-Ymodcap*np.sin(gamma[0:num_data])   # X model field in geographic coordinate\n",
    "Ymod = Xmodcap*np.sin(gamma[0:num_data])+Ymodcap*np.cos(gamma[0:num_data])  # Y model field in geographic coordinate\n",
    "Zmod = Zmodcap                           # Z model field in geographic coordinate\n",
    "\n",
    "XMAE = np.mean(abs(Xmod-Xres[0:num_data]))\n",
    "YMAE = np.mean(abs(Ymod-Yres[0:num_data]))\n",
    "ZMAE = np.mean(abs(Zmod-Zres[0:num_data]))\n",
    "\n",
    "# Calculate the Root Mean Square Deviation (RMSD) or RMS difference\n",
    "XRMSD = np.sqrt(np.mean((Xmod-Xres[0:num_data])**2))\n",
    "YRMSD = np.sqrt(np.mean((Ymod-Yres[0:num_data])**2))\n",
    "ZRMSD = np.sqrt(np.mean((Zmod-Zres[0:num_data])**2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<p style='text-align: justify;'>Calculates the SCHA model over the whole cap on a grid</p>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 137,
   "metadata": {},
   "outputs": [],
   "source": [
    "# As with the spherical harmonic code, want to calculate the three magnetic field components across a grid of the Indonesia region, but this grid needs to be within the cap.\n",
    "Xcoeff = []\n",
    "Ycoeff = []\n",
    "Zcoeff = []\n",
    "count = -1\n",
    "gridcolat = []\n",
    "gridlong = []\n",
    "gm = []\n",
    "\n",
    "for i in range(Latcentre-theta0,Latcentre+theta0+1):          # runs through latitudes\n",
    "    for j in range(Longcentre-theta0,Longcentre+theta0+1):    # runs through longitudes\n",
    "        dist = distance_on_unit_sphere(i,j,Latcentre,Longcentre)*180/math.pi\n",
    "        # calculates the distance in deg between the grid point and centre of cap\n",
    "        if dist<=theta0:                                    # eliminates any grid points that are outside the cap and eliminates cap longitude 90 and 270 deg points as these cause singularities\n",
    "            count+=1\n",
    "            cl = (90-i)*math.pi/180  # puts latitude into colat in radians\n",
    "            l = j*math.pi/180        # puts longitude into radians\n",
    "            gridcolat = gridcolat + [cl]    # makes vectors of the grid colats and longs (original) for later use\n",
    "            gridlong = gridlong + [l]\n",
    "            # colat in cap:\n",
    "            clc = sph_cosine_side(colatradc, cl, longradc-l)\n",
    "            # longitude in cap\n",
    "            if j==Longcentre and i>=Latcentre:\n",
    "                lc = 0\n",
    "            elif j==Longcentre and i<Latcentre:\n",
    "                lc = math.pi\n",
    "            else:\n",
    "                lc = sph_cosine_angle(colatradc, clc, cl)\n",
    "            # also need to rotate the X and Y values by angle gamma, such that:\n",
    "            if j==Longcentre and cl>=colatradc:\n",
    "                gm = gm + [0]\n",
    "            elif j==Longcentre and cl<colatradc:\n",
    "                gm = gm + [math.pi]\n",
    "            else:\n",
    "                gm = gm + [sph_cosine_angle(cl, clc, colatradc)]\n",
    "                \n",
    "            # need to determine which quadrant the longitude and rotation is in\n",
    "            if l>longradc:\n",
    "                lc = 2*math.pi-lc\n",
    "                gm[count] = 2*math.pi-gm[count]\n",
    "                \n",
    "            # now calculate sums for coefficients\n",
    "            P = assoc_Leg_func(n, A, clc)\n",
    "            dP = deriv_assoc_Leg_func(n, m, A, clc, P)\n",
    "            \n",
    "            # X coefficients are 1/r dU/d(colatcap)\n",
    "            Xcoeff = Xcoeff + [Bx(m, dP, lc, clc)]\n",
    "          \n",
    "            # Y coefficients are -1/rsin(colatcap) dU/d(longcap)\n",
    "            Ycoeff = Ycoeff + [By(n, m, P, lc, clc)]            \n",
    "            \n",
    "            # Z coefficients are du/dr\n",
    "            Zcoeff = Zcoeff + [Bz(n, m, P, lc, clc)]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<p style='text-align: justify;'>Now to calculate the model field values at each point is a matrix problem. (coeffs) x (gauss) = (field)</p>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 138,
   "metadata": {},
   "outputs": [],
   "source": [
    "Xcalccap = np.dot(np.array(Xcoeff),gauss)\n",
    "Ycalccap = np.dot(np.array(Ycoeff),gauss)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<p style='text-align: justify;'>Before plotting, need to revert the X and Y values back into geographic coordinates, not cap coords. Calculates the other components as well.</p>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 139,
   "metadata": {},
   "outputs": [],
   "source": [
    "Xcalc = Xcalccap*np.cos(gm)-Ycalccap*np.sin(gm)+data_grid_IGRF.iloc[0:len(Xcoeff), 3]\n",
    "Ycalc = Xcalccap*np.sin(gm)+Ycalccap*np.cos(gm)+data_grid_IGRF.iloc[0:len(Xcoeff), 4]\n",
    "Zcalc = np.dot(np.array(Zcoeff),gauss)+data_grid_IGRF.iloc[0:len(Xcoeff), 5]\n",
    "Dcalc = np.arctan(Ycalc/Xcalc)*180/np.pi\n",
    "Icalc = np.arctan(Zcalc/np.sqrt(Xcalc**2+Ycalc**2))*180/np.pi\n",
    "Fcalc = np.sqrt(Xcalc**2+Ycalc**2+Zcalc**2)\n",
    "Hcalc = np.sqrt(Xcalc**2+Ycalc**2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<p style='text-align: justify;'>Puts colatitude and longitude back to latitude and longitude in degrees</p>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 140,
   "metadata": {},
   "outputs": [],
   "source": [
    "gridlatd = 90-np.dot(gridcolat,180/np.pi)\n",
    "gridlongd = np.dot(gridlong,180/np.pi)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<p style='text-align: justify;'>Show G and H gauss coefficients</p>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 141,
   "metadata": {},
   "outputs": [],
   "source": [
    "kk = []\n",
    "for a in range(mmax + 1):\n",
    "    b = np.tile(a,(1,a+1))[0].tolist()\n",
    "    kk = kk + b\n",
    "\n",
    "gauss_g = np.zeros(len(kk))\n",
    "gauss_h = np.zeros(len(kk))\n",
    "q = np.zeros(len(kk)-mmax-1)\n",
    "g = 0\n",
    "h = 0\n",
    "p = -1\n",
    "\n",
    "for c in range(len(kk)):\n",
    "    gauss_g[c] = gauss[g]\n",
    "    if m[c]==0:\n",
    "        g = g+1\n",
    "    else:\n",
    "        g = g+2\n",
    "\n",
    "for c in range(len(kk)-mmax-1):\n",
    "    if m[c]==0:\n",
    "        h = h+3\n",
    "    else:\n",
    "        h = h+2\n",
    "    q[c] = h\n",
    "\n",
    "for c in range(len(kk)):\n",
    "    if m[c]==0:\n",
    "        gauss_h[c] = np.nan\n",
    "    else:\n",
    "        p = p+1\n",
    "        gauss_h[c] = gauss[np.int(q[p])]\n",
    "        \n",
    "df_coeff = pd.DataFrame(index=range(len(kk)), columns=['k','m','nkm','kmn','gmk','hmk'])\n",
    "df_coeff['k'] = kk\n",
    "df_coeff['m'] = m\n",
    "df_coeff['nkm'] = ['%.4f' % elem for elem in n]\n",
    "df_coeff['kmn'] = ['%.4f' % elem for elem in Kmn]\n",
    "df_coeff['gmk'] = ['%.4f' % elem for elem in gauss_g]\n",
    "df_coeff['hmk'] = ['%.4f' % elem for elem in gauss_h]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<p style='text-align: justify;'>Saving output file and print the statistics of the model.</p>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 142,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Statistical results of the SCHA model with kmax = 7\n",
      "\n",
      "\t\t\t\t\tX(nT)\tY(nT)\tZ(nT)\n",
      "Mean Absolute Error (MAE)\t\t  99.87\t  81.01\t 110.25\n",
      "Root Mean Square Difference (RMSD)\t 130.61\t 107.61\t 156.67\n",
      "\n",
      "\n",
      "\n",
      "Coefficients of the Spherical Cap Harmonic Analysis of the Indonesian Geomagnetic Data at the Epoch 2015 for spherical cap = 30°\n",
      "\n",
      "k = cap harmonic degree\n",
      "m = cap harmonic degree\n",
      "nkm = non-integer degree of the Legendre function\n",
      "kmn = Schmidt normalisation constant\n",
      "gmk, hmk = cap harmonic coefficients\n",
      "\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>k</th>\n",
       "      <th>m</th>\n",
       "      <th>nkm</th>\n",
       "      <th>kmn</th>\n",
       "      <th>gmk</th>\n",
       "      <th>hmk</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0.0000</td>\n",
       "      <td>1.0000</td>\n",
       "      <td>-484.4968</td>\n",
       "      <td></td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>4.0837</td>\n",
       "      <td>1.0000</td>\n",
       "      <td>-318.1156</td>\n",
       "      <td></td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>3.1196</td>\n",
       "      <td>2.5362</td>\n",
       "      <td>-863.6470</td>\n",
       "      <td>436.7040</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>2</td>\n",
       "      <td>0</td>\n",
       "      <td>6.8354</td>\n",
       "      <td>1.0000</td>\n",
       "      <td>292.3059</td>\n",
       "      <td></td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>2</td>\n",
       "      <td>1</td>\n",
       "      <td>6.8354</td>\n",
       "      <td>5.1775</td>\n",
       "      <td>1049.8816</td>\n",
       "      <td>-1087.9306</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>2</td>\n",
       "      <td>2</td>\n",
       "      <td>5.4928</td>\n",
       "      <td>6.1253</td>\n",
       "      <td>312.8332</td>\n",
       "      <td>-775.9311</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>3</td>\n",
       "      <td>0</td>\n",
       "      <td>10.0386</td>\n",
       "      <td>1.0000</td>\n",
       "      <td>-257.2242</td>\n",
       "      <td></td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>3</td>\n",
       "      <td>1</td>\n",
       "      <td>9.7121</td>\n",
       "      <td>7.2160</td>\n",
       "      <td>-1123.3533</td>\n",
       "      <td>1229.1932</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>3</td>\n",
       "      <td>2</td>\n",
       "      <td>9.3733</td>\n",
       "      <td>17.0110</td>\n",
       "      <td>-692.9836</td>\n",
       "      <td>616.9469</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>3</td>\n",
       "      <td>3</td>\n",
       "      <td>7.7524</td>\n",
       "      <td>15.4888</td>\n",
       "      <td>634.6269</td>\n",
       "      <td>-92.5589</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10</th>\n",
       "      <td>4</td>\n",
       "      <td>0</td>\n",
       "      <td>12.9083</td>\n",
       "      <td>1.0000</td>\n",
       "      <td>194.0154</td>\n",
       "      <td></td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>11</th>\n",
       "      <td>4</td>\n",
       "      <td>1</td>\n",
       "      <td>12.9083</td>\n",
       "      <td>9.4793</td>\n",
       "      <td>844.4350</td>\n",
       "      <td>-995.5895</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>12</th>\n",
       "      <td>4</td>\n",
       "      <td>2</td>\n",
       "      <td>12.3720</td>\n",
       "      <td>29.0690</td>\n",
       "      <td>860.3532</td>\n",
       "      <td>-523.9351</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>13</th>\n",
       "      <td>4</td>\n",
       "      <td>3</td>\n",
       "      <td>11.8074</td>\n",
       "      <td>53.3354</td>\n",
       "      <td>-608.5895</td>\n",
       "      <td>257.2099</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>14</th>\n",
       "      <td>4</td>\n",
       "      <td>4</td>\n",
       "      <td>9.9589</td>\n",
       "      <td>39.8614</td>\n",
       "      <td>312.4701</td>\n",
       "      <td>527.4549</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>15</th>\n",
       "      <td>5</td>\n",
       "      <td>0</td>\n",
       "      <td>16.0248</td>\n",
       "      <td>1.0000</td>\n",
       "      <td>-102.9499</td>\n",
       "      <td></td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>16</th>\n",
       "      <td>5</td>\n",
       "      <td>1</td>\n",
       "      <td>15.8215</td>\n",
       "      <td>11.5415</td>\n",
       "      <td>-493.3762</td>\n",
       "      <td>594.6422</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>17</th>\n",
       "      <td>5</td>\n",
       "      <td>2</td>\n",
       "      <td>15.6154</td>\n",
       "      <td>45.6896</td>\n",
       "      <td>-695.4600</td>\n",
       "      <td>293.1358</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>18</th>\n",
       "      <td>5</td>\n",
       "      <td>3</td>\n",
       "      <td>14.9180</td>\n",
       "      <td>105.9937</td>\n",
       "      <td>555.6899</td>\n",
       "      <td>-417.0814</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>19</th>\n",
       "      <td>5</td>\n",
       "      <td>4</td>\n",
       "      <td>14.1778</td>\n",
       "      <td>162.6264</td>\n",
       "      <td>-203.0030</td>\n",
       "      <td>-338.3575</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>20</th>\n",
       "      <td>5</td>\n",
       "      <td>5</td>\n",
       "      <td>12.1334</td>\n",
       "      <td>103.4865</td>\n",
       "      <td>-70.2786</td>\n",
       "      <td>-872.9003</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>21</th>\n",
       "      <td>6</td>\n",
       "      <td>0</td>\n",
       "      <td>18.9364</td>\n",
       "      <td>1.0000</td>\n",
       "      <td>38.5086</td>\n",
       "      <td></td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>22</th>\n",
       "      <td>6</td>\n",
       "      <td>1</td>\n",
       "      <td>18.9364</td>\n",
       "      <td>13.7460</td>\n",
       "      <td>172.7607</td>\n",
       "      <td>-202.5747</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>23</th>\n",
       "      <td>6</td>\n",
       "      <td>2</td>\n",
       "      <td>18.5830</td>\n",
       "      <td>64.1553</td>\n",
       "      <td>394.7865</td>\n",
       "      <td>-100.3626</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>24</th>\n",
       "      <td>6</td>\n",
       "      <td>3</td>\n",
       "      <td>18.2219</td>\n",
       "      <td>190.9258</td>\n",
       "      <td>-315.6057</td>\n",
       "      <td>358.1270</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>25</th>\n",
       "      <td>6</td>\n",
       "      <td>4</td>\n",
       "      <td>17.3910</td>\n",
       "      <td>364.9771</td>\n",
       "      <td>64.3862</td>\n",
       "      <td>372.7627</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>26</th>\n",
       "      <td>6</td>\n",
       "      <td>5</td>\n",
       "      <td>16.5041</td>\n",
       "      <td>486.5812</td>\n",
       "      <td>64.1954</td>\n",
       "      <td>748.2361</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>27</th>\n",
       "      <td>6</td>\n",
       "      <td>6</td>\n",
       "      <td>14.2861</td>\n",
       "      <td>270.0018</td>\n",
       "      <td>-361.1992</td>\n",
       "      <td>595.8797</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>28</th>\n",
       "      <td>7</td>\n",
       "      <td>0</td>\n",
       "      <td>22.0183</td>\n",
       "      <td>1.0000</td>\n",
       "      <td>-4.2517</td>\n",
       "      <td></td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>29</th>\n",
       "      <td>7</td>\n",
       "      <td>1</td>\n",
       "      <td>21.8702</td>\n",
       "      <td>15.8222</td>\n",
       "      <td>-34.6558</td>\n",
       "      <td>28.8938</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>30</th>\n",
       "      <td>7</td>\n",
       "      <td>2</td>\n",
       "      <td>21.7210</td>\n",
       "      <td>87.0682</td>\n",
       "      <td>-109.1371</td>\n",
       "      <td>18.2356</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>31</th>\n",
       "      <td>7</td>\n",
       "      <td>3</td>\n",
       "      <td>21.2461</td>\n",
       "      <td>300.1786</td>\n",
       "      <td>121.4010</td>\n",
       "      <td>-133.0427</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>32</th>\n",
       "      <td>7</td>\n",
       "      <td>4</td>\n",
       "      <td>20.7588</td>\n",
       "      <td>734.7575</td>\n",
       "      <td>13.7400</td>\n",
       "      <td>-166.3312</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>33</th>\n",
       "      <td>7</td>\n",
       "      <td>5</td>\n",
       "      <td>19.8121</td>\n",
       "      <td>1210.1830</td>\n",
       "      <td>156.0653</td>\n",
       "      <td>-397.3139</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>34</th>\n",
       "      <td>7</td>\n",
       "      <td>6</td>\n",
       "      <td>18.7979</td>\n",
       "      <td>1435.9065</td>\n",
       "      <td>454.2694</td>\n",
       "      <td>-105.5000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>35</th>\n",
       "      <td>7</td>\n",
       "      <td>7</td>\n",
       "      <td>16.4228</td>\n",
       "      <td>706.5894</td>\n",
       "      <td>52.9736</td>\n",
       "      <td>182.2293</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "    k  m      nkm        kmn         gmk         hmk\n",
       "0   0  0   0.0000     1.0000   -484.4968            \n",
       "1   1  0   4.0837     1.0000   -318.1156            \n",
       "2   1  1   3.1196     2.5362   -863.6470    436.7040\n",
       "3   2  0   6.8354     1.0000    292.3059            \n",
       "4   2  1   6.8354     5.1775   1049.8816  -1087.9306\n",
       "5   2  2   5.4928     6.1253    312.8332   -775.9311\n",
       "6   3  0  10.0386     1.0000   -257.2242            \n",
       "7   3  1   9.7121     7.2160  -1123.3533   1229.1932\n",
       "8   3  2   9.3733    17.0110   -692.9836    616.9469\n",
       "9   3  3   7.7524    15.4888    634.6269    -92.5589\n",
       "10  4  0  12.9083     1.0000    194.0154            \n",
       "11  4  1  12.9083     9.4793    844.4350   -995.5895\n",
       "12  4  2  12.3720    29.0690    860.3532   -523.9351\n",
       "13  4  3  11.8074    53.3354   -608.5895    257.2099\n",
       "14  4  4   9.9589    39.8614    312.4701    527.4549\n",
       "15  5  0  16.0248     1.0000   -102.9499            \n",
       "16  5  1  15.8215    11.5415   -493.3762    594.6422\n",
       "17  5  2  15.6154    45.6896   -695.4600    293.1358\n",
       "18  5  3  14.9180   105.9937    555.6899   -417.0814\n",
       "19  5  4  14.1778   162.6264   -203.0030   -338.3575\n",
       "20  5  5  12.1334   103.4865    -70.2786   -872.9003\n",
       "21  6  0  18.9364     1.0000     38.5086            \n",
       "22  6  1  18.9364    13.7460    172.7607   -202.5747\n",
       "23  6  2  18.5830    64.1553    394.7865   -100.3626\n",
       "24  6  3  18.2219   190.9258   -315.6057    358.1270\n",
       "25  6  4  17.3910   364.9771     64.3862    372.7627\n",
       "26  6  5  16.5041   486.5812     64.1954    748.2361\n",
       "27  6  6  14.2861   270.0018   -361.1992    595.8797\n",
       "28  7  0  22.0183     1.0000     -4.2517            \n",
       "29  7  1  21.8702    15.8222    -34.6558     28.8938\n",
       "30  7  2  21.7210    87.0682   -109.1371     18.2356\n",
       "31  7  3  21.2461   300.1786    121.4010   -133.0427\n",
       "32  7  4  20.7588   734.7575     13.7400   -166.3312\n",
       "33  7  5  19.8121  1210.1830    156.0653   -397.3139\n",
       "34  7  6  18.7979  1435.9065    454.2694   -105.5000\n",
       "35  7  7  16.4228   706.5894     52.9736    182.2293"
      ]
     },
     "execution_count": 142,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Save the file\n",
    "fileout = 'tmp\\\\'+str(year)+'_'+str(mmax)+'_'+str(theta0)+'.txt'\n",
    "f_out = open(fileout, 'w')\n",
    "\n",
    "f_out.write('Year    Long      Lat     X         Y        Z          H            D        I       F\\n')\n",
    "for s in range(len(gridlatd)):\n",
    "    out = '%4.0f %9.2f %9.2f %9.2f %9.2f %9.2f %9.2f %9.2f %9.2f %9.2f\\n' %(year,gridlongd[s],gridlatd[s],Xcalc[s],Ycalc[s],Zcalc[s],Hcalc[s],Dcalc[s],Icalc[s],Fcalc[s])\n",
    "    f_out.write(out)\n",
    "f_out.close()\n",
    "\n",
    "print('Statistical results of the SCHA model with kmax = '+str(mmax)+'\\n')\n",
    "print('\\t\\t\\t\\t\\tX(nT)\\tY(nT)\\tZ(nT)\\nMean Absolute Error (MAE)\\t\\t%7.2f\\t%7.2f\\t%7.2f\\nRoot Mean Square Difference (RMSD)\\t%7.2f\\t%7.2f\\t%7.2f\\n\\n\\n'%(XMAE,YMAE,ZMAE,XRMSD,YRMSD,ZRMSD))\n",
    "print('Coefficients of the Spherical Cap Harmonic Analysis of the Indonesian Geomagnetic Data at the Epoch '+str(year)+' for spherical cap = '+str(theta0)+'\\u00b0\\n')\n",
    "print('k = cap harmonic degree')\n",
    "print('m = cap harmonic degree')\n",
    "print('nkm = non-integer degree of the Legendre function')\n",
    "print('kmn = Schmidt normalisation constant')\n",
    "print('gmk, hmk = cap harmonic coefficients\\n')\n",
    "\n",
    "df_coeff.replace('nan', '', regex=True)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.7.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}