FEAT: Better handling of alt gc's. And add mat of first ghost cell.

Author: abukowski21

include/grid.h

@@ -102,10 +102,14 @@ class Grid
   arma_cube magQ_Corner;
   arma_cube magInvLat_Corner;
 
-  // Dipole grid has cells below the surface of earth. These variables hold masks
-  // to either access those cells or ignore them (with .elem()).
+  // Masks to either access the non-physical (ghost) cells, or ignore them - use with
+  // .elem()). Together they *should* hold the indices of all cells.
   arma::uvec isTooLowCell, isPhysicalCell;
-  arma_cube UseThisCell; // (bool values whether altitude is valid)
+  // (bool values whether altitude is valid)
+  arma_cube UseThisCell;
+  // Matrices whose elements denote the altitude index of the interiormost ghost cell
+  // in the k-up and k-down direction (altitude for geo grids, q for dipole).
+  arma_mat first_lower_gc, first_upper_gc;
 
   // These are the locations of the magnetic poles:
   //  ll -> lat, lon, radius independent

src/grid.cpp

@@ -12,6 +12,7 @@
 Grid::Grid(std::string gridtype) {
 
   // At this point, we only need 2 ghostcells.  Hardcode this:
+  // This is also (kinda?) set in sizes.h for the geo & mag grid independently
   nGCs = 2;
 
   Inputs::grid_input_struct grid_input = input.get_grid_inputs(gridtype);
@@ -272,6 +273,8 @@ Grid::Grid(std::string gridtype) {
 
   UseThisCell.set_size(nX, nY, nZ);
   UseThisCell.fill(true);
+  first_lower_gc.set_size(nX, nY);
+  first_upper_gc.set_size(nX, nY);
 
   cent_acc_vcgc = make_cube_vector(nLons, nLats, nAlts, 3);
 

src/init_geo_grid.cpp

@@ -151,8 +151,20 @@ void Grid::create_altitudes(Planets planet) {
   }
 
   // All cells on the geographic grid *should* be ok
-  isTooLowCell = find(geoAlt_scgc < 0.0);
-  isPhysicalCell = find(geoAlt_scgc > 0.0);
+  isTooLowCell = find(geoAlt_scgc < grid_input.alt_min * cKMtoM);
+  isPhysicalCell = find(geoAlt_scgc >= grid_input.alt_min * cKMtoM);
+  // get the ghost cell indices on each lat/lon point. 
+  // may be redundant can fill lower with nGCs-1, but this is here for now
+  arma::uvec theGCs;
+  for (iLon=0; iLon<nLons; iLon++){
+    for (iLat = 0; iLat<nLats; iLat++){
+      // find *last* cell below alt_min
+      theGCs = find(geoAlt_scgc.tube(iLon, iLat) < grid_input.alt_min * cKMtoM);
+      // Get the last element if the col-vec
+      first_lower_gc(iLon, iLat) = theGCs(theGCs.n_elem - 1);
+    }
+  }
+  first_upper_gc.fill(nAlts - nGCs * 2 - 1);
 
   report.exit(function);
   return;

src/init_mag_grid.cpp

@@ -506,10 +506,21 @@ bool Grid::init_dipole_grid(Quadtree quadtree_ion, Planets planet) {
   // Generate mask for physicsl cells //
   //////////////////////////////////////
 
-  isTooLowCell = find(geoAlt_scgc <= 0.0);
-  isPhysicalCell = find(geoAlt_scgc > 0.0);
+  isTooLowCell = find(geoAlt_scgc < grid_input.alt_min * cKMtoM);
+  isPhysicalCell = find(geoAlt_scgc >= grid_input.alt_min * cKMtoM);
   UseThisCell.elem(isTooLowCell).fill(false);
 
+  arma::uvec theGCs;
+  for (iLon=0; iLon<nLons; iLon++){
+    for (iLat = 0; iLat<nLats; iLat++){
+      // find *last* cell below alt_min
+      theGCs = find(geoAlt_scgc.tube(iLon, iLat) < grid_input.alt_min * cKMtoM);
+      // Get the last element if the col-vec
+      first_lower_gc(iLon, iLat) = theGCs(theGCs.n_elem - 1);
+    }
+  }
+  first_upper_gc.fill(nAlts - nGCs * 2 - 1);
+
   report.print(4, "Done altitude spacing for the dipole grid.");
 
   // Calculate magnetic field and magnetic coordinates: