Added rffit

dsmin.c

@@ -0,0 +1,124 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+
+#define TINY 1.0e-10
+#define NMAX 100000
+#define SWAP(a,b) {(a)+=(b);(b)=(a)-(b);(a)-=(b);}
+
+// Downhill Simplex Minimization
+int dsmin(double **p,double *y,int n,double ftol,double (*func)(double *))
+{
+  int i,j,nfunk=0;
+  int ihi,ilo,ise;
+  double *ptry,*pmid,*psum;
+  double tol,ytry,rtol,ysave;
+  double *vector_sum(double **,int);
+  double dsmod(double **,double *,double *,int,double (*func)(double *),int,double);
+
+  // Allocate memory
+  psum=(double *) malloc(sizeof(double) * n);
+
+  // Get function values
+  for (i=0;i<=n;i++) 
+    y[i]=func(p[i]);
+
+  // Sum vectors
+  psum=vector_sum(p,n);
+
+  // Start forever loop
+  for (;;) {
+    // Find high and low point
+    ilo=0;
+    ihi = (y[0]>y[1]) ? (ise=1,0) : (ise=0,1);
+    for (i=0;i<=n;i++) {
+      if (y[i]<=y[ilo]) ilo=i;
+      if (y[i]>y[ihi]) {
+      ise=ihi;
+      ihi=i;
+      } else if (y[i]>y[ise] && i!=ihi) ise=i;
+    }
+
+    // Compute fractional range from highest to lowest point
+    rtol=2.0*fabs(y[ihi]-y[ilo])/(fabs(y[ihi])+fabs(y[ilo])+TINY);
+
+    // Return if fractional tolerance is acceptable
+    if (rtol<ftol) 
+      break;
+
+    if (nfunk>=NMAX) {
+      printf("dsmin: NMAX exceeded!\n");
+      return -1;
+    }
+    nfunk+=2;
+
+    // Reflect simplex
+    ytry=dsmod(p,y,psum,n,func,ihi,-1.0);
+
+    if (ytry<=y[ilo]) // Goes right direction, extrapolate by factor 2
+      ytry=dsmod(p,y,psum,n,func,ihi,2.0);
+    else if (ytry>=y[ise]) { // 1D contraction
+      ysave=y[ihi];
+      ytry=dsmod(p,y,psum,n,func,ihi,0.5);
+      if (ytry>=ysave) {
+	for (i=0;i<=n;i++) {
+	  if (i!=ilo) {
+	    for (j=0;j<n;j++) 
+	      p[i][j]=psum[j]=0.5*(p[i][j]+p[ilo][j]);
+	    y[i]=(*func)(psum);
+	  }
+	}
+	nfunk+=n;
+	
+        psum=vector_sum(p,n);
+      }
+    } else --nfunk;
+  }
+  free(psum);
+
+  return nfunk;
+}
+
+// Sum vectors
+double *vector_sum(double **p,int n)
+{
+  int i,j;
+  double sum,*psum;
+
+  psum=(double *) malloc(sizeof(double) * n);
+
+  for (i=0;i<n;i++) {
+    sum=0.;
+    for (j=0;j<=n;j++)
+      sum+=p[j][i];
+    psum[i]=sum;
+  }
+
+  return psum;
+}
+
+// Simplex modification
+double dsmod(double **p,double *y,double *psum,int n,double (*func)(double *),int ihi,double fac)
+{
+  int i;
+  double fac1,fac2,ytry,*ptry;
+
+  ptry=(double *) malloc(sizeof(double) * n);
+  fac1=(1.0-fac)/(double) n;
+  fac2=fac1-fac;
+
+  for (i=0;i<n;i++)
+    ptry[i]=psum[i]*fac1-p[ihi][i]*fac2;
+  ytry=(*func)(ptry);
+  if (ytry<y[ihi]) {
+    y[ihi]=ytry;
+    for (i=0;i<n;i++) {
+      psum[i] += ptry[i]-p[ihi][i];
+      p[ihi][i]=ptry[i];
+    }
+  }
+
+  free(ptry);
+  
+  return ytry;
+}

frequencies.txt

@@ -280,3 +280,13 @@
 23533 2243.290
 23533 2244.314
 23533 2245.338
+23533 2243.653
+39630 2244.314
+40137 2242.404
+39086 2280.005
+36121 2282.160
+39765 2280.005
+26619 2282.407
+19822 2280.509
+37387 2280.008
+39012 2279.672

makefile

@@ -8,7 +8,10 @@ LFLAGS = -lcpgplot -lpgplot -lX11 -lpng -lm -lgsl -lgslcblas
 CC = gcc
 
 all:
-	make rfedit rfplot rffft rfpng rffind
+	make rfedit rfplot rffft rfpng rffind rffit
+
+rffit: rffit.o sgdp4.o satutl.o deep.o ferror.o dsmin.o simplex.o versafit.o
+	gfortran -o rffit rffit.o sgdp4.o satutl.o deep.o ferror.o dsmin.o simplex.o versafit.o $(LFLAGS)
 
 rfpng: rfpng.o rftime.o rfio.o rftrace.o sgdp4.o satutl.o deep.o ferror.o
 	gfortran -o rfpng rfpng.o rftime.o rfio.o rftrace.o sgdp4.o satutl.o deep.o ferror.o $(LFLAGS)

simplex.c

@@ -0,0 +1,28 @@
+// Creates Simplex
+#include <stdio.h>
+#include <stdlib.h>
+
+double **simplex(int n,double *a,double *da)
+{
+  int i,j;
+  double **p;
+
+  // Allocate pointers to rows
+  p=(double **) malloc(sizeof(double *) * (n+1));
+
+  // Allocate rows and set pointers
+  for (i=0;i<=n;i++)
+    p[i]=(double *) malloc(sizeof(double) * (n+1)*n);
+
+  // Fill simplex
+  for (i=0;i<=n;i++) {
+    for (j=0;j<n;j++) {
+      if (i<j) p[i][j]=a[j];
+      if (i==j) p[i][j]=a[j]+da[j];
+      if (i>j) p[i][j]=a[j]-da[j];
+    }
+  }  
+
+  return p;
+}
+

versafit.c

@@ -0,0 +1,168 @@
+// Versatile Fitting Routine
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include <string.h>
+
+int OUTPUT=1; // Print output on screen (1 = yes; 0 = no)
+int ERRCOMP=0; // Set reduced Chi-Squared to unity (1 = yes; 0 = no)
+
+int dsmin(double **,double *,int,double,double (*func)(double *));
+double **simplex(int,double *,double *);
+double parabolic_root(double,double,double,double);
+
+// Versafit fitting routine
+//
+// Inputs:
+//    m:      number of datapoints
+//    n:      number of parameters
+//    a:      parameters
+//    da:     expected spread in parameters
+//    func:   function to fit (Chi-squared function)
+//    dchisq  difference in Chi-squared
+//    tol:    tolerance
+//    opt:    options
+//            - n: no output
+void versafit(int m,int n,double *a,double *da,double (*func)(double *),double dchisq,double tol,char *opt)
+{
+  int i,j,k,l,nfunk,kmax=50;
+  double chisqmin;
+  double *b,*db;
+  double **p,*y;
+  double d[2],errcomp;
+
+  // Decode options
+  if (strchr(opt,'n')!=NULL) OUTPUT=0;
+  if (strchr(opt,'e')!=NULL) ERRCOMP=1;
+
+  // Intialize y
+  y=(double *) malloc(sizeof(double) * (n+1));
+
+  if (dchisq>=0.) {
+    // Compute simplex and minimize function
+    p=simplex(n,a,da);
+    nfunk=dsmin(p,y,n,tol,func);
+
+    // Average parameters
+    for (i=0;i<n;i++) {
+      a[i]=0.;
+      for (j=0;j<=n;j++) 
+	a[i]+=p[j][i];
+      a[i]/=(double) (n+1);
+    }
+
+    // Compute minimum
+    chisqmin=func(a);
+
+    // Compute error compensation
+    if (ERRCOMP) errcomp=sqrt(chisqmin/(double) (m-n));
+  }
+
+  // Basic Information
+  if (OUTPUT) {
+    printf("VersaFIT:\n");
+    if (m!=0)
+      printf("Number of datapoints: %i\n",m);
+    printf("Number of parameters: %i\n",n);
+    printf("Chi-squared: %14.5f\n",chisqmin);
+    if (m!=0)
+      printf("Reduced Chi-squared: %14.5f\n",chisqmin/(double) (m-n));
+    if (ERRCOMP) printf("Error compensation: %.4f\n",errcomp);
+    printf("Number of iterations: %i\n",nfunk);
+  
+    printf("\nParameters:\n");
+
+    // No error estimation
+    if (dchisq==0.) {
+      for (i=0;i<n;i++) 
+	printf(" a(%i): %12.5f\n",i+1,a[i]);
+    }
+  }
+
+  // With error estimation
+  if (dchisq!=0.) {
+    b=(double *) malloc(sizeof(double) * n);
+    db=(double *) malloc(sizeof(double) * n);
+
+    for (i=0;i<n;i++) {
+      if (da[i]!=0.) {
+	for (j=0;j<n;j++) {
+	  b[j]=a[j];
+	  db[j]=da[j];
+	}
+	d[0]=-da[i];
+	db[i]=0.;
+
+	for (k=0;k<kmax;k++) {
+	  b[i]=a[i]+d[0];
+
+	  // Minimize
+	  p=simplex(n,b,db);
+	  nfunk+=dsmin(p,y,n,tol,func);
+
+	  // Average parameters
+	  for (l=0;l<n;l++) {
+	    b[l]=0.;
+	    for (j=0;j<=n;j++) 
+	      b[l]+=p[j][l];
+	    b[l]/=(double) (n+1);
+	  }
+	  d[0]=parabolic_root(d[0],func(b),chisqmin,dchisq);
+
+	  if (fabs(chisqmin+dchisq-func(b))<tol) break;
+	}
+
+	d[1]=-d[0];
+	db[i]=0.;
+
+	for (k=0;k<kmax;k++) {
+	  b[i]=a[i]+d[1];
+
+	  // Minimize
+	  p=simplex(n,b,db);
+	  nfunk+=dsmin(p,y,n,tol,func);
+
+	  // Average parameters
+	  for (l=0;l<n;l++) {
+	    b[l]=0.;
+	    for (j=0;j<=n;j++) 
+	      b[l]+=p[j][l];
+	    b[l]/=(double) (n+1);
+	  }
+	  d[1]=parabolic_root(d[1],func(b),chisqmin,dchisq);
+
+	  if (fabs(chisqmin+dchisq-func(b))<tol) break;
+	}
+	da[i]=0.5*(fabs(d[0])+fabs(d[1]));
+	if (ERRCOMP) da[i]*=errcomp;
+      }
+    }
+
+    if (OUTPUT) 
+      for (i=0;i<n;i++) 
+	printf(" a(%i): %12.5f +- %9.5f\n",i+1,a[i],da[i]);
+  }
+  if (OUTPUT) printf("\nTotal number of iterations: %i\n",nfunk);
+
+  //  free(p);
+  //  free(y);
+  //  free(b);
+  //  free(db);
+
+  return;
+}
+
+// Compute root
+double parabolic_root(double x,double y,double y0,double dy)
+{
+  double a;
+
+  if (fabs(x)<1e-9) {
+    printf("Division by zero in function 'parabolic_root'\n");
+    x=1e-9;
+  }
+
+  a=(y-y0)/(x*x);
+  
+  return sqrt(fabs(dy/a))*x/fabs(x);
+}