Commit d06d9a71 authored by Giovanni La Mura's avatar Giovanni La Mura
Browse files

Merge branch 'fix_layered_spheres_bug' into 'master' 

Layered spheres bug fix

See merge request giacomo.mulas/np_tmcode!89
parents c8b6ff9b 47ee8b28

.gitlab-ci.yml

+5 −5
Original line number Diff line number Diff line
@@ -264,10 +264,10 @@ running_stage: 
      - mv c_TEDF c_TEDF_dev

      - mv c_TEDF.hd5 c_TEDF_dev.hd5

      - echo "Testing parallel execution of SPHERE"

      - ../../src/scripts/model_maker.py ../../test_data/sphere/config_181_xi.yml

      - ../../src/scripts/model_maker.py ../../test_data/sphere/config_two_layers.yml

      - echo "Running parallel np_sphere"

      - OMP_NUM_THREADS=2 ./np_sphere DEDFB_181_xi DSPH_181_xi .

      - mv c_OSPH c_OSPH_181_xi

      - OMP_NUM_THREADS=2 ./np_sphere DEDFB_two_layers DSPH_two_layers .

      - mv c_OSPH c_OSPH_two_layers

      - cd ../cluster

      - echo "Testing configurator for CLUSTER"

      - ../../src/scripts/model_maker.py ../../test_data/cluster/config_dev.yml
@@ -313,8 +313,8 @@ testing_stage: 
      - export FFILE=../../test_data/sphere/OSPH

      - python3 ../../src/scripts/pycompare.py --no-progress --ffile $FFILE --cfile c_OSPH_dev --html

      - echo "Comparing output of parallel SPHERE"

      - export FFILE=../../test_data/sphere/OSPH_181_xi

      - python3 ../../src/scripts/pycompare.py --no-progress --ffile $FFILE --cfile c_OSPH_181_xi --html

      - export FFILE=../../test_data/sphere/OSPH_two_layers

      - python3 ../../src/scripts/pycompare.py --no-progress --ffile $FFILE --cfile c_OSPH_two_layers --html

      - echo "Checking consistency among legacy and HDF5 configuration files"

      - ../testing/test_TEDF ../../test_data/sphere/DEDFB c_TEDF_dev c_TEDF_dev.hd5

      - cd ../inclusion

src/include/Configuration.h

+4 −0
Original line number Diff line number Diff line
@@ -312,6 +312,8 @@ protected: 
  double _wp;

  //! \brief Peak XI. QUESTION: correct?

  double _xip;

  //! \brief Maximum number of layers for the particle components.

  int _max_layers;

  //! \brief Flag to control whether to add an external layer.

  bool _use_external_sphere;
@@ -377,6 +379,8 @@ public: 
  const double& wp = _wp;

  //! \brief Read only view on peak XI.

  const double& xip = _xip;

  //! \brief Read only view on the maximum number of layers for the particle components.

  const int& max_layers = _max_layers;

  //! \brief Read only view on flag to control whether to add an external layer.

  const bool& use_external_sphere = _use_external_sphere;

src/libnptm/Configuration.cpp

+68 −57
Original line number Diff line number Diff line
@@ -356,10 +356,17 @@ ScattererConfiguration::ScattererConfiguration( 
  _iog_vec = iog_vector;

  _radii_of_spheres = ros_vector;

  _nshl_vec = nshl_vector;

  _use_external_sphere = is_external;

  _max_layers = 0;

  for (int li = 0; li < _configurations; li++) {

    if (_max_layers < _nshl_vec[li]) {

      _max_layers = _nshl_vec[li];

      if (li == 0 && _use_external_sphere) _max_layers += 1;

    }

  }

  _rcf = rcf_vector;

  _idfc = dielectric_func_type;

  _dc0_matrix = dc_matrix;

  _use_external_sphere = is_external;

  _exdc = ex;

  _wp = w;

  _xip = x;
@@ -391,12 +398,13 @@ ScattererConfiguration::ScattererConfiguration(const ScattererConfiguration& rhs 
  _exdc = rhs._exdc;

  _wp = rhs._wp;

  _xip = rhs._xip;

  _max_layers = rhs._max_layers;

  _iog_vec = new int[_number_of_spheres]();

  _radii_of_spheres = new double[_number_of_spheres]();

  _nshl_vec = new int[_configurations]();

  _rcf = new double*[_configurations];

  _scale_vec = new double[_number_of_scales]();

  _dc0_matrix = new dcomplex**[_configurations];

  _dc0_matrix = new dcomplex**[_max_layers];

  for (int si = 0; si < _number_of_scales; si++) _scale_vec[si] = rhs._scale_vec[si];

  for (int si = 0; si < _number_of_spheres; si++) {

    _iog_vec[si] = rhs._iog_vec[si];
@@ -406,11 +414,13 @@ ScattererConfiguration::ScattererConfiguration(const ScattererConfiguration& rhs 
    _radii_of_spheres[si] = rhs._radii_of_spheres[si];

    _nshl_vec[si] = rhs._nshl_vec[si];

    _rcf[si] = new double[_nshl_vec[si]]();

    _dc0_matrix[si] = new dcomplex*[_number_of_spheres];

    for (int sj = 0; sj < _nshl_vec[si]; sj++) _rcf[si][sj] = rhs._rcf[si][sj];

    for (int sj = 0; sj < _number_of_spheres; sj++) {

      _dc0_matrix[si][sj] = new dcomplex[dim3]();

      for (int sk = 0; sk < dim3; sk++) _dc0_matrix[si][sj][sk] = rhs._dc0_matrix[si][sj][sk];

  }

  for (int li = 0; li < _max_layers; li++) {

    _dc0_matrix[li] = new dcomplex*[_number_of_spheres];

    for (int lj = 0; lj < _number_of_spheres; lj++) {

      _dc0_matrix[li][lj] = new dcomplex[dim3]();

      for (int lk = 0; lk < dim3; lk++) _dc0_matrix[li][lj][lk] = rhs._dc0_matrix[li][lj][lk];

    }

  }

}
@@ -431,6 +441,7 @@ ScattererConfiguration::ScattererConfiguration(const mixMPI *mpidata) 
  itemp = sizeof(bool);

  char *ptemp = (char *) &_use_external_sphere;

  MPI_Bcast(ptemp, itemp, MPI_CHAR, 0, MPI_COMM_WORLD);

  MPI_Bcast(&_max_layers, 1, MPI_INT, 0, MPI_COMM_WORLD);

  MPI_Bcast(&_exdc, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);

  MPI_Bcast(&_wp, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);

  MPI_Bcast(&_xip, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
@@ -442,16 +453,18 @@ ScattererConfiguration::ScattererConfiguration(const mixMPI *mpidata) 
  MPI_Bcast(_nshl_vec, _configurations, MPI_INT, 0, MPI_COMM_WORLD);

  _rcf = new double*[_configurations];

  _scale_vec = new double[_number_of_scales]();

  _dc0_matrix = new dcomplex**[_configurations];

  _dc0_matrix = new dcomplex**[_max_layers];

  MPI_Bcast(_scale_vec, _number_of_scales, MPI_DOUBLE, 0, MPI_COMM_WORLD);

  int dim3 = (_idfc == 0) ? _number_of_scales : 1;

  for (int si = 0; si < _configurations; si++) {

    _rcf[si] = new double[_nshl_vec[si]]();

    MPI_Bcast(_rcf[si], _nshl_vec[si], MPI_DOUBLE, 0, MPI_COMM_WORLD);

    _dc0_matrix[si] = new dcomplex*[_number_of_spheres];

    for (int sj = 0; sj < _number_of_spheres; sj++) {

      _dc0_matrix[si][sj] = new dcomplex[dim3]();

      MPI_Bcast(_dc0_matrix[si][sj], dim3, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);

  }

  for (int li = 0; li < _max_layers; li++) {

    _dc0_matrix[li] = new dcomplex*[_number_of_spheres];

    for (int lj = 0; lj < _number_of_spheres; lj++) {

      _dc0_matrix[li][lj] = new dcomplex[dim3]();

      MPI_Bcast(_dc0_matrix[li][lj], dim3, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);

    }

  }

}
@@ -469,6 +482,7 @@ void ScattererConfiguration::mpibcast(const mixMPI *mpidata) { 
  itemp = sizeof(bool);

  char *ptemp = (char *) &_use_external_sphere;

  MPI_Bcast(ptemp, itemp, MPI_CHAR, 0, MPI_COMM_WORLD);

  MPI_Bcast(&_max_layers, 1, MPI_INT, 0, MPI_COMM_WORLD);

  MPI_Bcast(&_exdc, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);

  MPI_Bcast(&_wp, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);

  MPI_Bcast(&_xip, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
@@ -479,15 +493,17 @@ void ScattererConfiguration::mpibcast(const mixMPI *mpidata) { 
  int dim3 = (_idfc == 0) ? _number_of_scales : 1;

  for (int si = 0; si < _configurations; si++) {

    MPI_Bcast(_rcf[si], _nshl_vec[si], MPI_DOUBLE, 0, MPI_COMM_WORLD);

    for (int sj = 0; sj < _number_of_spheres; sj++) {

      MPI_Bcast(_dc0_matrix[si][sj], dim3, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);

  }

  for (int li = 0; li < _max_layers; li++) {

    for (int lj = 0; lj < _number_of_spheres; lj++) {

      MPI_Bcast(_dc0_matrix[li][lj], dim3, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);

    }

  }

}

#endif



ScattererConfiguration::~ScattererConfiguration() {

  for (int i = 0; i < _configurations; i++) {

  for (int i = 0; i < _max_layers; i++) {

    for (int j = 0; j < _number_of_spheres; j++) {

      delete[] _dc0_matrix[i][j];

    }
@@ -697,6 +713,7 @@ ScattererConfiguration* ScattererConfiguration::from_dedfb(const std::string& de 
  double *ros_vector = new double[configurations]();

  int *nshl_vector = new int[configurations]();

  double **rcf_vector = new double*[configurations];

  int max_layers = 0;

  last_configuration = 0;

  for (int i113 = 1; i113 <= fnsph; i113++) {

    if (iog_vector[i113 - 1] < i113) continue;
@@ -705,6 +722,9 @@ ScattererConfiguration* ScattererConfiguration::from_dedfb(const std::string& de 
    re = regex("[0-9]+");

    regex_search(str_target, m, re);

    nshl_vector[last_configuration - 1] = stoi(m.str());

    if (max_layers < nshl_vector[last_configuration - 1])

      max_layers = nshl_vector[last_configuration - 1];

    if (last_configuration == 1 && fies == 1) max_layers += 1;

    str_target = m.suffix().str();

    re = regex("-?[0-9]+\\.[0-9]+([eEdD][-+]?[0-9]+)?");

    regex_search(str_target, m, re);
@@ -725,10 +745,10 @@ ScattererConfiguration* ScattererConfiguration::from_dedfb(const std::string& de 
    } // ns112 loop

  } // i113 loop



  dcomplex ***dc0m = new dcomplex**[configurations];

  dcomplex *dc0 = new dcomplex[configurations]();

  dcomplex ***dc0m = new dcomplex**[max_layers];

  dcomplex *dc0 = new dcomplex[max_layers]();

  int dim3 = (fidfc == 0) ? fnxi : 1;

  for (int di = 0; di < configurations; di++) {

  for (int di = 0; di < max_layers; di++) {

    dc0m[di] = new dcomplex*[fnsph];

    for (int dj = 0; dj < fnsph; dj++) dc0m[di][dj] = new dcomplex[dim3]();

  } // di loop
@@ -763,21 +783,9 @@ ScattererConfiguration* ScattererConfiguration::from_dedfb(const std::string& de 
  delete[] dc0;



  ScattererConfiguration *config = new ScattererConfiguration(

							      fnsph,

							      configurations,

							      variable_vector,

							      fnxi,

							      variable_name,

							      iog_vector,

							      ros_vector,

							      nshl_vector,

							      rcf_vector,

							      fidfc,

							      dc0m,

							      (fies > 0),

							      fexdc,

							      fwp,

							      fxip

    fnsph, configurations, variable_vector, fnxi, variable_name,

    iog_vector, ros_vector, nshl_vector, rcf_vector, fidfc,

    dc0m, (fies > 0), fexdc, fwp, fxip

  );

  delete[] file_lines;

  delete[] variable_vector;
@@ -820,12 +828,17 @@ ScattererConfiguration* ScattererConfiguration::from_hdf5(const std::string& fil 
    ros_vector = new double[configuration_count]();

    rcf_vector = new double*[configuration_count];

    last_configuration = 0;

    int max_layers = 0;

    for (int i115 = 1; i115 <= configuration_count; i115++) {

      if (iog[i115 - 1] < i115) continue;

      else last_configuration++;

      str_name = "NSHL_" + to_string(last_configuration);

      str_type = "INT32_(1)";

      status = hdf_file->read(str_name, str_type, (nshl_vector + last_configuration - 1));

      if (max_layers < nshl_vector[last_configuration - 1]) {

	max_layers = nshl_vector[last_configuration - 1];

	if (last_configuration == 1 && ies == 1) max_layers += 1;

      }

      str_name = "ROS_" + to_string(last_configuration);

      str_type = "FLOAT64_(1)";

      status = hdf_file->read(str_name, str_type, (ros_vector + last_configuration - 1));
@@ -841,14 +854,14 @@ ScattererConfiguration* ScattererConfiguration::from_hdf5(const std::string& fil 
    status = hdf_file->read(str_name, str_type, xi_vec);

      

    int dim3 = (_idfc == 0) ? nxi : 1;

    int element_size = 2 * dim3 * nsph * configuration_count;

    int element_size = 2 * dim3 * nsph * max_layers;

    double *elements = new double[element_size]();

    str_name = "DC0M";

    str_type = "FLOAT64_(" + to_string(element_size) + ")";

    status = hdf_file->read(str_name, str_type, elements);

    dc0m = new dcomplex**[configuration_count];

    dc0m = new dcomplex**[max_layers];

    int dc_index = 0;

    for (int di = 0; di < configuration_count; di++) {

    for (int di = 0; di < max_layers; di++) {

      dc0m[di] = new dcomplex*[nsph];

      for (int dj = 0; dj < nsph; dj++) {

	dc0m[di][dj] = new dcomplex[dim3]();
@@ -864,21 +877,9 @@ ScattererConfiguration* ScattererConfiguration::from_hdf5(const std::string& fil 
    status = hdf_file->close();

    delete hdf_file;

    conf = new ScattererConfiguration(

				      nsph,

				      configuration_count,

				      xi_vec,

				      nxi,

				      "XIV",

				      iog,

				      ros_vector,

				      nshl_vector,

				      rcf_vector,

				      _idfc,

				      dc0m,

				      (ies == 1),

				      _exdc,

				      _wp,

				      _xip

      nsph, configuration_count, xi_vec, nxi, "XIV",

      iog, ros_vector, nshl_vector, rcf_vector, _idfc,

      dc0m, (ies == 1), _exdc, _wp, _xip

    );

    delete[] xi_vec;

  }
@@ -932,9 +933,14 @@ ScattererConfiguration* ScattererConfiguration::from_legacy(const std::string& f 
    for (int nsi = 0; nsi < nsh; nsi++)

      input.read(reinterpret_cast<char *>(&(_rcf_vec[last_configuration - 1][nsi])), sizeof(double));

  }

  _dc0m = new dcomplex**[configurations];

  int max_layers = 0;

  for (int li = 0; li < configurations; li++) {

    if (max_layers < _nshl_vec[li]) max_layers = _nshl_vec[li];

    if (li == 0 && _ies == 1) max_layers += 1;

  }

  _dc0m = new dcomplex**[max_layers];

  int dim3 = (_idfc == 0) ? _nxi : 1;

  for (int di = 0; di < configurations; di++) {

  for (int di = 0; di < max_layers; di++) {

    _dc0m[di] = new dcomplex*[_nsph];

    for (int dj = 0; dj < _nsph; dj++) _dc0m[di][dj] = new dcomplex[dim3]();

  } // di loop
@@ -995,7 +1001,7 @@ void ScattererConfiguration::print() { 
  for (int i = 0; i < _number_of_scales; i++) printf("\t%lg", _scale_vec[i]);

  printf("\t]\n");

  printf("DC0M  = [\n");

  for (int i = 0; i < _configurations; i++) {

  for (int i = 0; i < _max_layers; i++) {

    printf("         [\n");

    for (int j = 0; j < _number_of_spheres; j++) {

      printf("          [");
@@ -1079,7 +1085,7 @@ void ScattererConfiguration::write_hdf5(const std::string& file_name) { 
  }



  int dim3 = (idfc == 0) ? _number_of_scales : 1;

  int dc0m_size = 2 * dim3 * _number_of_spheres * _configurations;

  int dc0m_size = 2 * dim3 * _number_of_spheres * _max_layers;

  double *dc0m = new double[dc0m_size]();

  int dc0_index = 0;

  for (int jxi468 = 1; jxi468 <= _number_of_scales; jxi468++) {
@@ -1376,6 +1382,9 @@ bool ScattererConfiguration::operator ==(const ScattererConfiguration &other) { 
  if (_idfc != other._idfc) {

    return false;

  }

  if (_max_layers != other._max_layers) {

    return false;

  }

  if (_exdc != other._exdc) {

    return false;

  }
@@ -1411,13 +1420,15 @@ bool ScattererConfiguration::operator ==(const ScattererConfiguration &other) { 
	return false;

      }

    } // rj loop

  } // ri loop

  for (int li = 0; li < _max_layers; li++) {

    for (int dj = 0; dj < _number_of_spheres; dj++) {

      for (int dk = 0; dk < dim3; dk++) {

	if (_dc0_matrix[ri][dj][dk] != other._dc0_matrix[ri][dj][dk]) {

	if (_dc0_matrix[li][dj][dk] != other._dc0_matrix[li][dj][dk]) {

	  return false;

	}

      } // dk loop

    } // dj loop

  } // ri loop

  }

  return true;

}

src/scripts/config_gen.yml

0 → 100644
+93 −0
Original line number Diff line number Diff line 
system_settings:

  max_host_ram : 0

  max_gpu_ram  : 0



input_settings:

  input_folder : "."

  spheres_file : "DEDFB"

  geometry_file: "DCLU"

  

output_settings:

  output_folder: "."

  output_name  : "c_OCLU"

  formats      : [ "LEGACY", "HDF5" ]

  jwtm         : 1



particle_settings:

  application : "CLUSTER"

  n_spheres   : 8

  n_types     : 2

  sph_types   : [ 1, 1, 1, 1, 2, 2, 2, 2 ]

  n_layers    : [ 1, 1 ]

  radii       : [ 6.000e-08, 4.000e-8 ]

  rad_frac    : [ [ 1.0 ], [ 1.0 ] ]

  dielec_id   : [ [ 1 ], [ 1 ] ]



material_settings:

  diel_flag   : 0

  extern_diel : 1.00e+00

  dielec_path : "../../ref_data"

  dielec_file : [ "eps_ashok_C.csv" ]

  dielec_fmt  : [ "CSV" ]

  match_mode  : "GRID"

  diel_const  : [ ]



radiation_settings:

  polarization: "LINEAR"

  scale_start : 1.00e-06

  scale_end   : 2.00e-05

  scale_step  : 5.00e-09

  wp          : 2.99792e+08

  xip         : 1.00e+00

  step_flag   : 0

  scale_name  : "WAVELENGTH"



geometry_settings:

  li          : 4

  le          : 8

  npnt        : 149

  npntts      : 300

  iavm        : 0

  isam        : 0

  in_th_start : 0.0

  in_th_step  : 0.0

  in_th_end   : 0.0

  in_ph_start : 0.0

  in_ph_step  : 0.0

  in_ph_end   : 0.0

  sc_th_start : 0.0

  sc_th_step  : 0.0

  sc_th_end   : 0.0

  sc_ph_start : 0.0

  sc_ph_step  : 0.0

  sc_ph_end   : 0.0

  x_coords    : [

    8.00e-08,

    -8.00e-08,

    0.00e+00,

    0.00e+00,

    0.00e+00,

    0.00e+00,

    0.00e+00,

    0.00e+00

    ]

  y_coords    : [

    0.00e+00,

    0.00e+00,

    8.00e-08,

    -8.00e-08,

    0.00e+00,

    0.00e+00,

    0.00e+00,

    0.00e+00

    ]

  z_coords    : [

    0.00e+00,

    0.00e+00,

    0.00e+00,

    0.00e+00,

    8.00e-08,

    -8.00e-08,

    16.00e-08,

    -16.00e-08

    ]

src/scripts/generator.py

0 → 100644
+95 −0
Original line number Diff line number Diff line 
import math

import random

import yaml



#import pdb



from pathlib import PurePath

from sys import argv



seed = int(argv[1])



random.seed(seed)



config = {}

with open('config.yml', 'r') as stream:

    config = yaml.safe_load(stream)



nsph = int(config['particle_settings']['n_spheres'])



vec_thetas = [0.0 for i in range(nsph)]

vec_phis = [0.0 for i in range(nsph)]

vec_rads = [0.0 for i in range(nsph)]

vec_sph_x = [0.0 for i in range(nsph)]

vec_sph_y = [0.0 for i in range(nsph)]

vec_sph_z = [0.0 for i in range(nsph)]



sph_type_index = config['particle_settings']['sph_types'][0] - 1

vec_rads[0] = config['particle_settings']['radii'][sph_type_index]

max_rad = 20.0 * vec_rads[0]

placed_spheres = 1

attempts = 0

max_attempts = 100

for i in range(1, nsph):

    sph_type_index = config['particle_settings']['sph_types'][i] - 1

    vec_rads[i] = config['particle_settings']['radii'][sph_type_index]

    is_placed = False

    #breakpoint()

    while (not is_placed):

        if (attempts > max_attempts):

            print("WARNING: could not place sphere %d in allowed radius!"%i)

            break # while(not is_placed)

        vec_thetas[i] = math.pi * random.random()

        vec_phis[i] = 2.0 * math.pi * random.random()

        rho = vec_rads[0] + vec_rads[i]

        z = rho * math.cos(vec_thetas[i])

        y = rho * math.sin(vec_thetas[i]) * math.sin(vec_phis[i])

        x = rho * math.sin(vec_thetas[i]) * math.cos(vec_phis[i])

        j = 0

        while (j < i - 1):

            j += 1

            dx2 = (x - vec_sph_x[j]) * (x - vec_sph_x[j])

            dy2 = (y - vec_sph_y[j]) * (y - vec_sph_y[j])

            dz2 = (z - vec_sph_z[j]) * (z - vec_sph_z[j])

            dist2 = dx2 + dy2 + dz2

            rr2 = (vec_rads[i] + vec_rads[j]) * (vec_rads[i] + vec_rads[j])

            if (dist2 < rr2):

                # Spheres i and j are compenetrating.

                # Sphere i is moved out radially until it becomes externally

                # tangent to sphere j. Then the check is repeated, to verify

                # that no other sphere was penetrated. The process is iterated

                # until sphere i is placed or the maximum allowed radius is

                # reached.

                sinthi = math.sin(vec_thetas[i])

                sinthj = math.sin(vec_thetas[j])

                costhi = math.cos(vec_thetas[i])

                costhj = math.cos(vec_thetas[j])

                sinphi = math.sin(vec_phis[i])

                sinphj = math.sin(vec_phis[j])

                cosphi = math.cos(vec_phis[i])

                cosphj = math.cos(vec_phis[j])

                cosalpha = (

                    sinthi * cosphi * sinthj * cosphj

                    + sinthi * sinphi * sinthj * sinphj

                    + costhi * costhj

                )

                rho += 2.0 * vec_rads[j] * cosalpha

                z = rho * math.cos(vec_thetas[i])

                y = rho * math.sin(vec_thetas[i]) * math.sin(vec_phis[i])

                x = rho * math.sin(vec_thetas[i]) * math.cos(vec_phis[i])

                j = 0

                continue # while(j < i - 1)

        if (rho + vec_rads[i] > max_rad):

            # The current direction is filled. Try another one.

            attempts += 1

            continue # while(not is_placed)

        vec_sph_x[i] = x

        vec_sph_y[i] = y

        vec_sph_z[i] = z

        is_placed = True

        placed_spheres += 1

        attempts = 0



print(vec_sph_x)

print(vec_sph_y)

print(vec_sph_z)