[FEAT]: Support ghost cells in dipole grid

- I can't test in 1D since things break in the geographic grid. If the dipole grid is broken in 1D, let me know!

The altitude (along-field-line) ghost cells are treated as standard grid points.
- At low altitudes, points are close enough this is OK.
- At high altitudes, it's going to be hard to get ghost cells.
  - For example, if they extend along the field line, there will be points at +25 Re (L-shell of max_lat)
  - If they are scaled to lower radius, or just not as far along the field line, they will not be perpendicular to B anymore!

> Maybe specify where field lines should be separated btwn open/closed? Then closed field lines will overlap across equator, and open field lines will not be perpendicular to B in the ghost cells? IDK yet!

----

Ghost cells in latitude currently go over the poles. Discussion is happening on whether this is a good idea or not.

---

Author: abukowski21

src/init_mag_grid.cpp

@@ -59,6 +59,9 @@ void Grid::init_dipole_grid(Quadtree quadtree, Planets planet)
 
   Inputs::grid_input_struct grid_input = input.get_grid_inputs("ionGrid");
 
+  // Number of ghost cells:
+  int64_t nGCs = get_nGCs();
+
   // Get inputs
   report.print(3, "Setting inputs in dipole grid");
   precision_t min_apex = grid_input.min_apex * cKMtoM;
@@ -85,6 +88,11 @@ void Grid::init_dipole_grid(Quadtree quadtree, Planets planet)
   precision_t dlon = size_right_norm(0) * cPI / (nLons - 2 * nGCs);
   precision_t lon0 = lower_left_norm(0) * cPI;
   arma_vec lon1d(nLons);
+
+  // Quick 1D check for longitudes:
+  // Same as geo_grid here...
+  if (!HasXdim) dlon = 1.0 * cDtoR;
+
   for (iLon = 0; iLon < nLons; iLon++)
     lon1d(iLon) = lon0 + (iLon - nGCs + 0.5) * dlon;
 
@@ -96,7 +104,7 @@ void Grid::init_dipole_grid(Quadtree quadtree, Planets planet)
 
   geoLon_scgc = magLon_scgc;
 
-  // LATITUDES
+  // LATITUDES:
 
   // min_latitude calculated from min_lshell (min_apex input)
   precision_t min_lat = get_lat_from_r_and_lshell(1.0, min_lshell);
@@ -105,38 +113,44 @@ void Grid::init_dipole_grid(Quadtree quadtree, Planets planet)
   // some intermediates:
   precision_t del_lat, blat_min_, blat_max_, tmp_lat;
   // Integers for field-line loops:
-  // - nF=nLats/2 (so nLarts MUST be even)
-  // - nZ=nAlts
-  // - nZby2 = nZ/2
-  int64_t nF = nLats / 2, nZ = nAlts*2, nZby2 = nAlts;
+  // - nF=nLats/2 (so nLats MUST be even)
+  // - nZ=nAlts*2 - (we fill up HALF field lines)
+  // - nZby2 = nAlts
+  // -> Plus, experinemtal support for altitude ghost cells...
+  int64_t nF = (nLats) / 2 , nZ = (nAlts) * 2, nZby2 = (nAlts);
   // lShells and baseLats are currently set for southern hemisphere then mirrored
   arma_vec Lshells(nF), baseLats(nF);
 
+  blat_min_ = cos(pow(min_lat, 1.0 / LatStretch));
+  blat_max_ = cos(pow(max_lat, 1.0 / LatStretch));
+  del_lat = (blat_max_ - blat_min_) / (nF - nGCs * 2.0);
+
   // now make sure the user used 1 or an even number for nLats
   if (nLats % 2 != 0)
   {
-    if (nLats == 1)
+    if (!HasYdim) 
     {
-        report.print(0, ">> Running in 1D. Experinental!!");
+      del_lat = 1.0 * cDtoR;
+      report.print(0, "Running in single latitude dimension. Experinental!!");
       nF = 1;
     }
     else
       report.error("Cannot use odd nLats with dipole grid!");
   }
 
-  blat_min_ = cos(pow(min_lat, 1.0 / LatStretch));
-  blat_max_ = cos(pow(max_lat, 1.0 / LatStretch));
-  del_lat = (blat_max_ - blat_min_) / (nF - nGCs * 2.0);
-
+  // loop over all cells - everything including the ghost cells
+  // -> This means some points will go over the pole (baseLat > +- 90 degrees)
+  //    They're taken care of in the conversion to geographic coordinates.
   for (int i = 0; i < nF; i++)
   {
     // first put down "linear" spacing
     tmp_lat = blat_min_ + del_lat * (i - nGCs + 0.5);
-    // then scale it according to the exponent & acos it
+    // then scale it according to the exponent & acos
     tmp_lat = pow(acos(tmp_lat), LatStretch);
     // place values in array backwards, South pole -> equator.
     baseLats(nF - i - 1) = -tmp_lat;
   }
+  report.print(3, "Done setting base latitudes for dipole grid.");
 
   // Find L-Shell for each baseLat
   // using L=R/sin2(theta), where theta is from north pole
@@ -176,8 +190,7 @@ void Grid::init_dipole_grid(Quadtree quadtree, Planets planet)
     // inlo loop thru southern hemisphere, mirror in  north.
     for (int i_nZ = 0; i_nZ < nAlts; i_nZ++)
     {
-      // Should be using lshell_to_qn_qs, but it wasn't working right.
-      // This does the same, but won't work for offset dipoles.
+      // This won't work for offset dipoles.
       // Doesn't have any lat/lon dependence.
       delqp = (q_N - q_S) / (nZ + 1);
       qp0 = q_S + i_nZ * (delqp);
@@ -194,6 +207,7 @@ void Grid::init_dipole_grid(Quadtree quadtree, Planets planet)
       bLats(i_nF, i_nZ) = r_theta.second;
     }
   }
+  report.print(3, "Done generating q-spacing for dipole grid.");
 
   arma_vec rNorm1d(nAlts), lat1dAlong(nAlts);
   arma_cube r3d(nLons, nLats, nAlts);
@@ -218,11 +232,13 @@ void Grid::init_dipole_grid(Quadtree quadtree, Planets planet)
       magLat_scgc.tube(iLon, nLats - iLat - 1) = -lat1dAlong;
     }
   }
+  report.print(3, "Done generating symmetric latitude & altitude spacing in dipole.");
 
   geoLat_scgc = magLat_scgc;
   magAlt_scgc = r3d - planetRadius;
   geoAlt_scgc = magAlt_scgc;
 
+  report.print(4, "Beginning coordinate transformations of the dipole grid.");
   // Calculate the radius, of planet
   fill_grid_radius(planet);
 
@@ -256,6 +272,8 @@ void Grid::init_dipole_grid(Quadtree quadtree, Planets planet)
       gravity_vcgc[1] % gravity_vcgc[1] +
       gravity_vcgc[2] % gravity_vcgc[2]);
 
+  report.print(4, "Done gravity calculations for the dipole grid.");
+
   std::vector<arma_cube> llr, xyz, xyzRot1, xyzRot2;
   llr.push_back(magLon_scgc);
   llr.push_back(magLat_scgc);
@@ -275,11 +293,14 @@ void Grid::init_dipole_grid(Quadtree quadtree, Planets planet)
   geoLon_scgc = llr[0];
   geoLat_scgc = llr[1];
   geoAlt_scgc = llr[2] - planetRadius;
+  report.print(4, "Done dipole -> geographic transformations for the dipole grid.");
 
   calc_alt_grid_spacing();
+  report.print(4, "Done altitude spacing for the dipole grid.");
 
   // Calculate magnetic field and magnetic coordinates:
   fill_grid_bfield(planet);
+  report.print(4, "Done filling dipole grid with b-field!");
 
   report.exit(function);
   return;