source/*cc

source/ArrayModel.cc

+#include "Model.h"
+
+ArrayModel::ArrayModel(const Exposure &exposure, const CountsMap &counts, const Psf &psf, const Sources &sources, const Parameters &pars, const Skymap<char> &filter, unsigned int configure) :
+	BaseModel(counts,exposure,psf,sources,pars, filter, configure), fmodels(0)
+{
+}
+
+void ArrayModel::setModels( const std::vector<BaseModel*> & models) {
+	//Check for null pointers in the array
+	for (int i = 0; i < models.size(); ++i){
+		if (models[i] == 0) {
+			throw(ModelError("Null pointer in array model, aborting"));
+		}
+	}
+	//Clear up the variables
+	fvariables.clear();
+	
+	//Set the model array to the new models
+	fmodels = models;
+
+	//Join the variables
+	for (int i = 0; i < fmodels.size(); ++i){
+		fvariables.add(fmodels[i]->getVariables());
+	}
+}
+
+void ArrayModel::getMap(const Variables & vars, Skymap<double> &map) {
+	for (int i = 0; i < fmodels.size(); ++i){
+		fmodels[i]->getMap(vars, map);
+	}
+}
+
+BaseModel::gradMap ArrayModel::getComponents(const Variables & vars, const std::string &prefix){
+	gradMap output;
+	for (int i = 0; i < fmodels.size(); ++i){
+		gradMap tmp(fmodels[i]->getComponents(vars, prefix));
+		for (gradMap::iterator it = tmp.begin(); it != tmp.end(); ++it){
+		   //Check if something exists, if it does, add the new map
+		   //rather than replace it
+		   gradMap::iterator it2 = output.find((*it).first);
+		   if (it2 == output.end()) {
+			output[(*it).first] = (*it).second;
+		   } else {
+			output[(*it).first] += (*it).second;
+		   }
+		}
+	}
+	return output;
+}
+
+void ArrayModel::getGrads(const Variables & vars, const std::string & varName, Skymap<double> &map) {
+	for (int i = 0; i < fmodels.size(); ++i){
+		fmodels[i]->getGrads(vars,varName, map);
+	}
+}
+
+BaseModel::gradMap ArrayModel::getSecDer(const Variables & vars, const std::string & varName1, const std::string & varName2) {
+	gradMap output;
+	//Loop over the models and add to output
+	for (int i = 0; i < fmodels.size(); ++i){
+		gradMap tmp(fmodels[i]->getSecDer(vars,varName1,varName2));
+		for (gradMap::iterator it = tmp.begin(); it != tmp.end(); ++it){
+		   //Check if something exists, if it does, add the new map
+		   //rather than replace it
+		   gradMap::iterator it2 = output.find((*it).first);
+		   if (it2 == output.end()) {
+			output[(*it).first] = (*it).second;
+		   } else {
+			output[(*it).first] += (*it).second;
+		   }
+		}
+	}
+	return output;
+}
+
+void ArrayModel::modifySources(const Variables & vars, Sources &sources) const {
+	for (int i = 0; i < fmodels.size(); ++i) {
+		fmodels[i]->modifySources(vars,sources);
+	}
+}

source/BaseModel.cc

+#include "Goodness.h"
+#include "Exposure.h"
+#include "Counts.h"
+#include "Psf.h"
+#include "Sources.h"
+#include "Parameters.h"
+#include "Skymap.h"
+#include <string>
+#include <cmath>
+#include <iostream>
+#include <algorithm>
+#include <iterator>
+#include <valarray>
+#include "alm.h"
+#ifdef HAVE_ALM_MAP_TOOLS_H
+   #include "alm_map_tools.h"
+#endif
+#ifdef HAVE_ALM_HEALPIX_TOOLS_H
+   #include "alm_healpix_tools.h"
+#endif
+#include "xcomplex.h"
+#include "arr.h"
+#include "time.h"
+#include "Utils.h"
+#include <ostream>
+
+BaseModel::BaseModel(const CountsMap & counts, const Exposure &exposure, const Psf &psf, const Sources &sources, const Parameters &pars, const Skymap<char> &filter, unsigned int configure) : 
+   fcounts(counts), 
+   fexposure(exposure), 
+   fpsf(psf), 
+   fsources(sources), 
+   fparameters(pars), 
+   ffilter(filter),
+   fconfigure(configure),
+   EXPOSURE_CORRECT(1),
+   CONVOLVE(2) {
+	//Sanity check, to see that the energy range of the exposure and psf is
+	//greater than that of the counts.
+	if (counts.getEMin().min() < exposure.getEnergies().min()){
+		std::cout<<"Warning: "<<std::endl;
+		std::cout<<"Lower energy bound of counts is lower than the lowest energy of exposure, continuing with extrapolation"<<std::endl;
+		std::cout<<counts.getEMin().min()<<" < "<<exposure.getEnergies().min()<<std::endl;
+	}
+	if (counts.getEMin().min() < psf.getEnergies().front()){
+		std::cout<<"Warning: "<<std::endl;
+		std::cout<<"Lower energy bound of counts is lower than the lowest energy of psf, continuing with extrapolation"<<std::endl;
+		std::cout<<counts.getEMin().min()<<" < "<<psf.getEnergies().front()<<std::endl;
+	}
+	if (counts.getEMax().max() > exposure.getEnergies().max()){
+		std::cout<<"Warning: "<<std::endl;
+		std::cout<<"Higher energy bound of counts is higher than the highest energy of exposure, continuing with extrapolation"<<std::endl;
+		std::cout<<counts.getEMax().max()<<" > "<<exposure.getEnergies().max()<<std::endl;
+	}
+	if (counts.getEMax().max() > psf.getEnergies().back()){
+		std::cout<<"Warning: "<<std::endl;
+		std::cout<<"Higher energy bound of counts is higher than the highest energy of psf, continuing with extrapolation"<<std::endl;
+		std::cout<<counts.getEMax().max()<<" > "<<psf.getEnergies().back()<<std::endl;
+	}
+}
+
+const CountsMap & BaseModel::getCounts() const{
+	return fcounts;
+}
+
+const Skymap<char> & BaseModel::getFilter() const{
+   return ffilter;
+}
+
+std::valarray<double> BaseModel::spectraToCounts(std::valarray<double> &intensities, const std::valarray<double> &energies, const SM::Coordinate &co, bool psfRebin, std::valarray<double> *eMinArr, std::valarray<double> *eMaxArr) const {
+   //We use pointers here because of backward compatibility.  Passing a
+   //reference would be more natural.
+   //If psfRebin is given, make finer energy bins, compatible with psf
+   if (eMinArr == 0) eMinArr = new std::valarray<double>();
+   if (eMaxArr == 0) eMaxArr = new std::valarray<double>();
+   if (psfRebin){
+      getPsfEnergyBoundaries(*eMinArr, *eMaxArr);
+   }else{
+      fcounts.getCountsMap().getBoundaries(*eMinArr, *eMaxArr);
+   }
+
+   if (fconfigure & EXPOSURE_CORRECT) {
+      return fexposure.spectraToCounts(intensities, energies, *eMinArr, *eMaxArr, co);
+   } else {
+      //Do interpolation of intensities
+      std::valarray<double> output(eMinArr->size());
+      int il, iu;
+      for (int i = 0; i < eMinArr->size(); ++i){
+	 Utils::findIndexRange(energies, (*eMinArr)[i], (*eMaxArr)[i], il, iu);
+	 //Calculate the average spectrum over the entire bin
+      	 double flux = 0;
+	 for (int j = il; j < iu; ++j) {
+	    //Handle extrapolation gracefully
+	    double emin;
+      	    if (j == 0 && (*eMinArr)[i] < energies[j]){
+	       emin = (*eMinArr)[i];
+	    } else {
+	       emin = std::max((*eMinArr)[i], energies[j]);
+	    }
+	    double emax;
+	    if (j == energies.size() - 2 && (*eMaxArr)[i] > energies[j+1]){ 
+	       emax = (*eMaxArr)[i];
+	    } else {
+	       emax = std::min((*eMaxArr)[i], energies[j+1]);
+	    }
+	    if (intensities[j] > 0 && intensities[j+1] > 0) {
+	       double tmpInd = log(intensities[j+1]/intensities[j])/log(energies[j+1]/energies[j]);
+	       if (tmpInd != 1) {
+		  flux += intensities[j]*pow(energies[j], -tmpInd)/(tmpInd+1)*(pow(emax,tmpInd+1) - pow(emin,tmpInd+1));
+	       } else {
+		  flux += intensities[j]*pow(energies[j], -tmpInd)*(log(emax)-log(emin));
+	       }
+	    }
+	 } 
+       	 output[i] = flux/((*eMaxArr)[i]-(*eMinArr)[i]);
+	 if (isnan(output[i]) || isinf(output[i])){
+	    std::cout<<"NANs in spectraToCounts "<<intensities[iu]<<", "<<intensities[il]<<", "<<energies[iu]<<", "<<energies[il]<<std::endl;
+	 }
+      }
+      return output;
+   }
+
+}
+
+std::valarray<double> BaseModel::powerLawCounts(double prefactor, double index, double pivot, const SM::Coordinate &co, bool psfRebin, std::valarray<double> *eMinArr, std::valarray<double> *eMaxArr) const{
+   //We use pointers here because of backward compatibility.  Passing a
+   //reference would be more natural.
+   //If psfRebin is given, make finer energy bins, compatible with psf
+   if (eMinArr == 0) eMinArr = new std::valarray<double>();
+   if (eMaxArr == 0) eMaxArr = new std::valarray<double>();
+   if (psfRebin){
+      getPsfEnergyBoundaries(*eMinArr, *eMaxArr);
+   }else{
+      fcounts.getCountsMap().getBoundaries(*eMinArr, *eMaxArr);
+   }
+
+   if (fconfigure & EXPOSURE_CORRECT) {
+      return fexposure.spectraToCounts(index, prefactor, pivot, *eMinArr, *eMaxArr, co);
+   } else {
+      std::valarray<double> sp(eMinArr->size());
+      for (int i = 0; i < sp.size(); ++i){
+	 sp[i] = prefactor * pow(0.5*((*eMinArr)[i]+(*eMaxArr)[i])/pivot, index);
+      }
+      return sp;
+   }
+}
+
+Skymap<double> BaseModel::convolve(const Skymap<double> &skyMap, bool psfRebin) const{
+   Skymap<double> outMap;
+   if (fconfigure & CONVOLVE) {
+      outMap = fpsf.convolve(skyMap);
+   } else {
+      outMap = skyMap;
+   }
+   // Need to re-bin the outMap if psfRebin is true
+   if ( psfRebin ){
+      rebinPsfBoundariesToCounts(outMap);
+   } 
+   return outMap;
+}
+
+Skymap<double> BaseModel::fluxToCounts(const Skymap<double> &inputFluxMap, bool psfRebin) const{
+   // If psfRebin is true, rebin the input bins in such a way to split the bins
+   // in the count map into finer bins so we have very little energy gap for
+   // the psf correction.
+   std::valarray<double> eMinArr, eMaxArr;
+   if ( psfRebin ) {
+      getPsfEnergyBoundaries(eMinArr, eMaxArr);
+   }else{
+      fcounts.getCountsMap().getBoundaries(eMinArr,eMaxArr);
+   }
+   if (fconfigure & EXPOSURE_CORRECT) {
+      //Use the exposure method
+      return fexposure.mapToCounts(inputFluxMap, eMinArr, eMaxArr, fcounts.getCountsMap().Order());
+   } else {
+      //Rebin the map to conform to counts (eMinArr)
+      Skymap<double> output(fcounts.getCountsMap().Order(), eMinArr, eMaxArr);
+      Skymap<double> rebinned = inputFluxMap.rebin(fcounts.getCountsMap().Order());
+      for (int k = 0; k < output.Npix(); ++k) {
+	 int il, iu;
+	 for (int i = 0; i < eMinArr.size(); ++i){
+	    Utils::findIndexRange(inputFluxMap.getSpectra(), eMinArr[i], eMaxArr[i], il, iu);
+	    //Calculate the average spectrum over the entire bin
+	    double flux = 0;
+	    for (int j = il; j < iu; ++j) {
+	       //Handle extrapolation gracefully
+	       double emin;
+	       if (j == 0 && eMinArr[i] < inputFluxMap.getSpectra()[j]){
+		  emin = eMinArr[i];
+	       } else {
+		  emin = std::max(eMinArr[i], inputFluxMap.getSpectra()[j]);
+	       }
+	       double emax;
+	       if (j == inputFluxMap.nSpectra() - 2 && eMaxArr[i] > inputFluxMap.getSpectra()[j+1]){ 
+		  emax = eMaxArr[i];
+	       } else {
+		  emax = std::min(eMaxArr[i], inputFluxMap.getSpectra()[j+1]);
+	       }
+	       if (rebinned[k][j] > 0 && rebinned[k][j+1] > 0) {
+		  double tmpInd = log(rebinned[k][j+1]/rebinned[k][j])/log(inputFluxMap.getSpectra()[j+1]/inputFluxMap.getSpectra()[j]);
+		  if (tmpInd != 1) {
+		     flux += rebinned[k][j]*pow(inputFluxMap.getSpectra()[j], -tmpInd)/(tmpInd+1)*(pow(emax,tmpInd+1) - pow(emin,tmpInd+1));
+		  } else {
+		     flux += rebinned[k][j]*pow(inputFluxMap.getSpectra()[j], -tmpInd)*(log(emax)-log(emin));
+		  }
+	       }
+	    } 
+	    output[k][i] = flux/(eMaxArr[i]-eMinArr[i]);
+	    if (isnan(output[k][i]) || isinf(output[k][i]) ) {
+	       std::cout<<"NAN in fluxToCounts "<<rebinned[k][iu]<<", "<<rebinned[k][il]<<", "<<inputFluxMap.getSpectra()[iu]<<", "<<inputFluxMap.getSpectra()[il]<<std::endl;
+	    }
+	 }
+      }
+      return output;
+   }
+}
+
+Variables & BaseModel::getVariables(){
+	return fvariables;
+}
+const Variables & BaseModel::getVariables() const{
+	return fvariables;
+}
+
+void BaseModel::setVariables(const Variables &vars){
+	//Get the names of the current variables
+	std::vector<std::string> varNames = fvariables.getNames();
+	for (int i = 0; i < varNames.size(); ++i) {
+		fvariables[varNames[i]] = vars[varNames[i]];
+		fvariables.error(varNames[i]) = vars.error(varNames[i]);
+	}
+}
+
+Sources BaseModel::getSources(const Variables & vars) const {
+	Sources sources(fsources);
+	modifySources(vars, sources);
+	return sources;
+}
+
+void BaseModel::getPsfEnergyBoundaries(std::valarray<double> &eMin, std::valarray<double> &eMax) const{
+   std::vector<double> eMinVec, eMaxVec;
+   // First we need to get the energies at which the psf is defined
+   const std::vector<double> &psfE = fpsf.getEnergies();
+   // Then we get the boundaries from the counts map
+   const std::valarray<double> &cMin = fcounts.getEMin();
+   const std::valarray<double> &cMax = fcounts.getEMax();
+   // Now we need to loop over the counts energy boundaries and add the
+   // necessary splits to the energy arrays.  We split such that there is
+   // at least one psf defined in an energy bin, within the limits of the
+   // counts binning of course.
+   for (size_t i = 0; i < cMin.size(); ++i){
+      // Get the indexes of the psf within the energy bin
+      int lowIndex, highIndex;
+      // The routine brackets the energy bin, unless the psf energy range
+      // does not span the energy bin.
+      Utils::findIndexRange(psfE, cMin[i], cMax[i], lowIndex, highIndex);
+
+      // Unless the psf energy range does not extend beyond the minimum
+      // boundary, skip the first element
+      int j = 1; 
+      if ( psfE[lowIndex] > cMin[i] ) j = 0;
+      // Same for the last element
+      int k = 1;
+      if ( psfE[highIndex] < cMax[i] ) k = 0;
+      eMinVec.push_back(cMin[i]);
+      for (int l = lowIndex+j; l < highIndex-k; ++l){
+	 eMaxVec.push_back((psfE[l]+psfE[l+1])/2.);
+	 eMinVec.push_back((psfE[l]+psfE[l+1])/2.);
+      }
+      eMaxVec.push_back(cMax[i]);
+   }
+   eMin.resize(eMinVec.size());
+   eMax.resize(eMaxVec.size());
+   std::copy(eMinVec.begin(), eMinVec.end(), &eMin[0]);
+   std::copy(eMaxVec.begin(), eMaxVec.end(), &eMax[0]);
+}
+
+
+void BaseModel::rebinPsfBoundariesToCounts(Skymap<double> & inMap) const{
+   if (fcounts.getEMax().size() < inMap.nSpectra()) {
+      Skymap<double> tmpMap(fcounts.getCountsMap().Order(), fcounts.getEMin(), fcounts.getEMax(), RING);
+      //Get the boundaries of the input map
+      std::valarray<double> eMinArr, eMaxArr;
+      inMap.getBoundaries(eMinArr, eMaxArr);
+      //std::copy(&eMaxArr[0],&eMaxArr[0]+eMaxArr.size(), std::ostream_iterator<double>(std::cout, ", "));
+      //Loop over the count map bins and sum up the output map
+      int j = 0;
+      for (int i = 0; i < fcounts.getEMin().size(); ++i){
+	 if (fcounts.getEMin()[i] != eMinArr[j]) {
+	    std::cerr<<"Things might be wrong when rebinning in convolution.  Energy boundaries don't match"<<std::endl;
+	    std::cerr<<fcounts.getEMin()[i]<<" != "<<eMinArr[j]<<std::endl;
+	 }
+	 int cnt = 0;
+	 while (eMaxArr[j] <= fcounts.getEMax()[i] && j < eMaxArr.size()) {
+	    //std::cout<<j<<", ";
+	    for (int k = 0; k < inMap.Npix(); ++k){
+	       tmpMap[k][i] += inMap[k][j];
+	    }
+	    ++j;
+	    ++cnt;
+	 }
+	 if (!(fconfigure & EXPOSURE_CORRECT) && cnt != 0) {
+	    for (int k = 0; k < inMap.Npix(); ++k) {
+	       tmpMap[k][i] /= float(cnt);
+	    }
+	 }
+	 if (fcounts.getEMax()[i] != eMaxArr[j-1]){
+	    std::cerr<<"Things might be wrong when rebinning in convolution.  Energy boundaries don't match"<<std::endl;
+	    std::cerr<<fcounts.getEMax()[i]<<" != "<<eMaxArr[j-1]<<std::endl;
+	 }
+      }
+      inMap = tmpMap;
+   } 
+ }

source/CAT.cc

+using namespace std;
+// for slalib version which uses: if defined (__cplusplus) extern "C"
+// but this does not work, use version without this 
+// #define __cplusplus
+
+#include<iostream>
+#include"fitsio.h"
+#include"slalib.h"
+
+#include"CAT.h"
+
+
+
+//////////////////////////////////////////////////////
+  
+  
+///////////////////////////////////////////////////////
+  int CAT:: read (char *directory)
+{
+  char infile[500];
+
+  char err_text[100];
+  char comment[100];
+
+  fitsfile *fptr;
+  int status;
+  int hdunum,hdutype,total_hdus,colnum;
+  int n_columns,n_rows;
+  double nulval=0.;
+  int anynul   =0;
+  char nulstr[100];strcpy(nulstr,"null");
+  double dtr=acos(-1.0)/180.;
+  int i,j,k;
+
+
+  //char catalogue_file[]="gnrl_refr_cat_0020_reformatted.fits";
+  //char catalogue_file[]="isgriCat_2nd_Tony.fits";
+    char catalogue_file[]="isgri_2nd_cat_revised.fits";
+ 
+
+
+
+  cout<< "CAT::read  (" << directory<<")" <<endl;
+
+  
+  strcpy(infile,directory);
+  strcat(infile,catalogue_file);
+
+  
+  cout<<"   infile "<<infile<<endl;
+
+  status=0; // needed
+  fits_open_file(&fptr,infile,READONLY,&status) ;
+  cout<<"FITS read open status= "<<status<<" "<<infile<<endl;
+
+  fits_get_errstatus(status,err_text);
+  cout<<err_text<<endl;
+
+fits_get_num_hdus(fptr,&total_hdus,&status);
+  cout<<"total number of header units="<<total_hdus<<endl;
+
+  hdunum=1;
+
+  fits_movabs_hdu(fptr,hdunum+1,&hdutype,&status );
+  cout<<"FITS movabs_hdu hdunum="<<hdunum
+      <<" status= "<<status<<" hdutype="<<hdutype<<endl  ;
+
+   fits_read_key(fptr,TLONG,"TFIELDS",&n_columns,comment,&status);
+   cout<<" n_columns ="<<n_columns ;
+   fits_read_key(fptr,TLONG,"NAXIS2",&n_rows  ,comment,&status);
+   cout<<" n_rows    ="<<n_rows    <<endl;
+
+
+   n_sources=n_rows;
+
+   ra      =new double[n_rows];
+   dec     =new double[n_rows];
+   flux    =new double[n_rows];//AWS20060110
+   flux_err=new double[n_rows];//AWS20060110
+
+  fits_get_colnum(fptr,CASEINSEN,"RA_OBJ",&colnum,&status);
+  cout<<"RA_OBJ column ="<< colnum<<endl;
+  fits_read_col(fptr, TDOUBLE, colnum, 1,1,n_rows,&nulval,ra, &anynul,&status);
+
+  fits_get_colnum(fptr,CASEINSEN,"DEC_OBJ",&colnum,&status);
+  cout<<"DEC_OBJ column ="<< colnum<<endl;
+  fits_read_col(fptr, TDOUBLE, colnum, 1,1,n_rows,&nulval,dec, &anynul,&status);
+
+
+  // cout<<"ra :";for (i=0;i<n_rows;i++)cout<<" "<<ra [i];cout<<endl;
+  //cout<<"dec:";for (i=0;i<n_rows;i++)cout<<" "<<dec[i];cout<<endl;
+
+
+ fits_get_colnum(fptr,CASEINSEN,"NAME",&colnum,&status);
+ cout<<"NAME column ="<< colnum<<endl;
+ name=new char*[n_rows];
+ for (i=0;i<n_rows;i++)     name[i]=new char[100];
+ fits_read_col_str(fptr,  colnum, 1,1, n_rows,nulstr,  name, &anynul,&status);
+
+ // for (i=0;i<n_rows;i++)  cout<<   name[i] <<" ra="<<ra[i]<<" dec="<<dec[i]<<endl;
+
+  longitude =new double[n_rows];
+  latitude  =new double[n_rows];
+
+ for (i=0;i<n_rows;i++) 
+   {
+     // NB had to compile sla_c since undefined symbol using lib from osa4.2
+
+  slaEqgal(ra[i]*dtr, dec[i]*dtr,  &longitude[i],  &latitude[i]);
+  longitude[i]/=dtr;
+   latitude[i]/=dtr;
+   }
+
+
+
+  return 0;
+}
+///////////////////////////////////////////////////////
+  int CAT:: read (char *directory,char *filename) //AWS20060110
+{
+  char infile[500];
+
+  char err_text[100];
+  char comment[100];
+
+  fitsfile *fptr;
+  int status;
+  int hdunum,hdutype,total_hdus,colnum;
+  int n_columns,n_rows;
+  double nulval=0.;
+  int anynul   =0;
+  char nulstr[100];strcpy(nulstr,"null");
+  double dtr=acos(-1.0)/180.;
+  int i,j,k;
+
+
+  
+
+
+  cout<< "CAT::read  (" << directory<<"/"<<filename<<")" <<endl;
+
+  
+  strcpy(infile,directory);
+  strcat(infile,filename);
+
+  
+  cout<<"CAT::read   infile= "<<infile<<endl;
+
+  status=0; // needed
+  fits_open_file(&fptr,infile,READONLY,&status) ;
+  cout<<"FITS read open status= "<<status<<" "<<infile<<endl;
+
+  fits_get_errstatus(status,err_text);
+  cout<<err_text<<endl;
+
+  fits_get_num_hdus(fptr,&total_hdus,&status);
+  cout<<"total number of header units="<<total_hdus<<endl;
+
+  hdunum=total_hdus-1; // because spimodfit output has also grouping extension, make it thus more general assuming catalogue is last extension
+
+  fits_movabs_hdu(fptr,hdunum+1,&hdutype,&status );
+  cout<<"FITS movabs_hdu hdunum="<<hdunum
+      <<" status= "<<status<<" hdutype="<<hdutype<<endl  ;
+
+   fits_read_key(fptr,TLONG,"TFIELDS",&n_columns,comment,&status);
+   cout<<" n_columns ="<<n_columns ;
+   fits_read_key(fptr,TLONG,"NAXIS2",&n_rows  ,comment,&status);
+   cout<<" n_rows    ="<<n_rows    <<endl;
+
+
+   n_sources=n_rows;
+
+  ra      =new double[n_rows];
+  dec     =new double[n_rows];
+  flux    =new double[n_rows];
+