Commit ad938765 authored by Maude Le Jeune's avatar Maude Le Jeune
Browse files

mv gauss fit in dedicated header

parent 88a5aba4
// Copyright (C) 2007,2008 Marc Betoule
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// any later version.
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public License for more details.
// You should have received a copy of the GNU General Public License
// along with this program. If not, see <http://www.gnu.org/licenses/>.
// report bugs at betoule@apc.univ-paris7.fr
#include "healpix_base.h"
#include "Python.h"
#include "MapExt.h"
#include "rotmatrix.h"
#include "pointing.h"
#include <string>
#include <iostream>
#include <iterator>
#include "atlasmanager.h"
#include "catalog.h"
#include "math.h"
#include <sys/time.h>
using namespace std;
void gaussfit (Healpix_Map<double> &hmap, Healpix_Map<double> &smap, double* tabcat, int ns,int ncol, double fwhm_arcmin,double nfwhm ,double latcut,double maxxi2 ) {
bool sigmap = true;
Catalog cat(NULL,0); // todo load tabcat in cat
//# pix, glon, glat, A, sigmaA
//# theta phi A sigmaTheta sigmaPhi sigmaA chi2/dof
//cat.input(tabcat, ns, ncol);
// list<Source>::iterator s;
double chisqrdof;
double fwhm = fwhm_arcmin * degr2rad /60.;
double radius = nfwhm;
struct timeval startTime;
struct timeval endTime;
#pragma omp parallel
{
#pragma omp for
for (int i=0; i<ns; i++){
// for(s = cat.begin(); s != cat.end(); s++){
// pointing p = (*s).getPointing();
int icol = i*ncol;
double guess[3];
double err[3];
pointing p;
p.phi = tabcat[icol+1]* degr2rad;
p.theta = tabcat[icol+2];
p.theta = (90. - p.theta) * degr2rad;
guess[1] = p.theta;
guess[2] = p.phi;
guess[0] = tabcat[icol+3];
if (abs(90-(p.theta*180.0/pi)) < latcut)
chisqrdof = INFINITY;
else {
healpixGaussianFit(hmap, fwhm, guess, err, &chisqrdof,smap,radius);
}
// fit is ok, set to new value
if (chisqrdof < maxxi2 && guess[0]>0){
addGaussianPatternOnSky(hmap, pointing(guess[1],guess[2]), fwhm, guess[0], err[0]/10, true);
tabcat[icol+0] = guess[1];
tabcat[icol+1] = guess[2];
tabcat[icol+2] = guess[0];
tabcat[icol+3] = err[1];
tabcat[icol+4] = err[2];
tabcat[icol+5] = err[0];
tabcat[icol+6] = chisqrdof;
}
// fit is bad, keep old values
else {
tabcat[icol+0] = p.theta;
tabcat[icol+1] = p.phi;
tabcat[icol+2] = tabcat[icol+3];
tabcat[icol+3] = 0;
tabcat[icol+5] = tabcat[icol+4];
tabcat[icol+4] = 0;
tabcat[icol+6] = chisqrdof;
}
//i++;
}
}
}
......@@ -108,6 +108,9 @@ ostream &operator<< (ostream &os, const Catalog & catalog){
return os;
}
// void Catalog::output(ostream &os, double freq, string unit){
// os << "# pixel lon-deg lat-deg flux (Jy) fluxerror\n";
// double conversionFactor = (Unit("Jy/sr")/Unit(unit))(freq);
......
......@@ -315,6 +315,7 @@ public:
void input(istream &is);
void input(arr2<double> catarr);
void input(double* tabcat, int ns, int ncol);
private:
double nsigma;
AtlasManager * sky;
......
......@@ -487,7 +487,7 @@ void healpixGaussianFit(Healpix_Map<double> & m, double fwhm, double * guess, do
vector< int > listpix;
m.query_disc (p, radius, listpix);
d.n = listpix.size();
//cout << d.n <<endl;
//cout << "number of pixel = " << d.n <<endl;
d.sigma2 = fwhm*fwhm / (8. * M_LN2);
d.y = new double[d.n];
d.sigma = new double[d.n];
......@@ -529,6 +529,7 @@ void healpixGaussianFit(Healpix_Map<double> & m, double fwhm, double * guess, do
//fit
int status = GSL_CONTINUE;
for(int iter = 0;status == GSL_CONTINUE && iter < 500 ; iter++){
status = gsl_multifit_fdfsolver_iterate (s);
//printf ("status = %s\n", gsl_strerror (status));
......
......@@ -31,6 +31,8 @@
#include "atlasmanager.h"
#include "catalog.h"
#include "math.h"
#include <sys/time.h>
#include "psutils.h"
using namespace std;
......@@ -365,59 +367,11 @@ void _psgaussfit (double* tabmap, int npix, double* sigmamap, int npix2, double*
arr<double> arrmap(tabmap, npix);
MapExt<double> hmap(arrmap, RING);
// -------- GAUSSIAN FIT
bool sigmap = true;
arr<double> arrmap2(sigmamap, npix);
Healpix_Map <double> smap(arrmap2, RING);
sigmap = true;
Catalog cat(NULL,0); // todo load tabcat in cat
//# pix, glon, glat, A, sigmaA
//# theta phi A sigmaTheta sigmaPhi sigmaA chi2/dof
cat.input(tabcat, ns, ncol);
list<Source>::iterator s;
double guess[3];
double err[3];
double chisqrdof;
double fwhm = fwhm_arcmin * degr2rad /60.;
double radius = nfwhm;
int i=0;
gaussfit (hmap, smap, tabcat, ns, ncol, fwhm_arcmin, nfwhm , latcut, maxxi2 );
for(s = cat.begin(); s != cat.end(); s++){
pointing p = (*s).getPointing();
guess[1] = p.theta;
guess[2] = p.phi;
guess[0] = (*s).getAmplitude();
if (abs(90-(p.theta*180.0/pi)) < latcut)
chisqrdof = INFINITY;
else
healpixGaussianFit(hmap, fwhm, guess, err, &chisqrdof,smap,radius);
// fit is ok, set to new value
if (chisqrdof < maxxi2 && guess[0]>0){
addGaussianPatternOnSky(hmap, pointing(guess[1],guess[2]), fwhm, guess[0], err[0]/10, true);
tabcat[i*ncol+0] = guess[1];
tabcat[i*ncol+1] = guess[2];
tabcat[i*ncol+2] = guess[0];
tabcat[i*ncol+3] = err[1];
tabcat[i*ncol+4] = err[2];
tabcat[i*ncol+5] = err[0];
tabcat[i*ncol+6] = chisqrdof;
}
// fit is bad, keep old values
else {
tabcat[i*ncol+0] = p.theta;
tabcat[i*ncol+1] = p.phi;
tabcat[i*ncol+2] = (*s).getAmplitude();
tabcat[i*ncol+3] = 0;
tabcat[i*ncol+4] = 0;
tabcat[i*ncol+5] = (*s).getSigma();
tabcat[i*ncol+6] = chisqrdof;
}
i++;
}
arrmap = hmap.Map();
arrmap2.dealloc();
}
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment