/* Copyright 2024 INAF - Osservatorio Astronomico di Cagliari

   Licensed under the Apache License, Version 2.0 (the "License");
   you may not use this file except in compliance with the License.
   You may obtain a copy of the License at

       http://www.apache.org/licenses/LICENSE-2.0

   Unless required by applicable law or agreed to in writing, software
   distributed under the License is distributed on an "AS IS" BASIS,
   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
   See the License for the specific language governing permissions and
   limitations under the License.
 */

/*! \file Commons.cpp
 *
 * \brief Implementation of the common data structures.
 */
#ifndef INCLUDE_TYPES_H_
#include "../include/types.h"
#endif

#ifndef INCLUDE_CONFIGURATION_H_
#include "../include/Configuration.h"
#endif

#ifndef INCLUDE_COMMONS_H
#include "../include/Commons.h"
#endif

#ifdef USE_MPI
#include <mpi.h>
#endif

C1::C1(GeometryConfiguration *gconf, ScattererConfiguration *sconf) {
  lm = gconf->l_max;
  int li = gconf->li;
  int le = gconf->le;
  if (lm == 0) {
    lm = (li > le) ? li : le;
  }
  nsph = gconf->number_of_spheres;
  nlmmt = 2 * (lm * (lm + 2));

  vec_rmi = new dcomplex[nsph * lm]();
  vec_rei = new dcomplex[nsph * lm]();
  rmi = new dcomplex*[lm];
  rei = new dcomplex*[lm];
  for (int ri = 0; ri < lm; ri++) {
    rmi[ri] = &(vec_rmi[nsph * ri]);
    rei[ri] = &(vec_rei[nsph * ri]);
  }
  vec_w = new dcomplex[4 * nlmmt]();
  w = new dcomplex*[nlmmt];
  for (int wi = 0; wi < nlmmt; wi++) w[wi] = &(vec_w[4 * wi]);
  configurations = sconf->configurations;
  vint = new dcomplex[16]();
  vec_vints = new dcomplex[nsph * 16]();
  vints = new dcomplex*[nsph];
  rc = new double*[configurations];
  nshl = new int[configurations]();
  iog = new int[nsph]();
  int conf_index = 0;
  for (int vi = 0; vi < nsph; vi++) {
    vints[vi] = &(vec_vints[vi * 16]);
    iog[vi] = sconf->get_iog(vi);
    if (iog[vi] >= vi + 1) {
      nshl[conf_index] = sconf->get_nshl(conf_index);
      rc[conf_index] = new double[nshl[conf_index]]();
      conf_index++;
    }
  }
  fsas = new dcomplex[nsph]();
  sscs = new double[nsph]();
  sexs = new double[nsph]();
  sabs = new double[nsph]();
  sqscs = new double[nsph]();
  sqexs = new double[nsph]();
  sqabs = new double[nsph]();
  gcsv = new double[nsph]();
  rxx = new double[nsph]();
  ryy = new double[nsph]();
  rzz = new double[nsph]();
  if (nsph > 1) {
    for (int c1i = 0; c1i < nsph; c1i++) {
      rxx[c1i] = gconf->get_sph_x(c1i);
      ryy[c1i] = gconf->get_sph_y(c1i);
      rzz[c1i] = gconf->get_sph_z(c1i);
    }
  }
  ros = new double[configurations]();
  for (int ri = 0; ri < configurations; ri++) {
    ros[ri] = sconf->get_radius(ri);
  }

  sas = new dcomplex**[nsph];
  sas[0] = new dcomplex*[nsph*2];
  sas[0][0] = new dcomplex[nsph*2*2];
  sas[0][1] = sas[0][0]+2;
  for (int si = 1; si < nsph; si++) {
    sas[si] = sas[0]+2*si;
    sas[si][0] = sas[0][0]+4*si;
    sas[si][1] = sas[si][0]+2;
  }
}

C1::C1(const C1& rhs) {
  nlmmt = rhs.nlmmt;
  nsph = rhs.nsph;
  lm = rhs.lm;

  vec_rmi = new dcomplex[nsph * lm]();
  vec_rei = new dcomplex[nsph * lm]();
  for (long li=0; li<(nsph*lm); li++) {
    vec_rmi[li] = rhs.vec_rmi[li];
    vec_rei[li] = rhs.vec_rei[li];
  }
  rmi = new dcomplex*[lm];
  rei = new dcomplex*[lm];
  for (int ri = 0; ri < lm; ri++) {
    rmi[ri] = &(vec_rmi[nsph * ri]);
    rei[ri] = &(vec_rei[nsph * ri]);
    /*! The contents were already copied via vec_rmi and vec_rei */
  }
  vec_w = new dcomplex[4 * nlmmt]();
  for (long li=0; li<(4*nlmmt); li++) {
    vec_w[li] = rhs.vec_w[li];
  }    
  w = new dcomplex*[nlmmt];
  for (int wi = 0; wi < nlmmt; wi++) {
    w[wi] = &(vec_w[4 * wi]);
    /*! The contents were already copied via vec_w */
  }
  configurations = rhs.configurations;
  vint = new dcomplex[16]();
  for (int vi = 0; vi < 16; vi++) vint[vi] = rhs.vint[vi];
  vec_vints = new dcomplex[nsph * 16]();
  vints = new dcomplex*[nsph];
  rc = new double*[nsph];
  nshl = new int[nsph]();
  for (int ni = 0; ni < nsph; ni++) nshl[ni] = rhs.nshl[ni];
  iog = new int[nsph]();
  for (int ni = 0; ni < nsph; ni++) iog[ni] = rhs.iog[ni];
  int conf_index = 0;
  for (int vi = 0; vi < nsph; vi++) {
    vints[vi] = &(vec_vints[16 * vi]);
    for (int vj = 0; vj < 16; vj++) vints[vi][vj] = rhs.vints[vi][vj];
  }
  for (int ri=0; ri<configurations; ri++) {
    rc[ri] = new double[nshl[ri]]();
    for (int rj=0; rj<nshl[ri]; rj++) rc[ri][rj] = rhs.rc[ri][rj];
  }
  fsas = new dcomplex[nsph]();
  sscs = new double[nsph]();
  sexs = new double[nsph]();
  sabs = new double[nsph]();
  sqscs = new double[nsph]();
  sqexs = new double[nsph]();
  sqabs = new double[nsph]();
  gcsv = new double[nsph]();
  rxx = new double[nsph]();
  ryy = new double[nsph]();
  rzz = new double[nsph]();
  ros = new double[nsph]();

  sas = new dcomplex**[nsph];
  sas[0] = new dcomplex*[nsph*2];
  sas[0][0] = new dcomplex[nsph*2*2];
  sas[0][1] = sas[0][0]+2;
  for (int si = 1; si < nsph; si++) {
    sas[si] = sas[0]+2*si;
    sas[si][0] = sas[0][0]+4*si;
    sas[si][1] = sas[si][0]+2;
  }
  for (int si = 0; si < nsph; si++) {
    for (int sj=0; sj<2; sj++) {
      for (int sk=0; sk<2; sk++) sas[si][sj][sk] = rhs.sas[si][sj][sk];
    }
    fsas[si] = rhs.fsas[si];
    sscs[si] = rhs.sscs[si];
    sexs[si] = rhs.sexs[si];
    sabs[si] = rhs.sabs[si];
    sqscs[si] = rhs.sqscs[si];
    sqexs[si] = rhs.sqexs[si];
    sqabs[si] = rhs.sqabs[si];
    gcsv[si] = rhs.gcsv[si];
    rxx[si] = rhs.rxx[si];
    ryy[si] = rhs.ryy[si];
    rzz[si] = rhs.rzz[si];
  }
  for (int si = 0; si < configurations; si++) {
    ros[si] = rhs.ros[si];    
  }
}

#ifdef MPI_VERSION
C1::C1(const mixMPI *mpidata) {
  MPI_Bcast(&nlmmt, 1, MPI_INT, 0, MPI_COMM_WORLD);
  MPI_Bcast(&nsph, 1, MPI_INT, 0, MPI_COMM_WORLD);
  MPI_Bcast(&lm, 1, MPI_INT, 0, MPI_COMM_WORLD);
  vec_rmi = new dcomplex[nsph * lm]();
  vec_rei = new dcomplex[nsph * lm]();
  MPI_Bcast(vec_rmi, nsph*lm, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
  MPI_Bcast(vec_rei, nsph*lm, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
  rmi = new dcomplex*[lm];
  rei = new dcomplex*[lm];
  for (int ri = 0; ri < lm; ri++) {
    rmi[ri] = &(vec_rmi[nsph * ri]);
    rei[ri] = &(vec_rei[nsph * ri]);
    /*! The contents were copied already via vec_rmi and vec_rei */
  }
  vec_w = new dcomplex[4 * nlmmt]();
  MPI_Bcast(vec_w, 4*nlmmt, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
  w = new dcomplex*[nlmmt];
  for (int wi = 0; wi < nlmmt; wi++) {
    w[wi] = &(vec_w[4 * wi]);
    /*! The contents were copied already via vec_w */
  }
  MPI_Bcast(&configurations, 1, MPI_INT, 0, MPI_COMM_WORLD);
  vint = new dcomplex[16]();
  MPI_Bcast(vint, 16, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
  vec_vints = new dcomplex[nsph * 16]();
  MPI_Bcast(vec_vints, nsph*16, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
  vints = new dcomplex*[nsph];
  rc = new double*[nsph];
  nshl = new int[nsph]();
  MPI_Bcast(nshl, nsph, MPI_INT, 0, MPI_COMM_WORLD);
  iog = new int[nsph]();
  MPI_Bcast(iog, nsph, MPI_INT, 0, MPI_COMM_WORLD);
  for (int vi = 0; vi < nsph; vi++) {
    vints[vi] = &(vec_vints[16 * vi]);
    // The contents were copied already via vec_vints
  }
  for (int ri=0; ri<configurations; ri++) {
    rc[ri] = new double[nshl[ri]]();
    MPI_Bcast(rc[ri], nshl[ri], MPI_DOUBLE, 0, MPI_COMM_WORLD);
  }
  fsas = new dcomplex[nsph]();
  sscs = new double[nsph]();
  sexs = new double[nsph]();
  sabs = new double[nsph]();
  sqscs = new double[nsph]();
  sqexs = new double[nsph]();
  sqabs = new double[nsph]();
  gcsv = new double[nsph]();
  rxx = new double[nsph]();
  ryy = new double[nsph]();
  rzz = new double[nsph]();
  ros = new double[nsph]();
  MPI_Bcast(fsas, nsph, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
  MPI_Bcast(sscs, nsph, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(sexs, nsph, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(sabs, nsph, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(sqscs, nsph, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(sqexs, nsph, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(sqabs, nsph, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(gcsv, nsph, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(rxx, nsph, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(ryy, nsph, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(rzz, nsph, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(ros, configurations, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  sas = new dcomplex**[nsph];
  sas[0] = new dcomplex*[nsph*2];
  sas[0][0] = new dcomplex[nsph*2*2];
  sas[0][1] = sas[0][0]+2;
  for (int si = 1; si < nsph; si++) {
    sas[si] = sas[0]+2*si;
    sas[si][0] = sas[0][0]+4*si;
    sas[si][1] = sas[si][0]+2;
  }
  for (int si = 0; si < nsph; si++) {
    for (int sj=0; sj<2; sj++) {
      MPI_Bcast(sas[si][sj], 2, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
    }
  }
}

void C1::mpibcast(const mixMPI *mpidata) {
  MPI_Bcast(&nlmmt, 1, MPI_INT, 0, MPI_COMM_WORLD);
  MPI_Bcast(&nsph, 1, MPI_INT, 0, MPI_COMM_WORLD);
  MPI_Bcast(&lm, 1, MPI_INT, 0, MPI_COMM_WORLD);
  MPI_Bcast(vec_rmi, nsph*lm, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
  MPI_Bcast(vec_rei, nsph*lm, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
  MPI_Bcast(vec_w, 4*nlmmt, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
  MPI_Bcast(&configurations, 1, MPI_INT, 0, MPI_COMM_WORLD);
  MPI_Bcast(vint, 16, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
  MPI_Bcast(vec_vints, nsph*16, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
  MPI_Bcast(nshl, nsph, MPI_INT, 0, MPI_COMM_WORLD);
  MPI_Bcast(iog, nsph, MPI_INT, 0, MPI_COMM_WORLD);
  for (int ri=0; ri<configurations; ri++) {
    MPI_Bcast(rc[ri], nshl[ri], MPI_DOUBLE, 0, MPI_COMM_WORLD);
  }
  MPI_Bcast(fsas, nsph, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
  MPI_Bcast(sscs, nsph, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(sexs, nsph, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(sabs, nsph, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(sqscs, nsph, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(sqexs, nsph, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(sqabs, nsph, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(gcsv, nsph, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(rxx, nsph, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(ryy, nsph, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(rzz, nsph, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(ros, configurations, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  for (int si = 0; si < nsph; si++) {
    for (int sj=0; sj<2; sj++) {
      MPI_Bcast(sas[si][sj], 2, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
    }
  }
}
#endif

C1::~C1() {
  delete[] vec_rmi;
  delete[] vec_rei;
  delete[] rmi;
  delete[] rei;
  delete[] vec_w;
  delete[] w;
  int conf_index = 0;
  for (int ci = 0; ci < configurations; ci++) {
    delete[] rc[ci];
  }
  delete[] rc;
  delete[] vint;
  delete[] vec_vints;
  delete[] vints;
  delete[] sas[0][0];
  delete[] sas[0];
  delete[] sas;
  delete[] fsas;
  delete[] sscs;
  delete[] sexs;
  delete[] sabs;
  delete[] sqscs;
  delete[] sqexs;
  delete[] sqabs;
  delete[] gcsv;
  delete[] rxx;
  delete[] ryy;
  delete[] rzz;
  delete[] ros;
  delete[] iog;
  delete[] nshl;
}

C1_AddOns::C1_AddOns(C4 *c4) {
  nsph = c4->nsph;
  lmpo = c4->lmpo;
  nlemt = 2 * c4->nlem;
  litpo = c4->litpo;
  lmtpo = c4->lmtpo;
  litpos = c4->litpos;
  lmtpos = c4->lmtpos;
  nv3j = c4->nv3j;
  lm = c4->lm;
  vh = new dcomplex[(nsph * nsph - 1) * litpo]();
  vj0 = new dcomplex[nsph * lmtpo]();
  vj = new dcomplex[1]();
  vyhj = new dcomplex[(nsph * nsph - 1) * litpos]();
  vyj0 = new dcomplex[nsph * lmtpos]();
  am0v = new dcomplex[nlemt * nlemt]();
  am0m = new dcomplex*[nlemt];
  for (int ai = 0; ai < nlemt; ai++) {
    am0m[ai] = (am0v + nlemt * ai);
  }
  vintm = new dcomplex[16]();
  vintt = new dcomplex[16]();
  // This calculates the size of v3j0
  const int nv3j = c4->nv3j;
  v3j0 = new double[nv3j]();
  ind3j = new int*[lmpo];
  for (int ii = 0; ii < lmpo; ii++) ind3j[ii] = new int[c4->lm]();
  fsac = new dcomplex*[2];
  sac = new dcomplex*[2];
  fsacm = new dcomplex*[2];
  for (int fi = 0; fi < 2; fi++) {
    fsac[fi] = new dcomplex[2]();
    sac[fi] = new dcomplex[2]();
    fsacm[fi] = new dcomplex[2]();
  }
  scscp = new dcomplex[2]();
  ecscp = new dcomplex[2]();
  scscpm = new dcomplex[2]();
  ecscpm = new dcomplex[2]();
  sscs = new double[nsph]();
  ecscm = new double[2]();
  scscm = new double[2]();
  scsc = new double[2]();
  ecsc = new double[2]();
}

C1_AddOns::C1_AddOns(const C1_AddOns& rhs) {
  nsph = rhs.nsph;
  lmpo = rhs.lmpo;
  nlemt = rhs.nlemt;
  litpo = rhs.litpo;
  lmtpo = rhs.lmtpo;
  litpos = rhs.litpos;
  lmtpos = rhs.lmtpos;
  nv3j = rhs.nv3j;
  lm = rhs.lm;
  int vhsize=(nsph * nsph - 1) * litpo;
  vh = new dcomplex[vhsize]();
  for (int hi=0; hi<vhsize; hi++) vh[hi] = rhs.vh[hi];
  int vj0size=nsph * lmtpo;
  vj0 = new dcomplex[vj0size]();
  for (int hi=0; hi<vj0size; hi++) vj0[hi] = rhs.vj0[hi];
  vj = new dcomplex[1](); // QUESTION: is 1 really enough for a general case?
  vj[0] = rhs.vj[0];
  int vyhjsize=(nsph * nsph - 1) * litpos;
  vyhj = new dcomplex[vyhjsize]();
  for (int hi=0; hi<vyhjsize; hi++) vyhj[hi] = rhs.vyhj[hi];
  int vyj0size = nsph * lmtpos;
  vyj0 = new dcomplex[vyj0size]();
  for (int hi=0; hi<vyj0size; hi++) vyj0[hi] = rhs.vyj0[hi];
  am0v = new dcomplex[nlemt * nlemt]();
  am0m = new dcomplex*[nlemt];
  for (int ai = 0; ai < nlemt; ai++) {
    for (int aj = 0; aj < nlemt; aj++) am0v[nlemt * ai + aj] = rhs.am0v[nlemt * ai + aj];
    am0m[ai] = (am0v + nlemt * ai);
  }
  vintm = new dcomplex[16]();
  vintt = new dcomplex[16]();
  for (int hi = 0; hi < 16; hi++) {
    vintm[hi] = rhs.vintm[hi];
    vintt[hi] = rhs.vintt[hi];
  }
  fsac = new dcomplex*[2];
  sac = new dcomplex*[2];
  fsacm = new dcomplex*[2];
  for (int fi = 0; fi < 2; fi++) {
    fsac[fi] = new dcomplex[2]();
    sac[fi] = new dcomplex[2]();
    fsacm[fi] = new dcomplex[2]();
    for (int fj = 0; fj < 2; fj++) {
      fsac[fi][fj] = rhs.fsac[fi][fj];
      sac[fi][fj] = rhs.sac[fi][fj];
      fsacm[fi][fj] = rhs.fsacm[fi][fj];
    }
  }
  scscp = new dcomplex[2]();
  ecscp = new dcomplex[2]();
  scscpm = new dcomplex[2]();
  ecscpm = new dcomplex[2]();
  v3j0 = new double[nv3j]();
  for (int hi = 0; hi < nv3j; hi++) v3j0[hi] = rhs.v3j0[hi];
  ind3j = new int*[lmpo];
  for (int ii = 0; ii < lmpo; ii++) {
    ind3j[ii] = new int[lm]();
    for (int ij=0; ij<lm; ij++) ind3j[ii][ij] = rhs.ind3j[ii][ij];
  }
  sscs = new double[nsph]();
  for (int hi=0; hi<nsph; hi++) sscs[hi] = rhs.sscs[hi];
  ecscm = new double[2]();
  scscm = new double[2]();
  scsc = new double[2]();
  ecsc = new double[2]();
  for (int hi=0; hi<2; hi++) {
    scscp[hi] = rhs.scscp[hi];
    ecscp[hi] = rhs.ecscp[hi];
    scscpm[hi] = rhs.scscpm[hi];
    ecscpm[hi] = rhs.ecscpm[hi];
    ecscm[hi] = rhs.ecscm[hi];
    scscm[hi] = rhs.scscm[hi];
    scsc[hi] = rhs.scsc[hi];
    ecsc[hi] = rhs.ecsc[hi];
  }
}

C1_AddOns::~C1_AddOns() {
  delete[] sscs;
  delete[] vyj0;
  delete[] vyhj;
  for (int ii = lmpo - 1; ii > -1; ii--) delete[] ind3j[ii];
  delete[] ind3j;
  delete[] v3j0;
  delete[] vh;
  delete[] vj0;
  delete[] vj;
  delete[] am0m;
  delete[] am0v;
  delete[] vintm;
  delete[] vintt;
  for (int fi = 1; fi > -1; fi--) {
    delete[] fsac[fi];
    delete[] sac[fi];
    delete[] fsacm[fi];
  }
  delete[] fsac;
  delete[] sac;
  delete[] fsacm;
  delete[] scscp;
  delete[] ecscp;
  delete[] scscpm;
  delete[] ecscpm;
  delete[] ecscm;
  delete[] scscm;
  delete[] scsc;
  delete[] ecsc;
}

#ifdef MPI_VERSION
C1_AddOns::C1_AddOns(const mixMPI *mpidata) {
  MPI_Bcast(&nsph, 1, MPI_INT, 0, MPI_COMM_WORLD);
  MPI_Bcast(&lmpo, 1, MPI_INT, 0, MPI_COMM_WORLD);
  MPI_Bcast(&nlemt, 1, MPI_INT, 0, MPI_COMM_WORLD);
  MPI_Bcast(&litpo, 1, MPI_INT, 0, MPI_COMM_WORLD);
  MPI_Bcast(&lmtpo, 1, MPI_INT, 0, MPI_COMM_WORLD);
  MPI_Bcast(&litpos, 1, MPI_INT, 0, MPI_COMM_WORLD);
  MPI_Bcast(&lmtpos, 1, MPI_INT, 0, MPI_COMM_WORLD);
  MPI_Bcast(&nv3j, 1, MPI_INT, 0, MPI_COMM_WORLD);
  MPI_Bcast(&lm, 1, MPI_INT, 0, MPI_COMM_WORLD);
  int vhsize=(nsph * nsph - 1) * litpo;
  vh = new dcomplex[vhsize]();
  MPI_Bcast(vh, vhsize, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
  int vj0size=nsph * lmtpo;
  vj0 = new dcomplex[vj0size]();
  MPI_Bcast(vj0, vj0size, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
  vj = new dcomplex[1](); // QUESTION: is 1 really enough for a general case?
  MPI_Bcast(vj, 1, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
  int vyhjsize=(nsph * nsph - 1) * litpos;
  vyhj = new dcomplex[vyhjsize]();
  MPI_Bcast(vyhj, vyhjsize, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
  int vyj0size = nsph * lmtpos;
  vyj0 = new dcomplex[vyj0size]();
  MPI_Bcast(vyj0, vyj0size, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
  am0v = new dcomplex[nlemt * nlemt]();
  am0m = new dcomplex*[nlemt];
  MPI_Bcast(am0v, nlemt*nlemt, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
  for (int ai = 0; ai < nlemt; ai++) {
    am0m[ai] = (am0v + nlemt * ai);
  }
  vintm = new dcomplex[16]();
  vintt = new dcomplex[16]();
  MPI_Bcast(vintm, 16, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
  MPI_Bcast(vintt, 16, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
  fsac = new dcomplex*[2];
  sac = new dcomplex*[2];
  fsacm = new dcomplex*[2];
  for (int fi = 0; fi < 2; fi++) {
    fsac[fi] = new dcomplex[2]();
    sac[fi] = new dcomplex[2]();
    fsacm[fi] = new dcomplex[2]();
    MPI_Bcast(fsac[fi], 2, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
    MPI_Bcast(sac[fi], 2, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
    MPI_Bcast(fsacm[fi], 2, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
  }
  scscp = new dcomplex[2]();
  ecscp = new dcomplex[2]();
  scscpm = new dcomplex[2]();
  ecscpm = new dcomplex[2]();
  MPI_Bcast(scscp, 2, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
  MPI_Bcast(ecscp, 2, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
  MPI_Bcast(scscpm, 2, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
  MPI_Bcast(ecscpm, 2, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
  v3j0 = new double[nv3j]();
  MPI_Bcast(v3j0, nv3j, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  ind3j = new int*[lmpo];
  for (int ii = 0; ii < lmpo; ii++) {
    ind3j[ii] = new int[lm]();
    MPI_Bcast(ind3j[ii], lm, MPI_INT, 0, MPI_COMM_WORLD);
  }
  sscs = new double[nsph]();
  MPI_Bcast(sscs, nsph, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  ecscm = new double[2]();
  scscm = new double[2]();
  scsc = new double[2]();
  ecsc = new double[2]();
  MPI_Bcast(ecscm, 2, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(scscm, 2, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(scsc, 2, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(ecsc, 2, MPI_DOUBLE, 0, MPI_COMM_WORLD);
}

void C1_AddOns::mpibcast(const mixMPI *mpidata) {
  MPI_Bcast(&nsph, 1, MPI_INT, 0, MPI_COMM_WORLD);
  MPI_Bcast(&lmpo, 1, MPI_INT, 0, MPI_COMM_WORLD);
  MPI_Bcast(&nlemt, 1, MPI_INT, 0, MPI_COMM_WORLD);
  MPI_Bcast(&litpo, 1, MPI_INT, 0, MPI_COMM_WORLD);
  MPI_Bcast(&lmtpo, 1, MPI_INT, 0, MPI_COMM_WORLD);
  MPI_Bcast(&litpos, 1, MPI_INT, 0, MPI_COMM_WORLD);
  MPI_Bcast(&lmtpos, 1, MPI_INT, 0, MPI_COMM_WORLD);
  MPI_Bcast(&nv3j, 1, MPI_INT, 0, MPI_COMM_WORLD);
  MPI_Bcast(&lm, 1, MPI_INT, 0, MPI_COMM_WORLD);
  int vhsize=(nsph * nsph - 1) * litpo;
  MPI_Bcast(vh, vhsize, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
  int vj0size=nsph * lmtpo;
  MPI_Bcast(vj0, vj0size, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
  MPI_Bcast(vj, 1, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
  int vyhjsize=(nsph * nsph - 1) * litpos;
  MPI_Bcast(vyhj, vyhjsize, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
  int vyj0size = nsph * lmtpos;
  MPI_Bcast(vyj0, vyj0size, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
  MPI_Bcast(am0v, nlemt*nlemt, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
  MPI_Bcast(vintm, 16, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
  MPI_Bcast(vintt, 16, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
  for (int fi = 0; fi < 2; fi++) {
    MPI_Bcast(fsac[fi], 2, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
    MPI_Bcast(sac[fi], 2, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
    MPI_Bcast(fsacm[fi], 2, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
  }
  MPI_Bcast(scscp, 2, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
  MPI_Bcast(ecscp, 2, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
  MPI_Bcast(scscpm, 2, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
  MPI_Bcast(ecscpm, 2, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
  MPI_Bcast(v3j0, nv3j, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  for (int ii = 0; ii < lmpo; ii++) {
    MPI_Bcast(ind3j[ii], lm, MPI_INT, 0, MPI_COMM_WORLD);
  }
  MPI_Bcast(sscs, nsph, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(ecscm, 2, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(scscm, 2, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(scsc, 2, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(ecsc, 2, MPI_DOUBLE, 0, MPI_COMM_WORLD);
}
#endif

C2::C2(GeometryConfiguration *gconf, ScattererConfiguration *sconf) {
  nsph = gconf->number_of_spheres;
  int npnt = gconf->npnt;
  int npntts = gconf->npntts;
  int max_n = (npnt > npntts) ? npnt : npntts;
  nhspo = 2 * max_n - 1;
  nl = sconf->configurations;
  if (nsph == 1 && nl == 1) nl = 5;
  ris = new dcomplex[nhspo]();
  dlri = new dcomplex[nhspo]();
  vkt = new dcomplex[nsph]();
  dc0 = new dcomplex[nl]();
  vsz = new double[nsph]();
}

C2::C2(const C2& rhs) {
  nsph = rhs.nsph;
  nhspo = rhs.nhspo;
  nl = rhs.nl;
  ris = new dcomplex[nhspo]();
  dlri = new dcomplex[nhspo]();
  for (int ind=0; ind<nhspo; ind++) {
    ris[ind] = rhs.ris[ind];
    dlri[ind] = rhs.dlri[ind];
  }
  vkt = new dcomplex[nsph]();
  vsz = new double[nsph]();
  for (int ind=0; ind<nsph; ind++) {
    vkt[ind] = rhs.vkt[ind];
    vsz[ind] = rhs.vsz[ind];
  }
  dc0 = new dcomplex[nl]();
  for (int ind=0; ind<nl; ind++) dc0[ind] = rhs.dc0[ind];
}

C2::~C2() {
  delete[] ris;
  delete[] dlri;
  delete[] vkt;
  delete[] dc0;
  delete[] vsz;
}

#ifdef MPI_VERSION
C2::C2(const mixMPI *mpidata) {
  MPI_Bcast(&nsph, 1, MPI_INT, 0, MPI_COMM_WORLD);
  MPI_Bcast(&nhspo, 1, MPI_INT, 0, MPI_COMM_WORLD);
  MPI_Bcast(&nl, 1, MPI_INT, 0, MPI_COMM_WORLD);
  ris = new dcomplex[nhspo]();
  dlri = new dcomplex[nhspo]();
  MPI_Bcast(ris, nhspo, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
  MPI_Bcast(dlri, nhspo, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
  vkt = new dcomplex[nsph]();
  vsz = new double[nsph]();
  MPI_Bcast(vkt, nsph, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
  MPI_Bcast(vsz, nsph, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  dc0 = new dcomplex[nl]();
  MPI_Bcast(dc0, nl, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
}

void C2::mpibcast(const mixMPI *mpidata) {
  MPI_Bcast(&nsph, 1, MPI_INT, 0, MPI_COMM_WORLD);
  MPI_Bcast(&nhspo, 1, MPI_INT, 0, MPI_COMM_WORLD);
  MPI_Bcast(&nl, 1, MPI_INT, 0, MPI_COMM_WORLD);
  MPI_Bcast(ris, nhspo, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
  MPI_Bcast(dlri, nhspo, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
  MPI_Bcast(vkt, nsph, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
  MPI_Bcast(vsz, nsph, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(dc0, nl, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
}
#endif

C3::C3() {
  tsas = new dcomplex*[2];
  tsas[0] = new dcomplex[2];
  tsas[1] = new dcomplex[2];
  tfsas = 0.0 + 0.0 * I;
  gcs = 0.0;
  scs = 0.0;
  ecs = 0.0;
  acs = 0.0;
}

C3::C3(const C3& rhs) {
  tsas = new dcomplex*[2];
  for (int ti=0; ti<2; ti++) {
    tsas[ti] = new dcomplex[2];
    for (int tj=0; tj<2; tj++) tsas[ti][tj] = rhs.tsas[ti][tj];
  }
  tfsas = rhs.tfsas;
  gcs = rhs.gcs;
  scs = rhs.scs;
  ecs = rhs.ecs;
  acs = rhs.acs;
}

C3::~C3() {
  delete[] tsas[1];
  delete[] tsas[0];
  delete[] tsas;
}

#ifdef MPI_VERSION
C3::C3(const mixMPI *mpidata) {
  tsas = new dcomplex*[2];
  for (int ti=0; ti<2; ti++) {
    tsas[ti] = new dcomplex[2];
    MPI_Bcast(tsas[ti], 2, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
  }
  MPI_Bcast(&tfsas, 1, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
  MPI_Bcast(&gcs, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(&scs, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(&ecs, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(&acs, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
}

void C3::mpibcast(const mixMPI *mpidata) {
  for (int ti=0; ti<2; ti++) {
    MPI_Bcast(tsas[ti], 2, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
  }
  MPI_Bcast(&tfsas, 1, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
  MPI_Bcast(&gcs, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(&scs, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(&ecs, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(&acs, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
}
#endif

C4::C4(GeometryConfiguration *gconf) {
  li = gconf->li;
  le = gconf->le;
  lm = (li > le) ? li : le;
  nv3j = (lm * (lm  + 1) * (2 * lm + 7)) / 6;
  nsph = gconf->number_of_spheres;
  // The following is needed to initialize C1_AddOns
  litpo = li + li + 1;
  litpos = litpo * litpo;
  lmtpo = li + le + 1;
  lmtpos = lmtpo * lmtpo;
  nlim = li * (li + 2);
  nlem = le * (le + 2);
  lmpo = lm + 1;
}

C4::C4(const C4& rhs) {
    li = rhs.li;
    le = rhs.le;
    lm = rhs.lm;
    nv3j = rhs.nv3j;
    nsph = rhs.nsph;
    // The following is needed to initialize C1_AddOns
    litpo = rhs.litpo;
    litpos = rhs.litpos;
    lmtpo = rhs.lmtpo;
    lmtpos = rhs.lmtpos;
    nlim = rhs.nlim;
    nlem = rhs.nlem;
    lmpo = rhs.lmpo;
}

C4::~C4() {
}

#ifdef MPI_VERSION
C4::C4(const mixMPI *mpidata) {
  MPI_Bcast(&li, 1, MPI_INT, 0, MPI_COMM_WORLD);
  MPI_Bcast(&le, 1, MPI_INT, 0, MPI_COMM_WORLD);
  MPI_Bcast(&lm, 1, MPI_INT, 0, MPI_COMM_WORLD);
  MPI_Bcast(&nv3j, 1, MPI_INT, 0, MPI_COMM_WORLD);
  MPI_Bcast(&nsph, 1, MPI_INT, 0, MPI_COMM_WORLD);
  // The following is needed to initialize C1_AddOns
  MPI_Bcast(&litpo, 1, MPI_INT, 0, MPI_COMM_WORLD);
  MPI_Bcast(&litpos, 1, MPI_INT, 0, MPI_COMM_WORLD);
  MPI_Bcast(&lmtpo, 1, MPI_INT, 0, MPI_COMM_WORLD);
  MPI_Bcast(&lmtpos, 1, MPI_INT, 0, MPI_COMM_WORLD);
  MPI_Bcast(&nlim, 1, MPI_INT, 0, MPI_COMM_WORLD);
  MPI_Bcast(&nlem, 1, MPI_INT, 0, MPI_COMM_WORLD);
  MPI_Bcast(&lmpo, 1, MPI_INT, 0, MPI_COMM_WORLD);
}

void C4::mpibcast(const mixMPI *mpidata) {
  MPI_Bcast(&li, 1, MPI_INT, 0, MPI_COMM_WORLD);
  MPI_Bcast(&le, 1, MPI_INT, 0, MPI_COMM_WORLD);
  MPI_Bcast(&lm, 1, MPI_INT, 0, MPI_COMM_WORLD);
  MPI_Bcast(&nv3j, 1, MPI_INT, 0, MPI_COMM_WORLD);
  MPI_Bcast(&nsph, 1, MPI_INT, 0, MPI_COMM_WORLD);
  // The following is needed to initialize C1_AddOns
  MPI_Bcast(&litpo, 1, MPI_INT, 0, MPI_COMM_WORLD);
  MPI_Bcast(&litpos, 1, MPI_INT, 0, MPI_COMM_WORLD);
  MPI_Bcast(&lmtpo, 1, MPI_INT, 0, MPI_COMM_WORLD);
  MPI_Bcast(&lmtpos, 1, MPI_INT, 0, MPI_COMM_WORLD);
  MPI_Bcast(&nlim, 1, MPI_INT, 0, MPI_COMM_WORLD);
  MPI_Bcast(&nlem, 1, MPI_INT, 0, MPI_COMM_WORLD);
  MPI_Bcast(&lmpo, 1, MPI_INT, 0, MPI_COMM_WORLD);
}
#endif

C6::C6(int _lmtpo) {
  lmtpo = _lmtpo;
  rac3j = new double[lmtpo]();
}

C6::C6(const C6& rhs) {
  lmtpo = rhs.lmtpo;
  rac3j = new double[lmtpo]();
  for (int ri = 0; ri < lmtpo; ri++) rac3j[ri] = rhs.rac3j[ri];
}

C6::~C6() {
  delete[] rac3j;
}

#ifdef MPI_VERSION
C6::C6(const mixMPI *mpidata) {
  MPI_Bcast(&lmtpo, 1, MPI_INT, 0, MPI_COMM_WORLD);
  rac3j = new double[lmtpo]();
  MPI_Bcast(rac3j, lmtpo, MPI_DOUBLE, 0, MPI_COMM_WORLD);
}

void C6::mpibcast(const mixMPI *mpidata) {
  MPI_Bcast(&lmtpo, 1, MPI_INT, 0, MPI_COMM_WORLD);
  MPI_Bcast(rac3j, lmtpo, MPI_DOUBLE, 0, MPI_COMM_WORLD);
}
#endif

C9::C9(int ndi, int _nlem, int ndit, int _nlemt) {
  gis_size_0 = ndi;
  sam_size_0 = ndit;
  nlem = _nlem;
  nlemt = _nlemt;
  gis = new dcomplex*[ndi];
  gls = new dcomplex*[ndi];
  // allocate these arrays to be contiguous, C-style
  gis[0] = new dcomplex[ndi*nlem]();
  gls[0]  = new dcomplex[ndi*nlem]();  
  for (int gi = 1; gi < ndi; gi++) {
    gis[gi] = gis[0] + gi*nlem;
    gls[gi] = gls[0] + gi*nlem;
  }
  sam = new dcomplex*[ndit];
  sam[0] = new dcomplex[ndit*nlemt]();
  for (int si = 1; si < ndit; si++) sam[si] = sam[0] + si*nlemt;
}

C9::C9(const C9& rhs) {
  gis_size_0 = rhs.gis_size_0;
  sam_size_0 = rhs.sam_size_0;
  nlem = rhs.nlem;
  nlemt = rhs.nlemt;
  gis = new dcomplex*[gis_size_0];
  gls = new dcomplex*[gis_size_0];
  // allocate these arrays to be contiguous, C-style
  gis[0] = new dcomplex[gis_size_0*nlem]();
  gls[0]  = new dcomplex[gis_size_0*nlem]();  
  memcpy(gis[0], rhs.gis[0], gis_size_0*nlem*sizeof(dcomplex));
  memcpy(gls[0], rhs.gls[0], gis_size_0*nlem*sizeof(dcomplex));
  for (int gi = 1; gi < gis_size_0; gi++) {
    gis[gi] = gis[0] + gi*nlem;
    gls[gi] = gls[0] + gi*nlem;
  }
  sam = new dcomplex*[sam_size_0];
  sam[0] = new dcomplex[sam_size_0*nlemt]();
  memcpy(sam[0], rhs.sam[0], sam_size_0*nlemt*sizeof(dcomplex));
  for (int si = 1; si < sam_size_0; si++) {
    sam[si] = sam[0]+si*nlemt;
  }
}

C9::~C9() {
  delete[] gis[0];
  delete[] gls[0];
  delete[] gis;
  delete[] gls;
  delete[] sam[0];
  delete[] sam;
}

#ifdef MPI_VERSION
C9::C9(const mixMPI *mpidata) {
  MPI_Bcast(&gis_size_0, 1, MPI_INT, 0, MPI_COMM_WORLD);
  MPI_Bcast(&sam_size_0, 1, MPI_INT, 0, MPI_COMM_WORLD);
  MPI_Bcast(&nlem, 1, MPI_INT, 0, MPI_COMM_WORLD);
  MPI_Bcast(&nlemt, 1, MPI_INT, 0, MPI_COMM_WORLD);
  gis = new dcomplex*[gis_size_0];
  gls = new dcomplex*[gis_size_0];
  gis[0] = new dcomplex[gis_size_0*nlem]();
  gls[0]  = new dcomplex[gis_size_0*nlem]();  
  MPI_Bcast(gis[0], gis_size_0*nlem, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
  MPI_Bcast(gls[0], gis_size_0*nlem, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
  for (int gi = 1; gi < gis_size_0; gi++) {
    gis[gi] = gis[0] + gi*nlem;
    gls[gi] = gls[0] + gi*nlem;
  }
  sam = new dcomplex*[sam_size_0];
  sam[0] = new dcomplex[sam_size_0*nlemt]();
  MPI_Bcast(sam[0], sam_size_0*nlemt, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
  for (int si = 1; si < sam_size_0; si++) {
    sam[si] = sam[0]+si*nlemt;
  }
}

void C9::mpibcast(const mixMPI *mpidata) {
  MPI_Bcast(&gis_size_0, 1, MPI_INT, 0, MPI_COMM_WORLD);
  MPI_Bcast(&sam_size_0, 1, MPI_INT, 0, MPI_COMM_WORLD);
  MPI_Bcast(&nlem, 1, MPI_INT, 0, MPI_COMM_WORLD);
  MPI_Bcast(&nlemt, 1, MPI_INT, 0, MPI_COMM_WORLD);
  MPI_Bcast(gis[0], gis_size_0*nlem, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
  MPI_Bcast(gls[0], gis_size_0*nlem, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
  MPI_Bcast(sam[0], sam_size_0*nlemt, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
}
#endif

mixMPI::mixMPI() {
  mpirunning = 0;
  rank = 0;
  nprocs = 1;
}

#ifdef MPI_VERSION
mixMPI::mixMPI(MPI_Comm comm){
  mpirunning = 1;
  int ierr;
  // we should add some meaningful error checking and management here
  ierr = MPI_Comm_rank(MPI_COMM_WORLD, &rank);
  ierr = MPI_Comm_size(MPI_COMM_WORLD, &nprocs);
}
#endif

mixMPI::mixMPI(const mixMPI& rhs) {
  mpirunning = rhs.mpirunning;
  rank = rhs.rank;
  nprocs = rhs.nprocs;
}

mixMPI::~mixMPI() {
}

ClusterIterationData::ClusterIterationData(GeometryConfiguration *gconf, ScattererConfiguration *sconf, const mixMPI *mpidata) {
  c1 = new C1(gconf, sconf);
  c2 = new C2(gconf, sconf);
  c3 = new C3();
  c4 = new C4(gconf);
  c1ao = new C1_AddOns(c4);
  c6 = new C6(c4->lmtpo);
  const int ndi = c4->nsph * c4->nlim;
  const np_int ndit = 2 * ndi;
  c9 = new C9(ndi, c4->nlem, 2 * ndi, 2 * c4->nlem);
  gaps = new double[c4->nsph]();
  tqev = new double[3]();
  tqsv = new double[3]();
  tqse = new double*[2];
  tqspe = new dcomplex*[2];
  tqss = new double*[2];
  tqsps = new dcomplex*[2];
  tqce = new double*[2];
  tqcpe = new dcomplex*[2];
  tqcs = new double*[2];
  tqcps = new dcomplex*[2];
  for (int ti = 0; ti < 2; ti++) {
    tqse[ti] = new double[c4->nsph]();
    tqspe[ti] = new dcomplex[c4->nsph]();
    tqss[ti] = new double[c4->nsph]();
    tqsps[ti] = new dcomplex[c4->nsph]();
    tqce[ti] = new double[3]();
    tqcpe[ti] = new dcomplex[3]();
    tqcs[ti] = new double[3]();
    tqcps[ti] = new dcomplex[3]();
  }
  gapv = new double[3]();
  gapp = new dcomplex*[3];
  gappm = new dcomplex*[3];
  gap = new double*[3];
  gapm = new double*[3];
  for (int gi = 0; gi < 3; gi++) {
    gapp[gi] = new dcomplex[2]();
    gappm[gi] = new dcomplex[2]();
    gap[gi] = new double[2]();
    gapm[gi] = new double[2]();
  }
  u = new double[3]();
  us = new double[3]();
  un = new double[3]();
  uns = new double[3]();
  up = new double[3]();
  ups = new double[3]();
  unmp = new double[3]();
  unsmp = new double[3]();
  upmp = new double[3]();
  upsmp = new double[3]();
  argi = new double[1]();
  args = new double[1]();
  duk = new double[3]();
  cextlr = new double*[4];
  cext = new double*[4];
  cmullr = new double*[4];;
  cmul = new double*[4];
  for (int ci = 0; ci < 4; ci++) {
    cextlr[ci] = new double[4]();
    cext[ci] = new double[4]();
    cmullr[ci] = new double[4]();
    cmul[ci] = new double[4]();
  }
  zpv = new double***[c4->lm];
  for (int zi = 0; zi < c4->lm; zi++) {
    zpv[zi] = new double**[3];
    for (int zj = 0; zj < 3; zj++) {
      zpv[zi][zj] = new double*[2];
      for (int zk = 0; zk < 2; zk++) {
	zpv[zi][zj][zk] = new double[2]();
      }
    }
  }
  am_vector = new dcomplex[ndit * ndit]();
  am = new dcomplex*[ndit];
  for (int ai = 0; ai < ndit; ai++) {
    am[ai] = (am_vector + ai * ndit);
  }
  
  arg = 0.0 + 0.0 * I;
  // These are suspect initializations
  scan = 0.0;
  cfmp = 0.0;
  sfmp = 0.0;
  cfsp = 0.0;
  sfsp = 0.0;
  qsfi = 0.0;
  // End of suspect initializations
  wn = sconf->wp / 3.0e8;
  xip = sconf->xip;
  sqsfi = 1.0;
  vk = 0.0;
  number_of_scales = sconf->number_of_scales;
  xiblock = (int) ceil(((double) (sconf->number_of_scales-1))/((double) mpidata->nprocs));
  lastxi = ((mpidata->rank+1) * xiblock)+1;
  firstxi = lastxi-xiblock+1;
  if (lastxi > sconf->number_of_scales) lastxi = sconf->number_of_scales;
}

ClusterIterationData::ClusterIterationData(const ClusterIterationData& rhs) {
  c1 = new C1(*(rhs.c1));
  c2 = new C2(*(rhs.c2));
  c3 = new C3(*(rhs.c3));
  c4 = new C4(*(rhs.c4));
  c1ao = new C1_AddOns(*(rhs.c1ao));
  c6 = new C6(*(rhs.c6));
  const int ndi = c4->nsph * c4->nlim;
  const np_int ndit = 2 * ndi;
  c9 = new C9(*(rhs.c9));
  gaps = new double[c4->nsph]();
  for (int gi = 0; gi < c4->nsph; gi++) gaps[gi] = rhs.gaps[gi];
  tqev = new double[3]();
  tqsv = new double[3]();
  for (int ti = 0; ti < 3; ti++) {
    tqev[ti] = rhs.tqev[ti];
    tqsv[ti] = rhs.tqsv[ti];
  }
  tqse = new double*[2];
  tqspe = new dcomplex*[2];
  tqss = new double*[2];
  tqsps = new dcomplex*[2];
  tqce = new double*[2];
  tqcpe = new dcomplex*[2];
  tqcs = new double*[2];
  tqcps = new dcomplex*[2];
  for (int ti = 0; ti < 2; ti++) {
    tqse[ti] = new double[c4->nsph]();
    tqspe[ti] = new dcomplex[c4->nsph]();
    tqss[ti] = new double[c4->nsph]();
    tqsps[ti] = new dcomplex[c4->nsph]();
    for (int tj = 0; tj < c4->nsph; tj++) {
      tqse[ti][tj] = rhs.tqse[ti][tj];
      tqspe[ti][tj] = rhs.tqspe[ti][tj];
      tqss[ti][tj] = rhs.tqss[ti][tj];
      tqsps[ti][tj] = rhs.tqsps[ti][tj];
    }
    tqce[ti] = new double[3]();
    tqcpe[ti] = new dcomplex[3]();
    tqcs[ti] = new double[3]();
    tqcps[ti] = new dcomplex[3]();
    for (int tj = 0; tj < 3; tj++) {
      tqce[ti][tj] = rhs.tqce[ti][tj];
      tqcpe[ti][tj] = rhs.tqcpe[ti][tj];
      tqcs[ti][tj] = rhs.tqcs[ti][tj];
      tqcps[ti][tj] = rhs.tqcps[ti][tj];
    }
  }
  gapv = new double[3]();
  gapp = new dcomplex*[3];
  gappm = new dcomplex*[3];
  gap = new double*[3];
  gapm = new double*[3];
  for (int gi = 0; gi < 3; gi++) {
    gapv[gi] = rhs.gapv[gi];
    gapp[gi] = new dcomplex[2]();
    gappm[gi] = new dcomplex[2]();
    gap[gi] = new double[2]();
    gapm[gi] = new double[2]();
    for (int gj = 0; gj < 2; gj++) {
      gapp[gi][gj] = rhs.gapp[gi][gj];
      gappm[gi][gj] = rhs.gappm[gi][gj];
      gap[gi][gj] = rhs.gap[gi][gj];
      gapm[gi][gj] = rhs.gapm[gi][gj];
    }
  }
  u = new double[3]();
  us = new double[3]();
  un = new double[3]();
  uns = new double[3]();
  up = new double[3]();
  ups = new double[3]();
  unmp = new double[3]();
  unsmp = new double[3]();
  upmp = new double[3]();
  upsmp = new double[3]();
  duk = new double[3]();
  for (int ui = 0; ui < 3; ui++) {
    u[ui] = rhs.u[ui];
    us[ui] = rhs.us[ui];
    un[ui] = rhs.un[ui];
    uns[ui] = rhs.uns[ui];
    up[ui] = rhs.up[ui];
    ups[ui] = rhs.ups[ui];
    unmp[ui] = rhs.unmp[ui];
    unsmp[ui] = rhs.unsmp[ui];
    upmp[ui] = rhs.upmp[ui];
    upsmp[ui] = rhs.upsmp[ui];
    duk[ui] = rhs.duk[ui];
  }
  argi = new double[1]();
  args = new double[1]();
  argi[0] = rhs.argi[0];
  args[0] = rhs.args[0];
  cextlr = new double*[4];
  cext = new double*[4];
  cmullr = new double*[4];;
  cmul = new double*[4];
  for (int ci = 0; ci < 4; ci++) {
    cextlr[ci] = new double[4]();
    cext[ci] = new double[4]();
    cmullr[ci] = new double[4]();
    cmul[ci] = new double[4]();
    for (int cj = 0; cj < 4; cj++) {
      cextlr[ci][cj] = rhs.cextlr[ci][cj];
      cext[ci][cj] = rhs.cext[ci][cj];
      cmullr[ci][cj] = rhs.cmullr[ci][cj];
      cmul[ci][cj] = rhs.cmul[ci][cj];
    }
  }
  zpv = new double***[c4->lm];
  for (int zi = 0; zi < c4->lm; zi++) {
    zpv[zi] = new double**[3];
    for (int zj = 0; zj < 3; zj++) {
      zpv[zi][zj] = new double*[2];
      for (int zk = 0; zk < 2; zk++) {
	zpv[zi][zj][zk] = new double[2]();
	zpv[zi][zj][zk][0] = rhs.zpv[zi][zj][zk][0];
	zpv[zi][zj][zk][1] = rhs.zpv[zi][zj][zk][1];
      }
    }
  }
  am_vector = new dcomplex[ndit * ndit]();
  for (np_int ai = 0; ai < ndit * ndit; ai++) am_vector[ai] = rhs.am_vector[ai];
  am = new dcomplex*[ndit];
  for (np_int ai = 0; ai < ndit; ai++) {
    am[ai] = (am_vector + ai * ndit);
  }
  
  arg = rhs.arg;
  // These are suspect initializations
  scan = rhs.scan;
  cfmp = rhs.cfmp;
  sfmp = rhs.sfmp;
  cfsp = rhs.cfsp;
  sfsp = rhs.sfsp;
  qsfi = rhs.qsfi;
  // End of suspect initializations
  wn = rhs.wn;
  xip = rhs.xip;
  sqsfi = rhs.sqsfi;
  vk = rhs.vk;
  firstxi = rhs.firstxi;
  lastxi = rhs.lastxi;
  xiblock = rhs.xiblock;
  number_of_scales = rhs.number_of_scales;
}

#ifdef MPI_VERSION
ClusterIterationData::ClusterIterationData(const mixMPI *mpidata) {
  c1 = new C1(mpidata);
  c2 = new C2(mpidata);
  c3 = new C3(mpidata);
  c4 = new C4(mpidata);
  c1ao = new C1_AddOns(mpidata);
  c6 = new C6(mpidata);
  const int ndi = c4->nsph * c4->nlim;
  const np_int ndit = 2 * ndi;
  c9 = new C9(mpidata);
  gaps = new double[c4->nsph]();
  MPI_Bcast(gaps, c4->nsph, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  tqev = new double[3]();
  tqsv = new double[3]();
  MPI_Bcast(tqev, 3, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(tqsv, 3, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  tqse = new double*[2];
  tqspe = new dcomplex*[2];
  tqss = new double*[2];
  tqsps = new dcomplex*[2];
  tqce = new double*[2];
  tqcpe = new dcomplex*[2];
  tqcs = new double*[2];
  tqcps = new dcomplex*[2];
  for (int ti = 0; ti < 2; ti++) {
    tqse[ti] = new double[c4->nsph]();
    tqspe[ti] = new dcomplex[c4->nsph]();
    tqss[ti] = new double[c4->nsph]();
    tqsps[ti] = new dcomplex[c4->nsph]();
    MPI_Bcast(tqse[ti], c4->nsph, MPI_DOUBLE, 0, MPI_COMM_WORLD);
    MPI_Bcast(tqspe[ti], c4->nsph, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
    MPI_Bcast(tqss[ti], c4->nsph, MPI_DOUBLE, 0, MPI_COMM_WORLD);
    MPI_Bcast(tqsps[ti], c4->nsph, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
    tqce[ti] = new double[3]();
    tqcpe[ti] = new dcomplex[3]();
    tqcs[ti] = new double[3]();
    tqcps[ti] = new dcomplex[3]();
    MPI_Bcast(tqce[ti], 3, MPI_DOUBLE, 0, MPI_COMM_WORLD);
    MPI_Bcast(tqcpe[ti], 3, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
    MPI_Bcast(tqcs[ti], 3, MPI_DOUBLE, 0, MPI_COMM_WORLD);
    MPI_Bcast(tqcps[ti], 3, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
  }
  gapv = new double[3]();
  gapp = new dcomplex*[3];
  gappm = new dcomplex*[3];
  gap = new double*[3];
  gapm = new double*[3];
  MPI_Bcast(gapv, 3, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  for (int gi = 0; gi < 3; gi++) {
    gapp[gi] = new dcomplex[2]();
    gappm[gi] = new dcomplex[2]();
    gap[gi] = new double[2]();
    gapm[gi] = new double[2]();
    MPI_Bcast(gapp[gi], 2, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
    MPI_Bcast(gappm[gi], 2, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
    MPI_Bcast(gap[gi], 2, MPI_DOUBLE, 0, MPI_COMM_WORLD);
    MPI_Bcast(gapm[gi], 2, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  }
  u = new double[3]();
  us = new double[3]();
  un = new double[3]();
  uns = new double[3]();
  up = new double[3]();
  ups = new double[3]();
  unmp = new double[3]();
  unsmp = new double[3]();
  upmp = new double[3]();
  upsmp = new double[3]();
  duk = new double[3]();
  MPI_Bcast(u, 3, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(us, 3, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(un, 3, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(uns, 3, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(up, 3, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(ups, 3, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(unmp, 3, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(unsmp, 3, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(upmp, 3, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(upsmp, 3, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(duk, 3, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  argi = new double[1]();
  args = new double[1]();
  MPI_Bcast(argi, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(args, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  cextlr = new double*[4];
  cext = new double*[4];
  cmullr = new double*[4];;
  cmul = new double*[4];
  for (int ci = 0; ci < 4; ci++) {
    cextlr[ci] = new double[4]();
    cext[ci] = new double[4]();
    cmullr[ci] = new double[4]();
    cmul[ci] = new double[4]();
    MPI_Bcast(cextlr[ci], 4, MPI_DOUBLE, 0, MPI_COMM_WORLD);
    MPI_Bcast(cext[ci], 4, MPI_DOUBLE, 0, MPI_COMM_WORLD);
    MPI_Bcast(cmullr[ci], 4, MPI_DOUBLE, 0, MPI_COMM_WORLD);
    MPI_Bcast(cmul[ci], 4, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  }
  zpv = new double***[c4->lm];
  for (int zi = 0; zi < c4->lm; zi++) {
    zpv[zi] = new double**[3];
    for (int zj = 0; zj < 3; zj++) {
      zpv[zi][zj] = new double*[2];
      for (int zk = 0; zk < 2; zk++) {
	zpv[zi][zj][zk] = new double[2]();
	MPI_Bcast(zpv[zi][zj][zk], 2, MPI_DOUBLE, 0, MPI_COMM_WORLD);
      }
    }
  }
  am_vector = new dcomplex[ndit * ndit]();
  am = new dcomplex*[ndit];
  for (np_int ai = 0; ai < ndit; ai++) {
    am[ai] = (am_vector + ai * ndit);
    MPI_Bcast(am[ai], ndit, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
  }
  MPI_Bcast(&arg, 1, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
  MPI_Bcast(&scan, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(&cfmp, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(&sfmp, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(&cfsp, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(&sfsp, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(&qsfi, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(&wn, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(&xip, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(&sqsfi, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(&vk, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(&xiblock, 1, MPI_INT, 0, MPI_COMM_WORLD);
  MPI_Bcast(&number_of_scales, 1, MPI_INT, 0, MPI_COMM_WORLD);
  lastxi = ((mpidata->rank+1) * xiblock)+1;
  firstxi = lastxi-xiblock+1;
  if (lastxi > number_of_scales) lastxi = number_of_scales;
}

void ClusterIterationData::mpibcast(const mixMPI *mpidata) {
  c1->mpibcast(mpidata);
  c2->mpibcast(mpidata);
  c3->mpibcast(mpidata);
  c4->mpibcast(mpidata);
  c1ao->mpibcast(mpidata);
  c6->mpibcast(mpidata);
  const int ndi = c4->nsph * c4->nlim;
  const np_int ndit = 2 * ndi;
  c9->mpibcast(mpidata);
  MPI_Bcast(gaps, c4->nsph, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(tqev, 3, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(tqsv, 3, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  for (int ti = 0; ti < 2; ti++) {
    MPI_Bcast(tqse[ti], c4->nsph, MPI_DOUBLE, 0, MPI_COMM_WORLD);
    MPI_Bcast(tqspe[ti], c4->nsph, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
    MPI_Bcast(tqss[ti], c4->nsph, MPI_DOUBLE, 0, MPI_COMM_WORLD);
    MPI_Bcast(tqsps[ti], c4->nsph, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
    MPI_Bcast(tqce[ti], 3, MPI_DOUBLE, 0, MPI_COMM_WORLD);
    MPI_Bcast(tqcpe[ti], 3, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
    MPI_Bcast(tqcs[ti], 3, MPI_DOUBLE, 0, MPI_COMM_WORLD);
    MPI_Bcast(tqcps[ti], 3, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
  }
  MPI_Bcast(gapv, 3, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  for (int gi = 0; gi < 3; gi++) {
    MPI_Bcast(gapp[gi], 2, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
    MPI_Bcast(gappm[gi], 2, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
    MPI_Bcast(gap[gi], 2, MPI_DOUBLE, 0, MPI_COMM_WORLD);
    MPI_Bcast(gapm[gi], 2, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  }
  MPI_Bcast(u, 3, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(us, 3, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(un, 3, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(uns, 3, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(up, 3, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(ups, 3, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(unmp, 3, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(unsmp, 3, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(upmp, 3, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(upsmp, 3, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(duk, 3, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(argi, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(args, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  for (int ci = 0; ci < 4; ci++) {
    MPI_Bcast(cextlr[ci], 4, MPI_DOUBLE, 0, MPI_COMM_WORLD);
    MPI_Bcast(cext[ci], 4, MPI_DOUBLE, 0, MPI_COMM_WORLD);
    MPI_Bcast(cmullr[ci], 4, MPI_DOUBLE, 0, MPI_COMM_WORLD);
    MPI_Bcast(cmul[ci], 4, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  }
  for (int zi = 0; zi < c4->lm; zi++) {
    for (int zj = 0; zj < 3; zj++) {
      for (int zk = 0; zk < 2; zk++) {
	MPI_Bcast(zpv[zi][zj][zk], 2, MPI_DOUBLE, 0, MPI_COMM_WORLD);
      }
    }
  }
  // since MPI expects an int argument for the number of elements to transfer in one go, transfer am one row at a time
  for (int ai = 0; ai < ndit; ai++) {
    MPI_Bcast(am[ai], ndit, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
  }
  MPI_Bcast(&arg, 1, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
  MPI_Bcast(&scan, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(&cfmp, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(&sfmp, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(&cfsp, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(&sfsp, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(&qsfi, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(&wn, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(&xip, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(&sqsfi, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(&vk, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(&xiblock, 1, MPI_INT, 0, MPI_COMM_WORLD);
  MPI_Bcast(&number_of_scales, 1, MPI_INT, 0, MPI_COMM_WORLD);
}
#endif

ClusterIterationData::~ClusterIterationData() {
  const int nsph = c4->nsph;
  delete[] am_vector;
  delete[] am;
  for (int zi = c4->lm - 1; zi > -1; zi--) {
    for (int zj = 2; zj > -1; zj--) {
      delete[] zpv[zi][zj][1];
      delete[] zpv[zi][zj][0];
      delete[] zpv[zi][zj];
    }
    delete[] zpv[zi];
  }
  delete[] zpv;
  delete c1;
  delete c1ao;
  delete c2;
  delete c3;
  delete c4;
  delete c6;
  delete c9;
  delete[] gaps;
  for (int ti = 1; ti > -1; ti--) {
    delete[] tqse[ti];
    delete[] tqss[ti];
    delete[] tqspe[ti];
    delete[] tqsps[ti];
    delete[] tqce[ti];
    delete[] tqcpe[ti];
    delete[] tqcs[ti];
    delete[] tqcps[ti];
  }
  delete[] tqse;
  delete[] tqss;
  delete[] tqspe;
  delete[] tqsps;
  delete[] tqce;
  delete[] tqcpe;
  delete[] tqcs;
  delete[] tqcps;
  delete[] tqev;
  delete[] tqsv;
  for (int gi = 2; gi > -1; gi--) {
    delete[] gapp[gi];
    delete[] gappm[gi];
    delete[] gap[gi];
    delete[] gapm[gi];
  }
  delete[] gapp;
  delete[] gappm;
  delete[] gap;
  delete[] gapm;
  delete[] gapv;
  delete[] u;
  delete[] us;
  delete[] un;
  delete[] uns;
  delete[] up;
  delete[] ups;
  delete[] unmp;
  delete[] unsmp;
  delete[] upmp;
  delete[] upsmp;
  delete[] argi;
  delete[] args;
  delete[] duk;
  for (int ci = 3; ci > -1; ci--) {
    delete[] cextlr[ci];
    delete[] cext[ci];
    delete[] cmullr[ci];
    delete[] cmul[ci];
  }
  delete[] cextlr;
  delete[] cext;
  delete[] cmullr;
  delete[] cmul;
}

ScatteringAngles::ScatteringAngles(GeometryConfiguration *gconf) {
  int isam = gconf->isam;
  _th = gconf->in_theta_start;
  _thstp = gconf->in_theta_step;
  _thlst = gconf->in_theta_end;
  _ths = gconf->sc_theta_start;
  _thsstp = gconf->sc_theta_step;
  _thslst = gconf->sc_theta_end;
  _ph = gconf->in_phi_start;
  _phstp = gconf->in_phi_step;
  _phlst = gconf->in_phi_end;
  _phs = gconf->sc_phi_start;
  _phsstp = gconf->sc_phi_step;
  _phslst = gconf->sc_phi_end;
  double small = 1.0e-3;
  _nth = 0;
  _nph = 0;
  if (_thstp != 0.0) _nth = int((_thlst - _th) / _thstp + small);
  _nth++;
  if (_phstp != 0.0) _nph = int((_phlst - _ph) / _phstp + small);
  _nph++;
  _nths = 0;
  _nphs = 0;
  _thsca = 0.0;
  if (isam > 1) {
    _nths = 1;
    _thsca = _ths - _th;
  } else { // ISAM <= 1
    if (_thsstp == 0.0) _nths = 0;
    else _nths = int ((_thslst - _ths) / _thsstp + small);
    _nths++;
  }
  if (isam >= 1) {
    _nphs = 1;
  } else {
    if (_phsstp == 0.0) _nphs = 0;
    else _nphs = int((_phslst - _phs) / _phsstp + small);
    _nphs++;
  }
  _nk = nth * nph;
  _nks = nths * nphs;
  _nkks = nk * nks;
}

ScatteringAngles::ScatteringAngles(const ScatteringAngles &rhs) {
  _th = rhs._th;
  _thstp = rhs._thstp;
  _thlst = rhs._thlst;
  _ths = rhs._ths;
  _thsstp = rhs._thsstp;
  _thslst = rhs._thslst;
  _ph = rhs._ph;
  _phstp = rhs._phstp;
  _phlst = rhs._phlst;
  _phs = rhs._phs;
  _phsstp = rhs._phsstp;
  _phslst = rhs._phslst;
  _nth = rhs._nth;
  _nph = rhs._nph;
  _nths = rhs._nths;
  _nphs = rhs._nphs;
  _thsca = rhs._thsca;
  _nk = rhs._nk;
  _nks = rhs._nks;
  _nkks = rhs._nkks;
}

#ifdef MPI_VERSION
ScatteringAngles::ScatteringAngles(const mixMPI *mpidata) {
  MPI_Bcast(&_th, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(&_thstp, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(&_thlst, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(&_ths, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(&_thsstp, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(&_thslst, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(&_ph, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(&_phstp, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(&_phlst, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(&_phs, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(&_phsstp, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(&_phslst, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(&_nth, 1, MPI_INT, 0, MPI_COMM_WORLD);
  MPI_Bcast(&_nph, 1, MPI_INT, 0, MPI_COMM_WORLD);
  MPI_Bcast(&_nths, 1, MPI_INT, 0, MPI_COMM_WORLD);
  MPI_Bcast(&_nphs, 1, MPI_INT, 0, MPI_COMM_WORLD);
  MPI_Bcast(&_thsca, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(&_nk, 1, MPI_INT, 0, MPI_COMM_WORLD);
  MPI_Bcast(&_nks, 1, MPI_INT, 0, MPI_COMM_WORLD);
  MPI_Bcast(&_nkks, 1, MPI_INT, 0, MPI_COMM_WORLD);
}

void ScatteringAngles::mpibcast(const mixMPI *mpidata) {
  MPI_Bcast(&_th, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(&_thstp, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(&_thlst, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(&_ths, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(&_thsstp, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(&_thslst, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(&_ph, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(&_phstp, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(&_phlst, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(&_phs, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(&_phsstp, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(&_phslst, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(&_nth, 1, MPI_INT, 0, MPI_COMM_WORLD);
  MPI_Bcast(&_nph, 1, MPI_INT, 0, MPI_COMM_WORLD);
  MPI_Bcast(&_nths, 1, MPI_INT, 0, MPI_COMM_WORLD);
  MPI_Bcast(&_nphs, 1, MPI_INT, 0, MPI_COMM_WORLD);
  MPI_Bcast(&_thsca, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(&_nk, 1, MPI_INT, 0, MPI_COMM_WORLD);
  MPI_Bcast(&_nks, 1, MPI_INT, 0, MPI_COMM_WORLD);
  MPI_Bcast(&_nkks, 1, MPI_INT, 0, MPI_COMM_WORLD);
}
#endif
