galaxy.cc

/**
 * @file galaxy.cc
 * @author Luca Maccione, Daniele Gaggero
 * @email luca.maccione@desy.de
 * @email daniele.gaggero@sissa.it
 * @brief In this file all the classes related to the model of the galaxy are implemented.
 */

#include "geometry.h"
#include "galaxy.h"
#include "grid.h"
#include "gas.h"
#include "input.h"
#include "nucleilist.h"
#include "sources.h"
#include "eloss.h"
#include "fitsio.h"
#include "errorcode.h"
#include "bfield.h"
#include "diffusion.h"

#include <fstream>
#include <iostream>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_spline.h>

#define DIFFTHRESHOLD 0.2

using namespace std;

TConvectionVelocity::TConvectionVelocity(TGrid* coord, TGeometry* geom, Input* in, TSource* SourceTerm) {
    
  nrn_sn = SourceTerm->GetSource(in->xobs,in->yobs,in->zobs);
    
  if (in->feedback >1) cout << "Called ConvectionVelocity"  << endl;
  dvdz = in->dvdz;
  conv_index_radial = in->conv_index_radial;
  set_profile_conv = in->set_profile_conv;
  conv_threshold = in->conv_threshold;
  vector<double> zgrid = coord->GetZ();
  dimz = zgrid.size();
  double rem_vel[dimz];
  double velocity=0;
  
  char buff[1000];
  sprintf(buff,"ASCII_spectra/%s/convection.dat",in->run_id.c_str());
  ofstream datafile;
  if(in->write_flag) datafile.open(buff);

  if (coord->GetType() == "2D") {
    vector<double> rgrid = coord->GetR();
    dimr = rgrid.size();
        
    for (int ir=0; ir<dimr; ir++) {
      double radius = rgrid[ir];

      for (int iz = 0; iz<dimz; iz++) {
	double zeta = zgrid[iz];
                
	if(fabs(zeta)<=in->z_k) velocity = ((((in->v0)-(in->vb))*pow(zeta,2.0)/pow(in->z_k,2.0))+(in->vb))* GetProfile(radius,0,zeta,SourceTerm);
	if(fabs(zeta)>in->z_k) velocity = ((in->v0) + dvdz*(fabs(zeta)-in->z_k)) * GetProfile(radius,0,zeta,SourceTerm);

	if((in->vb)>(in->v0)) cerr << "WARNING: vb > v0!" << endl;
         	
	//smoothen drop to avoid non-numerical values in fluxes
	if(set_profile_conv == Radial && conv_index_radial>0.){
	  if(ir==dimr-1) velocity=0.; //last bin  = zero
	  if(ir==dimr-3) rem_vel[iz]=velocity;
	  if(ir==dimr-2)velocity=0.5*rem_vel[iz]; //second to last bin is set to 0.5 * third to last bin
	}
      
	if(in->write_flag) datafile << radius << " " << zeta << " " << velocity/km/Myr*kpc << endl;

	vc.push_back(velocity);
      }
    }

    //MW130705: CN coefficients now here for 2D, too
    for (unsigned int i = 0; i < dimr; ++i) {
      for (unsigned int k = 0; k < dimz; ++k) {

	double vCk = 0.0; // vC(i)
	double vCk1 = 0.0; // vC(i+1)
	double vC1k = 0.0; // vC(i-1)

	if ( coord->GetZ().at(k) > 0 ) {
	  vCk1 = 0.0;
	  vCk  = vc[conv_index(i,k)]/coord->GetDeltaZ_down(k);
	  vC1k = vc[conv_index(i,k-1)]/coord->GetDeltaZ_down(k);
	}
	else if ( coord->GetZ().at(k) < 0) {
	  vCk  = vc[conv_index(i,k)]/coord->GetDeltaZ_up(k);
	  vC1k = 0.0;
	  vCk1 = vc[conv_index(i,k+1)]/coord->GetDeltaZ_up(k);
	}
	else {
	  vCk  = vc[conv_index(i,k)]*2/(coord->GetDeltaZ_up(k) + coord->GetDeltaZ_down(k));
	  vC1k = -0.5*vc[conv_index(i,k-1)]/coord->GetDeltaZ_down(k);
	  vCk1 = -0.5*vc[conv_index(i,k+1)]/coord->GetDeltaZ_up(k);
	}

	CNconv_alpha1_z.push_back( vC1k );
	CNconv_alpha2_z.push_back( vCk );
	CNconv_alpha3_z.push_back( vCk1 );

#ifdef DEBUGMODE
	cout<<"[MW-DEBUG-CONV]"<<i<<" "<<k<<" "<<conv_index(i,k)<<" | ";
	cout<<vc[conv_index(i,k)]<<" | "<<vC1k<<" "<<vCk<<" "<<vCk1<<" | ";
	cout<<coord->GetDeltaZ_up(k)<<" "<<coord->GetDeltaZ_down(k)<<" "<<coord->GetDeltaZ(k)<<endl;
#endif
      }
    }
  } // end 2D
  else {
    vector<double> xgrid = coord->GetX();
    dimx = xgrid.size();
    vector<double> ygrid = coord->GetY();
    dimy = ygrid.size();
        
    for (int ix=0; ix<dimx; ix++) {
      double x = xgrid[ix];
      for (int iy=0; iy<dimy; iy++) {
	double y = ygrid[iy];
	for (int iz = 0; iz<dimz; iz++) {
	  double z = zgrid[iz];
                    
	  if(fabs(z)<=in->z_k) velocity = ((((in->v0)-(in->vb))*pow(z,2.0)/pow(in->z_k,2.0))+(in->vb))* GetProfile(x,y,z,SourceTerm);
	  if(fabs(z)>in->z_k) velocity = ((in->v0) + dvdz*(fabs(z)-in->z_k)) * GetProfile(x,y,z,SourceTerm);
        
	  if((in->vb)>(in->v0)) cerr << "WARNING: vb > v0!" << endl;

	  velocity*= max( min( pow(geom->GetPattern(ix,iy,iz), in->SA_convec), in->SA_cut_convec), 1./in->SA_cut_convec );
	  velocity *= pow( in->LB_convec, coord->IsInLocalBubble(xgrid[ix],ygrid[iy],zgrid[iz]) );

	  vc.push_back(velocity);
	}
      }
    }

    if(set_profile_conv == Radial && conv_index_radial>0.){
      for (int ix=0; ix<dimx; ix++) for (int iy=0; iy<dimy; iy++) for (int iz = 0; iz<dimz; iz++){
	    //smoothen drop to avoid non-numerical values in fluxes, SK 06/13
	    if(ix==dimx-2 && iy==dimy-2)  vc[conv_index(ix,iy,iz)]=0.5*vc[conv_index(ix-1,iy,iz)];//second to last bin is set to 0.5 * third to last bin
	    if(iy==dimy-2 && ix==dimx-2)  vc[conv_index(ix,iy,iz)]=0.5*vc[conv_index(ix,iy-1,iz)];//second to last bin is set to 0.5 * third to last bin
	    if(ix==1 && iy==1)  vc[conv_index(ix,iy,iz)]=0.5*vc[conv_index(ix+1,iy,iz)];//second to last bin is set to 0.5 * third to last bin
	    if(iy==1 && ix==1)  vc[conv_index(ix,iy,iz)]=0.5*vc[conv_index(ix,iy+1,iz)];//second to last bin is set to 0.5 * third to last bin
	    if(ix==dimx-1 && iy==dimy-1) vc[conv_index(ix,iy,iz)]=0.; //border bins==0
	    if(ix==0 && iy==0) vc[conv_index(ix,iy,iz)]=0.; //border bins==0
	    if(ix==dimx-1 && iy==0) vc[conv_index(ix,iy,iz)]=0.; //border bins==0
	    if(ix==0 && iy==dimy-1) vc[conv_index(ix,iy,iz)]=0.; //border bins==0
	  }
    }

    //MW130624: CN coefficients are now here, the Evolutor just calls these vectors
    for (unsigned int i = 0; i < dimx; ++i) {
      for (unsigned int j = 0; j < dimy; ++j) {
	for (unsigned int k = 0; k < dimz; ++k) {

	  double vCk = 0.0; // vC(i)
	  double vCk1 = 0.0; // vC(i+1)
	  double vC1k = 0.0; // vC(i-1)

	  if ( coord->GetZ().at(k) > 0 ) {
	    vCk1 = 0.0;
	    vCk  = vc[conv_index(i,j,k)]/coord->GetDeltaZ_down(k);
	    vC1k = vc[conv_index(i,j,k-1)]/coord->GetDeltaZ_down(k);
	  }
	  else if ( coord->GetZ().at(k) < 0) {
	    vCk  = vc[conv_index(i,j,k)]/coord->GetDeltaZ_up(k);
	    vC1k = 0.0;
	    vCk1 = vc[conv_index(i,j,k+1)]/coord->GetDeltaZ_up(k);
	  }
	  else {
	    vCk  = vc[conv_index(i,j,k)]*2/(coord->GetDeltaZ_up(k) + coord->GetDeltaZ_down(k));
	    vC1k = -0.5*vc[conv_index(i,j,k-1)]/coord->GetDeltaZ_down(k);
	    vCk1 = -0.5*vc[conv_index(i,j,k+1)]/coord->GetDeltaZ_up(k);
	  }

	  CNconv_alpha1_z.push_back( vC1k );
	  CNconv_alpha2_z.push_back( vCk );
	  CNconv_alpha3_z.push_back( vCk1 );

	}
      }            
    }
  }

  if(in->write_flag) datafile.close();

}

double TConvectionVelocity::GetProfile(double x, double y, double zeta, TSource* SourceTerm) {
	
  double radial = 1;
  double result = 1.;

  switch(set_profile_conv) {
  case Constant :
    return 1.0;
    break;
    
  case Radial : //MW130621: Qtau is doing the same as Radial from Convection_new, just not evaluating it again, but instead taking the value of the SourceTerm.

    //MW130711: z-dependency should not be accounted for!
    zeta = 0;
    
    //MW130711: for comparison with DRAGON-KIT, don't normalize.
    radial = SourceTerm->GetSource(x,y,zeta);
    if (radial < conv_threshold) radial=conv_threshold;
    return pow(radial, conv_index_radial);
    
  case Qtau :

    radial = SourceTerm->GetSource(x,y,zeta);
    radial /= nrn_sn;
    
    if (radial < conv_threshold) radial=conv_threshold;
    
    return pow(radial, conv_index_radial);
    break;

  default :
    return -1;
  }
  
}

TReaccelerationCoefficient::TReaccelerationCoefficient(vector<double> pp, TDiffusionCoefficient* dperp, TGeometry* geom, Input* in) {
  
  double a;
  double Dpp_constant[pp.size()];
  
  for (unsigned int i = 0; i < pp.size(); ++i){
    if(pp[i] < in->rho_b)
      a = (in->DiffT == Anisotropic) ? in->DeltaPar : dperp->GetDelta(); //MW130711: integrate Anisotropic Diffusion
    else
      a = (in->DiffT == Anisotropic) ? in->DeltaPar : dperp->GetDelta_h();
    
    Dpp_constant[i]= 1.0/(a*(4.-a)*(4.-a*a));       // Ptuskin-2003
    if(in->diff_reacc == 1) Dpp_constant[i] *= 4.0/3.0; // Seo & Ptuskin
  }
  
  vector<double> DiffSpectrum = dperp->GetSpectrum();
  for (unsigned int i = 0; i < pp.size(); ++i) sp.push_back(Dpp_constant[i]*in->vAlfven*in->vAlfven*pp[i]*pp[i]/DiffSpectrum[i]);
  
#ifdef DEBUGMODE
  for (unsigned int i = 0; i < pp.size(); ++i){
    cout<<"[MW-DEBUG REACC G] "<<" "<<i<<" "<<Dpp_constant[i]<<" "<<in->vAlfven<<" "<<pp[i]<<" ";
    cout<<DiffSpectrum[i]<<" | " <<in->DiffT<<" "<<in->DeltaPar<<" "<<" "<<dperp->GetDelta()<<" "<<dperp->GetDelta_h()<<" "<<endl;
  }
#endif

  dimr = dperp->GetDimR();
  dimx = dperp->GetDimX();
  dimy = dperp->GetDimY();
  dimz = dperp->GetDimZ();

  vector<double> DiffProfile = (in->DiffT == Anisotropic) ? dperp->GetDPar() : dperp->GetDiffusionCoefficient();
    
  if(dperp->GetCoord()->GetType() == "3D")
    {
      int index = 0;
      for (vector<double>::iterator i = DiffProfile.begin(); i != DiffProfile.end(); ++i)
        {
	  int ix = dperp->GetCoord()->GetXFromIndexD_3D(index);
	  int iy = dperp->GetCoord()->GetYFromIndexD_3D(index);
	  int iz = dperp->GetCoord()->GetZFromIndexD_3D(index);

	  double xx = dperp->GetCoord()->GetX()[ix];
	  double yy = dperp->GetCoord()->GetY()[iy];
	  double zz = dperp->GetCoord()->GetZ()[iz];

	  double reacc_spatial = 1.0/(*i);

	  double spiral_factor_dperp = max( min( pow(geom->GetPattern(ix,iy,iz), in->SA_diff), in->SA_cut_diff), 1./in->SA_cut_diff );
	  double spiral_factor_dpp = max( min( spiral_factor_dperp * pow(geom->GetPattern(ix,iy,iz), 2*in->SA_vA), in->SA_cut_vA), 1./in->SA_cut_vA );

	  reacc_spatial *= spiral_factor_dpp; //mw 130422
	  if (dperp->GetCoord()->IsInLocalBubble(xx,yy,zz)) reacc_spatial *= pow(in->LB_vA, 2*dperp->GetCoord()->IsInLocalBubble(xx,yy,zz)) * pow(in->LB_diff, dperp->GetCoord()->IsInLocalBubble(xx,yy,zz));
	  dpp.push_back(reacc_spatial);

	  index++;
        }
    }
  else
    {
      for (vector<double>::iterator i = DiffProfile.begin(); i != DiffProfile.end(); ++i)
        {
	  double reacc_spatial = 1.0/(*i);
	  dpp.push_back(reacc_spatial);
        }
    }
}

//************************************************************
//************** THE GALAXY CONSTRUCTOR **********************
//************************************************************

Galaxy::Galaxy(Input* inputStructure_, TNucleiList* nucleiList_) {
  
  if (inputStructure_ == NULL) {
    cerr<<"No Input specified!"<<endl;
    return ;
  }
  
  inputStructure = inputStructure_;
  const int feedback = inputStructure_->feedback;
  
  if (feedback >1) cout<<"Welcome to the Galaxy constructor!"<<endl;
  
  ifstream infile(inputStructure->sourcedata.c_str(),ios::in);
  
  if ( !infile.is_open() ){
    cerr<<"... file "<<inputStructure->sourcedata<<" does not exist. Using config_files/template.source.param!"<<endl;
    infile.open("config_files/template.source.param",ios::in);
    if ( !infile.is_open() ){
      cerr<<"WARNING! config_files/template.source.param does not exist, either! Exiting."<<endl;
      exit(NOSOURCEDATA);
    }
  }
    
  int particle_ID;
  map<int,double> abundances_map;
  map<int, vector<double> > inj_indexes;
  map<int, vector<double> > break_positions;

  //the code reads the .source.param. First column: nucleus ID; second column: abundance; 
  //other columns: inj_slope - break rigidity - inj slope - break rigidity - (...) - highest energy inj_slope; 
  //the number of breaks is arbitrary
  //reads the .source.param to a table	
  //DG29.09.2013
  
  typedef vector<double> Row;
  vector<Row> table;	 
  while (infile) {
    string line; getline(infile, line);
    istringstream temp_string(line);
    Row row;
    while (temp_string) {
      double data;
      temp_string >> data;
      row.push_back(data);
    }
    if (!infile.eof())
      table.push_back(row);
  }
  
  //fills abundances_map inj_indexes and break_positions for each nucleus i using the table
  if (feedback >0) cout<<"... reading source.param table with "<<table.size()<<" rows!"<<endl;
  
  for (unsigned int i=0; i<table.size(); i++) {
    Row row; row = table[i];
    if (feedback >1) cout<<"... reading line in .source.param of size: "<<row.size()<<endl;
    int nid = 0;	
    for (unsigned j=0; j<row.size()-1; j++) {
      if (feedback>1) cout<<table[i][j]<<", ";
      if (j==0) {
	nid = table[i][0]; //cout << nid << endl;
      }
      if (j==1) abundances_map[nid] = table[i][1];
      if (j>0 && j%2==0)
	inj_indexes[nid].push_back(table[i][j]);	
      if (j>1 && j%2!=0)
	break_positions[nid].push_back(table[i][j]);	
    }
    if (feedback>1) cout<<endl;
  }    
  
  vector<int> list = nucleiList_->GetList();
  
  // Remove the nucleus not included in the .source.param
  /*for (vector<int>::iterator it_current_nucleus = list.begin(); it_current_nucleus != list.end(); ++it_current_nucleus) {
    map<int,double>::iterator it_current_nucleus_abundance = abundances_map.find(*it_current_nucleus);
    if ( it_current_nucleus_abundance == abundances_map.end() ){
    if (feedback>0) cout<<"Remove nucleus "<<*it_current_nucleus<<" from the list."<<endl;
    nucleiList_->DeleteNucleusFromList( *it_current_nucleus );
    }
    }*/
  
  list.clear();
  
  list = nucleiList_->GetList();
  
  if (feedback>1) cout<<"Size nuclei list: "<<list.size()<<endl;
  
  if (feedback>1) cout<<"Setting abundances, inj slopes and break positions for each nucleus in the list"<<endl;
  
  for (vector<int>::iterator it_current_nucleus = list.begin(); it_current_nucleus != list.end(); ++it_current_nucleus) {
    
    // loop over NucleiList from nucleilist.cc
    
    map<int,double>::iterator it_current_nucleus_abundance = abundances_map.find(*it_current_nucleus);
    
    if ( it_current_nucleus_abundance != abundances_map.end() ){
      
      _fSourceAbundances[*it_current_nucleus] = (*it_current_nucleus_abundance).second;
      if (feedback>1){
	cout<<"Nucleus id -> "<<(*it_current_nucleus_abundance).first<<". Abundance found in .source.param -> ";
	cout<<_fSourceAbundances[*it_current_nucleus]<<endl;
      }
      
      if ( inputStructure->UseInjectionIndexAllNuclei == false ){	
	if (feedback>1) cout<<"Slopes are NOT specified in the XML and are taken from .source.param file!"<<endl;
	_fInjSpectrum_rho[*it_current_nucleus]   = break_positions[*it_current_nucleus];	
	_fInjSpectrum_alpha[*it_current_nucleus] = inj_indexes[*it_current_nucleus];
      }
      else {
	
	if (feedback>1){
	  cout<<"Break positions and slopes are taken from xml file!"<<endl;
	  cout<<"Number of slopes: "<<inputStructure->inp_inj_indexes.size()<<endl;
	  cout<<"Number of breaks: "<<inputStructure->inp_break_positions.size()<<endl;
	}

	for (int j=0; j<inputStructure->inp_break_positions.size(); j++)
	  _fInjSpectrum_rho[*it_current_nucleus] = inputStructure->inp_break_positions;
	for (int j=0; j<inputStructure->inp_inj_indexes.size(); j++)
	  _fInjSpectrum_alpha[*it_current_nucleus] = inputStructure->inp_inj_indexes;
      }	
    }
    else {
      
      if (feedback>1) cout<<"Nucleus id -> "<<(*it_current_nucleus_abundance).first<<". Abundance NOT found in .source.param!"<<endl;
      
      _fSourceAbundances[*it_current_nucleus] = 0.0;
      _fInjSpectrum_rho[*it_current_nucleus].push_back(1.);
      _fInjSpectrum_alpha[*it_current_nucleus].push_back(0.);
      _fInjSpectrum_alpha[*it_current_nucleus].push_back(0.);
    }
  }
  
  _fSourceAbundances[-1000] = 1.0;
  _fSourceAbundances[1000]  = 0.0;
  _fSourceAbundances[-999]  = 0.0;
  _fSourceAbundances[-1998] = 0.0;
  
  TESTMODE = inputStructure->TESTMODE;
  
  if (inputStructure->MOVING == false) MOVING = false;
  else MOVING = true;		
  
  if (MOVING)  {
    source_x0 = inputStructure->source_x0;
    source_y0 = inputStructure->source_y0;
    source_z0 = inputStructure->source_z0;
    source_vx = inputStructure->source_vx;
    source_vy = inputStructure->source_vy;
    source_vz = inputStructure->source_vz;
  } 
  
  if (inputStructure->MOVING_CLUMP == false) MOVING_CLUMP = false;
  else MOVING_CLUMP = true;		
  
  if (MOVING_CLUMP) { 
    clump_x0 = inputStructure->clump_x0;
    clump_y0 = inputStructure->clump_y0;
    clump_z0 = inputStructure->clump_z0;
    clump_vx = inputStructure->clump_vx;
    clump_vy = inputStructure->clump_vy;
    clump_vz = inputStructure->clump_vz;
    clump_deltat = inputStructure->clump_deltat;
  } 

#ifdef DEBUGMODE
  for (vector<int>::iterator it = list.begin(); it != list.end(); ++it){
    cout<<"injection "<<*it<<" "<<_fSourceAbundances[*it]<<" "<<_fInjSpectrum_rho_0[*it]<<" "<<_fInjSpectrum_rho_1[*it]<<" ";
    cout<<_fInjSpectrum_rho_2[*it]<<" "<<_fInjSpectrum_alpha_0[*it]<<" "<<_fInjSpectrum_alpha_1[*it]<<" "<<_fInjSpectrum_alpha_2[*it]<<" "<<_fInjSpectrum_alpha_3[*it]<<endl;
  }
#endif
    
  if (feedback>1) cout<<"Preparing the grid."<<endl; 
  
  if (inputStructure->gridtype == "2D") _fCoordinates =  new TGrid2D(inputStructure_);
  else _fCoordinates = new TGrid3D(inputStructure_);
  if (feedback>0) cout<<"... grid done!"<<endl;
  
  if (feedback>1) cout<<"Preparing the geometry."<<endl;
  if (inputStructure->SA_type == "None") _fGeometry = new TUniformGeometry(_fCoordinates, inputStructure_);
  else _fGeometry = new TSpiralGeometry(_fCoordinates, inputStructure_);
  if (feedback>0) cout<<"... geometry done!"<<endl;
  
  if (feedback>1) cout<<"Preparing the gas... "<<endl;
  _fGas.push_back(new TH2Gas(_fCoordinates, inputStructure_, _fGeometry));
  _fGas.push_back(new THIGas(_fCoordinates, inputStructure_, _fGeometry));
  _fGas.push_back(new THIIGas(_fCoordinates, inputStructure_, _fGeometry));
  
  _fTotalGas = new TGas(_fCoordinates, inputStructure_);
  if (feedback>1) cout<<"[MW] Sum of TotalGas vector before filling(should be zero): "<<_fTotalGas->GetTotalContent()<<endl;
  //MW 130429: Construct TotalGas from other components
  *_fTotalGas += *_fGas[0];
  *_fTotalGas += *_fGas[1];
  *_fTotalGas += *_fGas[2];
  if (feedback>1) cout<<"[MW] Sum of TotalGas vector after filling (should be "<<_fGas[0]->GetTotalContent()+_fGas[1]->GetTotalContent()+_fGas[2]->GetTotalContent()<<"): "<<_fTotalGas->GetTotalContent()<<endl;
  
  // SETTING GAS ABUNDANCES
  
  _fGasAbundances[1001] = 1.;
  _fGasAbundances[2004] = 0.11;
  _fGasAbundances[6012] = 0.05;

  if (feedback>0) cout<<"... gas done!"<<endl;
    
  _fDMSource = new TDMSource(_fCoordinates, inputStructure_);
  if (feedback>1) cout<<"Creating astrophysical source... "<<endl;
  _fSource   = new TAstrophysicalSource(_fCoordinates, inputStructure_, _fGeometry, inputStructure->SNR_model);
  _fSourceExtra = (inputStructure->prop_extracomp) ? new TAstrophysicalSource(_fCoordinates, inputStructure_, _fGeometry, inputStructure->SNR_model_Extra) : NULL; //CAREFUL! hard coded model = model_extra
  if (feedback>0) cout<<"... all astrophysical sources done!"<<endl;
  
  switch (inputStructure->BM) {
  case Pshirkov:
    _fB = new TPshirkovField(inputStructure->B0disk, 5., 1., 10., inputStructure->B0halo, 8., 1.3, inputStructure->B0turb, 8.5, inputStructure->zt, inputStructure->robs, _fCoordinates, _fGeometry);
    break;
  case Farrar:
    _fB = new TFarrarField(inputStructure->betaFarrar, _fCoordinates, _fGeometry);
    break;
  case Uniform:
    _fB = new TUniformField(inputStructure->B0turb, _fCoordinates, _fGeometry);
    break;
  case Simple:
    _fB = new TSimpleField(inputStructure->b0, inputStructure->robs, _fCoordinates, _fGeometry);
    break;
  case ToyModel:
    if (feedback >1) cout << "ToyModel mag field was specified!" <<endl;
    _fB = new ToyModelField(inputStructure->bx, inputStructure->by, inputStructure->bz, inputStructure->bturb, _fCoordinates, _fGeometry); 
    break;
  default :
    _fB = NULL;
  }
    
  _fDperp = NULL;
  _fDpp = NULL;
  _fDperpEl = NULL;
  _fDppEl = NULL;
    
  if (inputStructure->DiffT != Anisotropic) {
    
    if (inputStructure->gridtype == "2D") _fDperp = new TDiffusionCoefficient2D(_fCoordinates, inputStructure_, _fSource, _fB);
    else _fDperp = new TDiffusionCoefficient3D(_fCoordinates, inputStructure_, _fSource, _fB, _fGeometry, 0, 0);
    
    if (feedback >0) cout<<"... diffusion coefficient done!"<<endl;
    
    _fDpp = (inputStructure->REACC) ? new TReaccelerationCoefficient(_fCoordinates->GetMomentum(), _fDperp, _fGeometry, inputStructure_) : NULL;
  }
    
  _fVC = (inputStructure->CONVECTION) ? new TConvectionVelocity(_fCoordinates, _fGeometry, inputStructure_, _fSource) : NULL;
    
  if (inputStructure->prop_lep || inputStructure->prop_extracomp || inputStructure->prop_DMel) {
    
    if (inputStructure->DiffT != Anisotropic) {
            
      if (inputStructure->gridtype == "2D") _fDperpEl = new TDiffusionCoefficient2D(_fCoordinates, inputStructure_, _fSource, _fB, 1);
      else _fDperpEl = new TDiffusionCoefficient3D(_fCoordinates, inputStructure_, _fSource, _fB, _fGeometry, 0, 0, 1);
            
      _fDppEl = (inputStructure->REACC) ? new TReaccelerationCoefficient(_fCoordinates->GetMomentumEl(), _fDperpEl, _fGeometry, inputStructure_) : NULL;
    }
    _fISRF = new TISRF(_fCoordinates, ISRFfile, _fGeometry, inputStructure_);
  }
  else {
    _fDperpEl = NULL;
    _fDppEl = NULL;
    _fISRF = NULL;
  }
}

void Galaxy::Delete() {
  if (_fTotalGas) delete _fTotalGas;
  for (vector<TGas*>::iterator i = _fGas.begin(); i != _fGas.end(); ++i) {
    if (*i) delete *i;
  }
  _fGas.clear();
  if (_fISRF) delete _fISRF;
  if (_fB) delete _fB;
  if (_fGeometry) delete _fGeometry;
}

Galaxy::~Galaxy() {
  if (_fCoordinates) delete _fCoordinates;
  if (_fSource) delete _fSource;
  if (_fSourceExtra) delete _fSourceExtra;
  if (_fDMSource) delete _fDMSource;
  if (_fDperp) delete _fDperp;
  if (_fDpp) delete _fDpp;
  if (_fDperpEl) delete _fDperpEl;
  if (_fDppEl) delete _fDppEl;
  if (_fVC) delete _fVC;
  if (_fISRF) delete _fISRF;
  if (_fB) delete _fB;
  if (_fTotalGas) delete _fTotalGas;
  for (vector<TGas*>::iterator i = _fGas.begin(); i != _fGas.end(); ++i) {
    if (*i) delete *i;
  }
  _fGas.clear();

  _fSourceAbundances.clear();
  _fInjSpectrum_rho.clear();
  _fInjSpectrum_alpha.clear();
}