Misc changes
Cees Bassa
parent
6a9a8514bf
commit
83ca5c21e0
|
|
@ -0,0 +1,133 @@
|
|||
#include <stdlib.h>
|
||||
#include <stdio.h>
|
||||
#include <gsl/gsl_rng.h>
|
||||
#include <gsl/gsl_randist.h>
|
||||
#include <gsl/gsl_matrix.h>
|
||||
#include <gsl/gsl_vector.h>
|
||||
#include <gsl/gsl_blas.h>
|
||||
#include <gsl/gsl_multifit_nlinear.h>
|
||||
|
||||
struct data {
|
||||
size_t n;
|
||||
double *x,*y,*sy,*w;
|
||||
};
|
||||
|
||||
int gauss_f(const gsl_vector *q,void *d,gsl_vector *f)
|
||||
{
|
||||
int i,n=((struct data *) d)->n;
|
||||
double *x=((struct data *) d)->x;
|
||||
double *y=((struct data *) d)->y;
|
||||
double a[5],ym,arg,ex,fac;
|
||||
|
||||
// Extract parameters
|
||||
for (i=0;i<5;i++)
|
||||
a[i]=gsl_vector_get(q,i);
|
||||
|
||||
// Compute values
|
||||
for (i=0;i<n;i++) {
|
||||
arg=(x[i]-a[0])/a[1];
|
||||
ex=exp(-arg*arg);
|
||||
fac=2.0*a[2]*ex*arg;
|
||||
ym=a[2]*ex+a[3]+a[4]*x[i];
|
||||
gsl_vector_set(f,i,ym-y[i]);
|
||||
}
|
||||
|
||||
return GSL_SUCCESS;
|
||||
}
|
||||
|
||||
int gauss_df(const gsl_vector *q,void *d,gsl_matrix *J)
|
||||
{
|
||||
int i,n=((struct data *) d)->n;
|
||||
double *x=((struct data *) d)->x;
|
||||
double *y=((struct data *) d)->y;
|
||||
double a[5],ym,arg,ex,fac;
|
||||
|
||||
// Extract parameters
|
||||
for (i=0;i<5;i++)
|
||||
a[i]=gsl_vector_get(q,i);
|
||||
|
||||
// Compute values
|
||||
for (i=0;i<n;i++) {
|
||||
arg=(x[i]-a[0])/a[1];
|
||||
ex=exp(-arg*arg);
|
||||
fac=2.0*a[2]*ex*arg;
|
||||
gsl_matrix_set(J,i,0,fac/a[1]);
|
||||
gsl_matrix_set(J,i,1,fac*arg/a[1]);
|
||||
gsl_matrix_set(J,i,2,ex);
|
||||
gsl_matrix_set(J,i,3,1.0);
|
||||
gsl_matrix_set(J,i,4,x[i]);
|
||||
}
|
||||
|
||||
return GSL_SUCCESS;
|
||||
}
|
||||
|
||||
int fit_gaussian(double *x,double *y,double *sy,int n,double *q,double *sq,int m,double *chisq)
|
||||
{
|
||||
int i;
|
||||
const gsl_multifit_nlinear_type *T=gsl_multifit_nlinear_trust;
|
||||
gsl_multifit_nlinear_workspace *w;
|
||||
gsl_multifit_nlinear_fdf fdf;
|
||||
gsl_multifit_nlinear_parameters fdf_params=gsl_multifit_nlinear_default_parameters();
|
||||
gsl_vector *f;
|
||||
gsl_matrix *J;
|
||||
struct data d;
|
||||
gsl_rng *r;
|
||||
int status,info;
|
||||
const double xtol=1e-8;
|
||||
const double gtol=1e-8;
|
||||
const double ftol=0.0;
|
||||
|
||||
// Fill struct
|
||||
d.n=n;
|
||||
d.x=(double *) malloc(sizeof(double)*d.n);
|
||||
d.y=(double *) malloc(sizeof(double)*d.n);
|
||||
d.sy=(double *) malloc(sizeof(double)*d.n);
|
||||
d.w=(double *) malloc(sizeof(double)*d.n);
|
||||
for (i=0;i<d.n;i++) {
|
||||
d.x[i]=x[i];
|
||||
d.y[i]=y[i];
|
||||
d.sy[i]=sy[i];
|
||||
d.w[i]=1.0/(d.sy[i]*d.sy[i]);
|
||||
}
|
||||
|
||||
gsl_vector_view a=gsl_vector_view_array(q,m);
|
||||
gsl_vector_view wts=gsl_vector_view_array(d.w,d.n);
|
||||
gsl_matrix *covar=gsl_matrix_alloc(m,m);
|
||||
|
||||
// Function definition
|
||||
fdf.f=gauss_f;
|
||||
fdf.df=gauss_df;
|
||||
fdf.fvv=NULL;
|
||||
fdf.n=d.n;
|
||||
fdf.p=m;
|
||||
fdf.params=&d;
|
||||
|
||||
// Allocate and initialize
|
||||
w=gsl_multifit_nlinear_alloc(T,&fdf_params,d.n,m);
|
||||
gsl_multifit_nlinear_winit(&a.vector,&wts.vector,&fdf,w);
|
||||
f=gsl_multifit_nlinear_residual(w);
|
||||
|
||||
// Solve and compute covariance and chi-squared
|
||||
status=gsl_multifit_nlinear_driver(20,xtol,gtol,ftol,NULL,NULL,&info,w);
|
||||
J=gsl_multifit_nlinear_jac(w);
|
||||
gsl_multifit_nlinear_covar(J,0.0,covar);
|
||||
gsl_blas_ddot(f,f,chisq);
|
||||
|
||||
// Extract parameters
|
||||
for (i=0;i<m;i++) {
|
||||
q[i]=gsl_vector_get(w->x,i);
|
||||
sq[i]=sqrt(gsl_matrix_get(covar,i,i));
|
||||
}
|
||||
|
||||
// Free
|
||||
gsl_multifit_nlinear_free (w);
|
||||
gsl_matrix_free (covar);
|
||||
gsl_rng_free (r);
|
||||
free(d.x);
|
||||
free(d.y);
|
||||
free(d.sy);
|
||||
free(d.w);
|
||||
|
||||
return status;
|
||||
}
|
||||
|
||||
4
makefile
4
makefile
|
|
@ -25,8 +25,8 @@ rffind: rffind.o rfio.o rftime.o
|
|||
rfedit: rfedit.o rfio.o rftime.o
|
||||
$(CC) -o rfedit rfedit.o rfio.o rftime.o -lm
|
||||
|
||||
rfplot: rfplot.o rftime.o rfio.o rftrace.o sgdp4.o satutl.o deep.o ferror.o
|
||||
gfortran -o rfplot rfplot.o rftime.o rfio.o rftrace.o sgdp4.o satutl.o deep.o ferror.o $(LFLAGS)
|
||||
rfplot: rfplot.o rftime.o rfio.o rftrace.o sgdp4.o satutl.o deep.o ferror.o gaussfit.o
|
||||
gfortran -o rfplot rfplot.o rftime.o rfio.o rftrace.o sgdp4.o satutl.o deep.o ferror.o gaussfit.o $(LFLAGS)
|
||||
|
||||
rfspec: rfspec.o rftime.o rfio.o
|
||||
gfortran -o rfspec rfspec.o rftime.o rfio.o $(LFLAGS)
|
||||
|
|
|
|||
3
rffit.c
3
rffit.c
|
|
@ -1,6 +1,7 @@
|
|||
#include <stdio.h>
|
||||
#include <string.h>
|
||||
#include <stdlib.h>
|
||||
#include <ctype.h>
|
||||
#include <math.h>
|
||||
#include <cpgplot.h>
|
||||
#include <getopt.h>
|
||||
|
|
@ -83,7 +84,7 @@ struct site get_site(int site_id)
|
|||
file=fopen(filename,"r");
|
||||
if (file==NULL) {
|
||||
printf("File with site information not found!\n");
|
||||
return;
|
||||
exit(0);
|
||||
}
|
||||
while (fgets(line,LIM,file)!=NULL) {
|
||||
// Skip
|
||||
|
|
|
|||
124
rfplot.c
124
rfplot.c
|
|
@ -2,6 +2,7 @@
|
|||
#include <string.h>
|
||||
#include <stdlib.h>
|
||||
#include <math.h>
|
||||
#include <ctype.h>
|
||||
#include <cpgplot.h>
|
||||
#include <getopt.h>
|
||||
#include "rftime.h"
|
||||
|
|
@ -22,9 +23,11 @@ void time_axis(double *mjd,int n,float xmin,float xmax,float ymin,float ymax);
|
|||
void usage(void);
|
||||
void plot_traces(struct trace *t,int nsat,float fcen);
|
||||
struct trace fit_trace(struct spectrogram s,struct select sel,int site_id,int graves);
|
||||
struct trace fit_gaussian_trace(struct spectrogram s,struct select sel,int site_id,int graves);
|
||||
void convolve(float *y,int n,float *w,int m,float *z);
|
||||
float gauss(float x,float w);
|
||||
void quadfit(float x[],float y[],int n,float a[]);
|
||||
int fit_gaussian(double *x,double *y,double *sy,int n,double *q,double *sq,int m,double *chisq);
|
||||
|
||||
// Fit trace
|
||||
struct trace locate_trace(struct spectrogram s,struct trace t,int site_id)
|
||||
|
|
@ -629,6 +632,13 @@ int main(int argc,char *argv[])
|
|||
continue;
|
||||
}
|
||||
|
||||
// Fit
|
||||
if (c=='F') {
|
||||
tf=fit_gaussian_trace(s,sel,site_id,graves);
|
||||
tf.site=site_id;
|
||||
continue;
|
||||
}
|
||||
|
||||
// Identify
|
||||
if (c=='i') {
|
||||
if (graves==0)
|
||||
|
|
@ -1243,6 +1253,120 @@ struct trace fit_trace(struct spectrogram s,struct select sel,int site_id,int gr
|
|||
return t;
|
||||
}
|
||||
|
||||
// Fit trace
|
||||
struct trace fit_gaussian_trace(struct spectrogram s,struct select sel,int site_id,int graves)
|
||||
{
|
||||
int i,j,k,l,m,sn;
|
||||
int i0,i1,j0,j1,jmax,nw,status;
|
||||
float x,y,s1,s2,z,za,zs,zm,sigma;
|
||||
double f,*px,*py,*psy,a[5],sa[5],chisq;
|
||||
struct trace t;
|
||||
FILE *file;
|
||||
|
||||
// Allocate
|
||||
nw=2*(int) sel.w;
|
||||
px=(double *) malloc(sizeof(double)*nw);
|
||||
py=(double *) malloc(sizeof(double)*nw);
|
||||
psy=(double *) malloc(sizeof(double)*nw);
|
||||
|
||||
// Set up trace
|
||||
t.satno=99999;
|
||||
t.n=(int) ceil(sel.x[sel.n-1]-sel.x[0]);
|
||||
t.mjd=(double *) malloc(sizeof(double)*t.n);
|
||||
t.freq=(double *) malloc(sizeof(double)*t.n);
|
||||
t.za=(float *) malloc(sizeof(float)*t.n);
|
||||
|
||||
// Open file
|
||||
file=fopen("out.dat","w");
|
||||
|
||||
// Loop over selected regions
|
||||
for (k=0,l=0;k<sel.n-1;k++) {
|
||||
for (x=sel.x[k];x<=sel.x[k+1];x+=1.0) {
|
||||
y=(x-sel.x[k])/(sel.x[k+1]-sel.x[k])*(sel.y[k+1]-sel.y[k])+sel.y[k];
|
||||
i=(int) floor(x);
|
||||
j0=(int) floor(y-sel.w);
|
||||
j1=(int) floor(y+sel.w);
|
||||
|
||||
// Keep in range
|
||||
if (j0<0)
|
||||
j0=0;
|
||||
if (j1>=s.nchan)
|
||||
j1=s.nchan;
|
||||
|
||||
// Fill arrays
|
||||
for (j=j0,m=0;j<j1;j++,m++) {
|
||||
px[m]=(float) j;
|
||||
py[m]=s.z[i+s.nsub*j];
|
||||
psy[m]=1.0;
|
||||
}
|
||||
|
||||
// Find maximum and significance
|
||||
zm=0.0;
|
||||
jmax=0;
|
||||
s1=0.0;
|
||||
s2=0.0;
|
||||
sn=0;
|
||||
for (j=j0;j<j1;j++) {
|
||||
z=s.z[i+s.nsub*j];
|
||||
s1+=z;
|
||||
s2+=z*z;
|
||||
sn++;
|
||||
if (z>zm) {
|
||||
zm=z;
|
||||
jmax=j;
|
||||
}
|
||||
}
|
||||
|
||||
// Estimate parameters
|
||||
a[0]=jmax;
|
||||
a[1]=1.0;
|
||||
a[2]=zm;
|
||||
a[3]=0.0;
|
||||
a[4]=0.0;
|
||||
|
||||
// Fit gaussian
|
||||
status=fit_gaussian(px,py,psy,m,a,sa,5,&chisq);
|
||||
|
||||
|
||||
printf("%d %f %f %f %f %f %f %d %f\n",status,chisq,a[0],a[1],a[2],a[3],a[4],jmax,zm);
|
||||
|
||||
// Remove
|
||||
s1-=zm;
|
||||
s2-=zm*zm;
|
||||
sn--;
|
||||
za=s1/(float) sn;
|
||||
zs=sqrt(s2/(float) sn-za*za);
|
||||
sigma=(zm-za)/zs;
|
||||
|
||||
// Store
|
||||
if (sigma>sel.sigma && s.mjd[i]>1.0) {
|
||||
f=s.freq-0.5*s.samp_rate+(double) jmax*s.samp_rate/(double) s.nchan;
|
||||
if (graves==0)
|
||||
fprintf(file,"%lf %lf %f %d\n",s.mjd[i],f,sigma,site_id);
|
||||
else
|
||||
fprintf(file,"%lf %lf %f %d 9999\n",s.mjd[i],f,sigma,site_id);
|
||||
cpgpt1((float) i,(float) jmax,17);
|
||||
cpgpt1((float) i,(float) a[0],24);
|
||||
t.mjd[l]=s.mjd[i];
|
||||
t.freq[l]=f;
|
||||
t.za[l]=0.0;
|
||||
l++;
|
||||
}
|
||||
}
|
||||
}
|
||||
t.n=l;
|
||||
|
||||
// Close file
|
||||
fclose(file);
|
||||
|
||||
// Free
|
||||
free(px);
|
||||
free(py);
|
||||
free(psy);
|
||||
|
||||
return t;
|
||||
}
|
||||
|
||||
// Discrete convolution
|
||||
void convolve(float *y,int n,float *w,int m,float *z)
|
||||
{
|
||||
|
|
|
|||
Loading…
Reference in New Issue