From 83ca5c21e00cc921d2963e7ae844b034038701df Mon Sep 17 00:00:00 2001
From: Cees Bassa <cbassa@gmail.com>
Date: Tue, 16 Jan 2018 19:50:14 +0100
Subject: [PATCH] Misc changes

---
 gaussfit.c | 133 +++++++++++++++++++++++++++++++++++++++++++++++++++++
 makefile   |   4 +-
 rffit.c    |   3 +-
 rfplot.c   | 124 +++++++++++++++++++++++++++++++++++++++++++++++++
 4 files changed, 261 insertions(+), 3 deletions(-)
 create mode 100644 gaussfit.c

diff --git a/gaussfit.c b/gaussfit.c
new file mode 100644
index 0000000..6226da5
--- /dev/null
+++ b/gaussfit.c
@@ -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;
+}
+
diff --git a/makefile b/makefile
index 4a06324..393f1d6 100644
--- a/makefile
+++ b/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)
diff --git a/rffit.c b/rffit.c
index bfc7903..5c07fd2 100644
--- a/rffit.c
+++ b/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
diff --git a/rfplot.c b/rfplot.c
index bfb9144..94c3b7c 100644
--- a/rfplot.c
+++ b/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)
 {