sheets.h

#ifndef COSMO_SHEET_H
#define COSMO_SHEET_H

#include "../../utils/RK4Register.h"
#include "../../utils/Array.h"
#include "../../utils/Timer.h"
#include "../../cosmo_types.h"
#include "../bssn/bssn.h"
#include "../../utils/TriCubicInterpolator.h"
#include "../Lambda/lambda.h"
#include <cmath>
#include <algorithm>
#include <iostream>
#include <string>
#include "zlib.h"
#include <limits>
#include <sstream>
#include <iomanip>
#include <fstream>
#include <vector>

namespace cosmo
{


class Sheet
{
public:
  // Simulation information
  idx_t nx, ny, nz;
  idx_t ns1, ns2, ns3; 
  real_t lx, ly, lz; 
  real_t dx, dy, dz; 
  
  register_t Dx, Dy, Dz; 
  register_t vx, vy, vz; 
  
  arr_t tmp; 

  bool follow_null_geodesics;
  real_t rescale_sheet;
  real_t ray_bundle_epsilon, det_g_obs;

  idx_t step;

  // internal types
  enum carrierCountScheme { per_dx = 0, per_ds = 1};
  carrierCountScheme carrier_count_scheme;

  enum depositScheme { CIC = 0, PCS = 1, CINT = 2 };
  depositScheme deposit;

  idx_t carriers_per_dx,
        carriers_per_dy,
        carriers_per_dz;

  Sheet();
  ~Sheet();

  void setDt(real_t dt);

  real_t _S1IDXtoX0(idx_t s1) { return (real_t)s1*lx/ns1; }
  real_t _S2IDXtoY0(idx_t s2) { return (real_t)s2*ly/ns2; }
  real_t _S3IDXtoZ0(idx_t s3) { return (real_t)s3*lz/ns3; }
  real_t _S1IDXtoX0(real_t s1) {  return (real_t)s1*lx/(real_t)ns1; }
  real_t _S2IDXtoY0(real_t s2) {  return (real_t)s2*ly/(real_t)ns2; }
  real_t _S3IDXtoZ0(real_t s3) {  return (real_t)s3*lz/(real_t)ns3; }

  void _MassDeposit(real_t weight, real_t x_idx, real_t y_idx, real_t z_idx,
                    arr_t &rho);

  void _CICDeposit(real_t weight, real_t x_idx, real_t y_idx, real_t z_idx,
                   arr_t &rho);

  void _PCSDeposit(real_t weight, real_t x_idx, real_t y_idx, real_t z_idx,
                   arr_t &rho);

  void _CINTDeposit(real_t weight, real_t x_idx, real_t y_idx, real_t z_idx,
                   arr_t &rho);

  void _deconvolve(arr_t &field);

  void _pushSheetMassToRho(idx_t s1, idx_t s2, idx_t s3);

  void _setMetricPotentials();

  void _RK4Calc();

  void _stepInit();

  void _K1Finalize();

  void _K2Finalize();

  void _K3Finalize();

  void _K4Finalize();

  real_t _getXRangeInSVoxel(register_t & DX, idx_t s1_idx, idx_t s2_idx,
                            idx_t s3_idx, real_t X0_lower, real_t X0_upper);

  void addBSSNSource(BSSN *bssn, real_t tot_mass);

  void rescaleFieldPerturbations(arr_t & field, real_t multiplier);
  void rescaleVelocityPerturbations(arr_t & ux, arr_t & uy, arr_t & uz, real_t multiplier);
  void rescalePositionPerturbations(arr_t & dx, arr_t & dy, arr_t & dz, real_t multiplier);
  void rescaleAllFieldPerturbations(BSSN *bssn, real_t multiplier);

  void RKStep(BSSN *bssn);

  void stepInit();

  void K1Finalize();
  void K2Finalize();
  void K3Finalize();
  void K4Finalize();
  
  arr_t d1alpha_a, d2alpha_a, d3alpha_a;

  arr_t d1gammai11_a, d1gammai22_a, d1gammai33_a,
    d1gammai12_a, d1gammai13_a, d1gammai23_a;

  arr_t d2gammai11_a, d2gammai22_a, d2gammai33_a,
    d2gammai12_a, d2gammai13_a, d2gammai23_a;

  arr_t d3gammai11_a, d3gammai22_a, d3gammai33_a,
    d3gammai12_a, d3gammai13_a, d3gammai23_a;

  arr_t d1beta1_a, d1beta2_a, d1beta3_a;
  arr_t d2beta1_a, d2beta2_a, d2beta3_a;
  arr_t d3beta1_a, d3beta2_a, d3beta3_a;

  std::vector<real_t> getgammaiIJ(idx_t s1, idx_t s2, idx_t s3, BSSN *bssnSim);
  std::vector<real_t> getgammaIJ(idx_t s1, idx_t s2, idx_t s3, BSSN *bssnSim);
  std::vector<real_t> getRayDataAtS(idx_t s, BSSN *bssnSim, Lambda * lambda);

};

real_t dot_cov_spatial_vectors(real_t * v1, real_t * v2, std::vector<real_t> gammaiIJ);
real_t mag_cov_spatial_vector(real_t * v1, std::vector<real_t> gammaiIJ);
real_t dot_cont_spatial_vectors(real_t * v1, real_t * v2, std::vector<real_t> gammaIJ);
real_t mag_cont_spatial_vector(real_t * v1, std::vector<real_t> gammaIJ);
real_t dot_4_vectors_vvg(real_t v1[4], real_t v2[4], real_t g[4][4]);
real_t dot_4_vectors_vv(real_t v1[4], real_t v2[4]);


} //namespace cosmo
#endif // include guard