Commit 76ddc292 authored by Danny Boxhoorn's avatar Danny Boxhoorn
Browse files

Version from 2006/11/16

parent ec780632
#include <stdio.h>
#include <math.h>
void base_func(float sigma, int m, int n, double *ker_x, double *ker_y,
int kerhw, int do_norm)
{
double sig, norm_x, norm_y ;
int i, idx;
sig = (double)(2*sigma*sigma);
norm_x = norm_y = 0.0;
for(i=-kerhw;i<=kerhw;i++) {
idx = kerhw + i;
ker_y[idx] = ker_x[idx] = exp(-i*i/sig);
if( m > 0 )
ker_x[idx] *= pow(i,m);
if( n > 0 )
ker_y[idx] *= pow(i,n);
norm_x += ker_x[idx];
norm_y += ker_y[idx];
}
if( do_norm ) {
if( norm_x != 0.0 && norm_y != 0.0 ) {
for(i=0;i<=2*kerhw;i++) {
ker_x[i] /= norm_x;
ker_y[i] /= norm_y;
}
}
else {
fprintf(stderr,"base_func: warning: normalization to 0.0 ignored!\n");
}
return ;
}
return ;
}
#include <stdlib.h>
#include "defs.h"
#include "funcs.h"
/*
Given image im(nx,ny) and centroid (cx,cy) estimates background
taking data from annulus (r1,r2). Sigma clipping with n_sig until
nothing changes
*/
float bkg(float *im, int cx, int cy, PAR_STRUCT *par)
{
int i, j, npix, ir, *idx, delta, nx, ny;
float r1, r2, rr, rr1, rr2, bglev, absdev, mean, *tab, n_sig;
n_sig = par->nsig_bkg;
r1 = par->anrad1;
r2 = par->anrad2;
nx = par->nx0;
ny = par->ny0;
npix = 3.14*((r2+1.0)*(r2+1.0) - r1*r1);
tab = (float *) malloc(npix*sizeof(float));
idx = (int *) malloc(npix*sizeof(int));
ir = r2 + 1;
rr1 = r1*r1;
rr2 = r2*r2;
/* load pixels within annulus (r1,r2) into tab */
npix = 0;
for(j=-ir;j<=ir;j++) {
for(i=-ir;i<=ir;i++) {
rr = i*i + j*j;
if( rr1<rr && rr<rr2 )
tab[npix++] = im[cx+i+nx*(cy+j)];
}
}
/*
take median as first background level
use little trick to fool stupid 1,n indexing of numerical recipes
*/
indexx(npix,tab-1,idx-1);
bglev = tab[idx[npix/2]-1];
/*
reject pixels deviating by more than n_sig*sigma from current bglev
subsequent bglev will be means of the remaining pixels
quit when no more rejections
*/
delta = 1;
while(delta>0) {
absdev=0.0;
for(i=0;i<npix;i++)
absdev += (bglev - tab[i])*(bglev - tab[i]);
absdev *= n_sig*n_sig/(float)npix;
j=0;
mean=0.0;
for(i=0;i<npix;i++) {
if((bglev - tab[i])*(bglev - tab[i]) < absdev) {
tab[j++] = tab[i];
mean += tab[i];
}
}
if(j>1) {
delta = npix - j;
npix = j - 1;
mean /= npix;
bglev = mean;
}
else
return 0.0;
}
free(tab);
free(idx);
return (bglev);
}
typedef struct {
int nx0;
int ny0;
int psfhw;
int psfn;
int bkg_mode;
int verbose;
float aprad;
float fitrad;
float normrad;
float anrad1;
float anrad2;
float nsig_bkg;
float sat_level;
float min_level;
float bad_value;
float gain;
float exptime;
float err;
float readnoise;
} PAR_STRUCT;
typedef struct {
float p_flux;
float a_flux;
float p2_flux;
float a2_flux;
float p_err;
float a_err;
float p2_err;
float a2_err;
float corr;
float chi2_n;
float ker_chi2;
float bg;
float bg2;
float fwhm;
int nbad;
} STAR;
typedef struct
{
int hw;
int ndeg_spat;
int ndeg_local;
int ngauss;
int recenter;
double *buf;
double cos;
double sin;
double ax;
double ay;
double sigma_inc;
double sigma_mscale;
double *vec;
double *local_vec;
double fitrad;
double normrad;
double gain;
double x_orig;
double y_orig;
} PSF_STRUCT;
typedef struct {
int nsec_x;
int nsec_y;
int nx;
int ny;
int hw;
int wdeg;
int ncomp;
int nvecs;
int ntot;
int nwxy;
float *sig;
int *deg;
double *chi2_n;
double **vec;
double **vecs;
} KER_STRUCT;
extern void get_phot(float *, float *, float *, double, double, int, int,
double *, STAR *, PAR_STRUCT *);
extern void get_params(char *, PAR_STRUCT *);
extern void read_psf(char *, PSF_STRUCT *);
extern void init_psf(PSF_STRUCT *, double, double);
extern double psf_core(PSF_STRUCT *, double, double);
extern void read_kernel(char *, KER_STRUCT *, int);
extern void spatial_coeffs(KER_STRUCT *, double, double, double *);
extern void make_vectors(KER_STRUCT *);
extern void make_kernel(KER_STRUCT *, double *, double *, int);
extern void im_convolve(double *, double *, int, int, double *, int);
extern float *rfits_f(char *, char *, int *, int *,
float *, float *, double *, double *, long int *);
extern void wfits_f(char *, char *, long int, float *);
extern void make_psf(PSF_STRUCT *, double, double, int, int, double *);
extern float bkg(float *, int, int, PAR_STRUCT *);
extern void indexx(int, float *, int *);
extern int get_sector(double *, double *, double *, double *, KER_STRUCT *);
extern void base_func(float, int, int, double *, double *, int, int);
extern double get_fwhm(double *, double, double, int, int,
PAR_STRUCT *, double *);
#include <math.h>
#include "defs.h"
double get_fwhm(double *psfim, double x, double y, int cx, int cy,
PAR_STRUCT *par, double *ratio)
{
int i, j, kx, ky, irad, psfn, psfhw;
double rr, rad2, psf_pix, sxx, sxy, syy, sum, delta;
double sigma_x, sigma_y, dx, dy, fwhm;
rad2 = par->fitrad*par->fitrad;
irad = (int)par->fitrad + 2;
psfhw = par->psfhw;
psfn = par->psfn;
sxx = sxy = syy = sum = 0.0;
for(i=cx-irad;i<=cx+irad;i++) {
for(j=cy-irad;j<=cy+irad;j++) {
kx = i - cx;
ky = j - cy;
dx = i - x;
dy = j - y;
rr = dx*dx + dy*dy;
psf_pix = psfim[psfhw+kx+psfn*(psfhw+ky)];
if( rr <= rad2 ) {
sxx += psf_pix*dx*dx;
sxy += psf_pix*dx*dy;
syy += psf_pix*dy*dy;
sum += psf_pix;
}
}
}
if( sum != 0.0) {
sxx /= sum;
sxy /= sum;
syy /= sum;
}
delta = sxx*sxx - 2.0*sxx*syy + syy*syy + 4.0*sxy*sxy;
sigma_x = (sxx + syy + sqrt(delta))*0.5;
sigma_y = (sxx + syy - sqrt(delta))*0.5;
if( sigma_y > 0.0 && sigma_x > 0.0 ) {
if( sigma_y > sigma_x )
*ratio = sqrt(sigma_x/sigma_y);
else
*ratio = sqrt(sigma_y/sigma_x);
fwhm = sqrt(4.0*log(2.0)*(sigma_x + sigma_y));
if( fwhm > 5.0 )
fwhm = 5.0;
} else {
*ratio = 0.0;
fwhm = 5.0;
}
return (fwhm);
}
#include <stdio.h>
#include "defs.h"
void get_params(char *parfname, PAR_STRUCT *par)
{
FILE *inpf;
inpf = fopen(parfname,"r");
fscanf(inpf,"%*s = %f %*s = %f",&par->aprad,&par->fitrad);
fscanf(inpf,"%*s = %f %*s = %f",&par->normrad,&par->anrad1);
fscanf(inpf,"%*s = %f %*s = %d",&par->anrad2,&par->bkg_mode);
fscanf(inpf,"%*s = %f %*s = %f",&par->nsig_bkg,&par->sat_level);
fscanf(inpf,"%*s = %f %*s = %f",&par->min_level,&par->bad_value);
fscanf(inpf,"%*s = %f %*s = %f",&par->gain,&par->exptime);
fscanf(inpf,"%*s = %f %*s = %d",&par->err,&par->verbose);
fscanf(inpf,"%*s = %f\n",&par->readnoise);
fclose(inpf);
return;
}
#include <math.h>
#include "defs.h"
void get_phot(float *difim, float *im, float *refim, double x, double y,
int cx, int cy, double *psfim, STAR *obj, PAR_STRUCT *par)
{
int i, j, kx, ky, irad, count, nbad, psfn, psfhw, nx0, idx;
double rr, norm, p_sig2, a_sig2, p_phot, a_phot,p2_phot, a2_phot;
double a, b, c, sum_a, sum_b, sum_c, p_sum2, rad2;
double bg,bg2, sat, min, gain, time, pix, pix_psf, weight, pix_ref;
double readnoise2=0;
rad2 = par->fitrad*par->fitrad;
gain = par->gain;
time = par->exptime;
sat = par->sat_level;
min = par->min_level;
bg = obj->bg;
bg2 = obj->bg2;
readnoise2 = (par->readnoise*par->readnoise);
irad = (int)par->fitrad + 2;
psfhw = par->psfhw;
psfn = par->psfn;
nx0 = par->nx0;
count = 0;
nbad = 0;
sum_a = 0.0;
sum_b = 0.0;
sum_c = 0.0;
p_sum2 = 0.0;
p_phot = 0.0;
p_sig2 = 0.0;
a_phot = 0.0;
a_sig2 = 0.0;
norm = 0.0;
p2_phot = 0.0;
a2_phot = 0.0;
for(i=0;i<psfn*psfn;i++)
norm += psfim[i];
for(i=cx-irad;i<=cx+irad;i++) {
for(j=cy-irad;j<=cy+irad;j++) {
kx = i - cx;
ky = j - cy;
rr = (i-x)*(i-x) + (j-y)*(j-y);
idx = i + nx0*j;
pix_psf = psfim[psfhw+kx+psfn*(psfhw+ky)];
pix = difim[idx];
weight = im[idx] + (readnoise2*gain);
pix_ref = refim[idx];
if( rr <= rad2 ) {
if( pix != par->err && weight<sat && pix_ref<sat &&
weight>=min && pix_ref>=min ) {
a = pix_psf*(pix - bg);
b = pix_psf*pix_psf;
c = (pix - bg)*(pix - bg);
sum_a += a;
sum_b += b;
sum_c += c;
a_phot += pix - bg;
a_sig2 += weight;
p_phot += a/(weight*norm);
p_sig2 += b/(weight*norm*norm);
p_sum2 += c/weight;
count++;
}
else
nbad++;
}
}
}
if( p_sig2 != 0.0 && count > 0 ) {
obj->p_flux = (float)(p_phot/(p_sig2*time));
obj->p_err = (float)(1.0/(sqrt(p_sig2*gain)*time));
obj->chi2_n = (float)((p_sum2 - p_phot*p_phot/p_sig2)*gain/count);
obj->corr = (float)(sum_a/sqrt(sum_b*sum_c));
obj->a_flux = (float)(a_phot/time);
obj->a_err = (float)(sqrt(a_sig2/gain)/time);
obj->nbad = nbad;
/* obj->a2_flux = (float)(a2_phot/time); */
/* obj->a2_err = (float)(sqrt(a_sig2/gain)/time); */
/* obj->p2_flux = (float)(p2_phot/(p_sig2*time)); */
/* obj->p2_err = (float)(1.0/(sqrt(p_sig2*gain)*time)); */
}
else {
obj->p_flux = par->bad_value;
obj->p_err = par->bad_value;
obj->chi2_n = par->bad_value;
obj->corr = par->bad_value;
obj->a_flux = par->bad_value;
obj->a_err = par->bad_value;
obj->nbad = -1;
}
return;
}
void im_convolve(double *im, double *cnvim, int nx, int ny,
double *kerim, int kerhw)
{
int i, j, m, n, kern, idx;
kern = 2*kerhw + 1;
for(j=kerhw;j<ny-kerhw;j++) {
for(i=kerhw;i<nx-kerhw;i++) {
idx = i + nx*j;
cnvim[idx] = 0.0;
for(n=-kerhw;n<=kerhw;n++)
for(m=-kerhw;m<=kerhw;m++)
cnvim[idx] += im[m+idx+nx*n]*kerim[kerhw-m+kern*(kerhw-n)];
}
}
return;
}
void indexx(int n, float *arrin, int *indx)
{
int l,j,ir,indxt,i;
float q;
for (j=1;j<=n;j++) indx[j]=j;
l=(n >> 1) + 1;
ir=n;
for (;;)
{
if (l > 1)
q=arrin[(indxt=indx[--l])];
else
{
q=arrin[(indxt=indx[ir])];
indx[ir]=indx[1];
if (--ir == 1)
{
indx[1]=indxt;
return;
}
}
i=l;
j=l << 1;
while (j <= ir && j>0)
{
if (j < ir && arrin[indx[j]] < arrin[indx[j+1]]) j++;
if (q < arrin[indx[j]])
{
indx[i]=indx[j];
j += (i=j);
if(i<1 || i>(n+1) || j<1 || j>(n+1)) break;
}
else j=ir+1;
}
indx[i]=indxt;
}
}
<
#include<stdio.h>
#include<math.h>
#include "defs.h"
void init_psf(PSF_STRUCT *psf, double xpsf, double ypsf)
{
int m, n, icomp, ncomp, itot, sdeg, ldeg;
double a1, a2, *vec, *local_vec, x_orig, y_orig;
x_orig = psf->x_orig;
y_orig = psf->y_orig;
sdeg = psf->ndeg_spat;
ldeg = psf->ndeg_local;
ncomp = psf->ngauss*(ldeg+1)*(ldeg+2)/2;
local_vec = psf->local_vec;
vec = psf->vec;
for(icomp=0;icomp<ncomp;icomp++)
local_vec[icomp] = 0.0;
itot = 0;
a1 = 1.0;
for(m=0;m<=sdeg;m++) {
a2 = 1.0;
for(n=0;n<=sdeg-m;n++) {