From 6a99239b6f30f5ea9e8e2c179f256346e516c08c Mon Sep 17 00:00:00 2001 From: Cees Bassa Date: Wed, 3 Dec 2014 21:17:34 +0100 Subject: [PATCH] Added rffit --- dsmin.c | 124 ++++ frequencies.txt | 10 + makefile | 5 +- rffit.c | 1711 +++++++++++++++++++++++++++++++++++++++++++++++ simplex.c | 28 + versafit.c | 168 +++++ 6 files changed, 2045 insertions(+), 1 deletion(-) create mode 100644 dsmin.c create mode 100644 rffit.c create mode 100644 simplex.c create mode 100644 versafit.c diff --git a/dsmin.c b/dsmin.c new file mode 100644 index 0000000..7222508 --- /dev/null +++ b/dsmin.c @@ -0,0 +1,124 @@ +#include +#include +#include + +#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=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 +#include +#include +#include +#include +#include +#include "sgdp4h.h" + +#define LIM 80 +#define NMAX 1024 +#define D2R M_PI/180.0 +#define R2D 180.0/M_PI +#define XKMPER 6378.137 // km +#define KG 0.07436680 // earth radii, earth masses, minutes +#define C 299792.458 // km/s +#define FLAT (1.0/298.257) + +struct site { + int id; + double lat,lng; + float alt; + char observer[64]; +} site; +struct point { + char timestamp[24]; + double mjd,freq,v; + float t,f; + float flux; + int flag,site_id,rsite_id; + struct site s,r; +}; +struct data { + int n; + struct point *p; + int fitfreq; + double mjdmin,mjdmax,mjd0; + float freqmin,freqmax,fluxmin,fluxmax,f0,ffit; + char satname[LIM]; +} d; +orbit_t orb; +int fgetline(FILE *file,char *s,int lim); +struct data read_data(char *filename); +double date2mjd(int year,int month,double day); +void mjd2nfd(double mjd,char *nfd); +struct point decode_line(char *line); +double modulo(double,double); +double gmst(double); +double dgmst(double); +void obspos_xyz(double,struct site,xyz_t *,xyz_t *); +int velocity(orbit_t orb,double mjd,struct site s,double *v,double *azi,double *alt); +void deselect_inside(float x0,float y0,float x,float y); +void highlight(float x0,float y0,float x,float y,int flag); +void deselect_outside(float xmin,float ymin,float xmax,float ymax); +void deselect_nearest(float x,float y,float xmin,float ymin,float xmax,float ymax); +void save_data(float xmin,float ymin,float xmax,float ymax,char *filename); +void equatorial2horizontal(double mjd,struct site s,double ra,double de,double *azi,double *alt); +double chisq(double a[]); +void versafit(int m,int n,double *a,double *da,double (*func)(double *),double dchisq,double tol,char *opt); +double compute_rms(void); +void mjd2date(double mjd,int *year,int *month,double *day); +void print_tle(orbit_t orb,char *filename); +void search(void); +double fit_curve(orbit_t orb,int *ia); +double mjd2doy(double mjd,int *yr); +double doy2mjd(int year,double doy); + +// Get observing site +struct site get_site(int site_id) +{ + int i=0,status; + char line[LIM]; + FILE *file; + int id; + double lat,lng; + float alt; + char abbrev[3],observer[64]; + struct site s; + char *env,filename[LIM]; + + env=getenv("ST_DATADIR"); + sprintf(filename,"%s/data/sites.txt",env); + + file=fopen(filename,"r"); + if (file==NULL) { + printf("File with site information not found!\n"); + return; + } + while (fgets(line,LIM,file)!=NULL) { + // Skip + if (strstr(line,"#")!=NULL) + continue; + + // Strip newline + line[strlen(line)-1]='\0'; + + // Read data + status=sscanf(line,"%4d %2s %lf %lf %f", + &id,abbrev,&lat,&lng,&alt); + strcpy(observer,line+38); + + // Change to km + alt/=1000.0; + + // Copy site + if (id==site_id) { + s.lat=lat; + s.lng=lng; + s.alt=alt; + s.id=id; + strcpy(s.observer,observer); + } + + } + fclose(file); + + return s; +} + +// Select diagonal +void diagonal_select(float x0,float y0,float x1,float y1,int flag) +{ + int i; + float v; + float ymin,ymax; + printf("%f %f %f %f\n",x0,y0,x1,y1); + + for (i=0;i0.0 && v<1.0 && d.p[i].f>ymin && d.p[i].f1e-9) { + sprintf(bstar,"%11.4e",10*orb.bstar); + sbstar[0] = bstar[0]; sbstar[1] = bstar[1]; sbstar[2] = bstar[3]; sbstar[3] = bstar[4]; + sbstar[4] = bstar[5]; sbstar[5] = bstar[6]; sbstar[6] = bstar[8]; sbstar[7] = bstar[10]; sbstar[8] = '\0'; + } + // Print lines + sprintf(line1,"1 %05dU %2d%012.8f .00000000 00000-0 %8s 0 0",orb.satno,orb.ep_year-2000,orb.ep_day,sbstar); + sprintf(line2,"2 %05d %8.4f %8.4f %07.0f %8.4f %8.4f %11.8f 0",orb.satno,DEG(orb.eqinc),DEG(orb.ascn),1E7*orb.ecc,DEG(orb.argp),DEG(orb.mnan),orb.rev); + + // Compute checksums + for (i=0,csum=0;i0) { + printf("Identified %d candidate(s), best fitting satellite is %05d.\n",i,satno); + read_twoline(fp,satno,&orb); + rms=fit_curve(orb,ia); + } else { + printf("No candidates found.\n"); + satno=-1; + } + fclose(fp); + + return satno; +} + +void usage() +{ + printf("dpplot -d -c [tle catalog] -i [satno] -h\n\ndata file: Tabulated doppler curve\ntle catalog: Catalog with TLE's (optional)\nsatno: Satellite to load from TLE catalog (optional)\n\n"); + + return; +} + +int main(int argc,char *argv[]) +{ + int i,j,flag,redraw=1,plot_curve=1,plot_type=1,residuals=0,elset=0; + int imode,year,style,color; + char xlabel[64],ylabel[32],text[64]; + int site_id=4171; + float xmin,xmax,ymin,ymax; + float xminsel,xmaxsel,yminsel,ymaxsel; + float x0,y0,x,y; + double mjd,v,v1,azi,alt,rms=0.0,day,mjdtca=56658.0; + float t,f,vtca; + char c,nfd[32]="2014-01-01T00:00:00"; + int mode=0,posn=0,click=0; + char *catalog,*datafile,filename[64],string[64],bstar[10]=" 00000-0"; + int arg=0,nobs=0; + FILE *fp,*std; + char line0[72],line1[72],line2[72]; + int ia[]={0,0,0,0,0,0}; + float dx[]={0.1,0.1,0.4,0.4,0.7,0.7},dy[]={0.0,-0.25,0.0,-0.25,0.0,-0.25}; + int satno=-1,status; + struct site s0,s1; + int site_number[16],nsite=0; + + // Decode options + while ((arg=getopt(argc,argv,"d:c:i:hs:"))!=-1) { + switch(arg) { + case 'd': + datafile=optarg; + break; + + case 'c': + catalog=optarg; + break; + + case 'i': + satno=atoi(optarg); + break; + + case 'h': + usage(); + return 0; + break; + + case 's': + site_id=atoi(optarg); + break; + + default: + usage(); + return 0; + } + } + + // Read data + d=read_data(datafile); + d.fitfreq=1; + + // Count number of sites and assign colors + for (i=0;i=16) { + printf("Too many observing sites.\n"); + return 0; + } + } + + // Set default observing site + site=get_site(site_id); + + // Read TLE + if (satno>=0) { + fp=fopen(catalog,"rb"); + if (fp==NULL) + fatal_error("File open failed for reading %s\n",catalog); + status=read_twoline(fp,satno,&orb); + if (status==-1) { + printf("No elements found for %5d\n",satno); + satno=-1; + } + fclose(fp); + } + + freopen("/tmp/stderr.txt","w",stderr); + + cpgopen("/xs"); + cpgask(0); + + xmin=d.mjdmin-d.mjd0-0.005; + xmax=d.mjdmax-d.mjd0+0.005; + ymin=-12.0*d.f0/C; + ymax=12*d.f0/C; + // ymin=d.freqmin; + // ymax=d.freqmax; + + day=mjd2doy(d.mjd0,&year); + sprintf(xlabel,"MJD - %5.0f (%02d%03.0f)",d.mjd0,year-2000,floor(day)); + sprintf(ylabel,"Frequency - %.0f kHz",d.f0); + + // For ever loop + for (;;) { + if (redraw==1) { + // Plot buttons + cpgpage(); + cpgsvp(0.1,0.95,0.0,0.2); + cpgswin(0.0,1.0,-0.5,0.5); + + // Buttons + cpgtext(0.12,-0.05,"Inclination"); + + cpgtext(0.42,-0.05,"Eccentricity"); + cpgpt1(0.4,0.0,19); + cpgtext(0.72,-0.05,"Mean Anomaly"); + cpgpt1(0.7,0.0,19); + cpgtext(0.12,-0.3,"Ascending Node"); + cpgpt1(0.1,-0.25,19); + cpgtext(0.42,-0.3,"Arg. of Perigee"); + cpgpt1(0.4,-0.25,19); + cpgtext(0.72,-0.3,"Mean Motion"); + cpgpt1(0.7,-0.25,19); + + // Toggles + for (i=0;i<6;i++) { + cpgpt1(dx[i],dy[i],19); + if (ia[i]==1) { + cpgsci(2); + cpgpt1(dx[i],dy[i],16); + cpgsci(1); + } + } + + // Sky plot + cpgsvp(0.715,0.95,0.5,0.99); + cpgwnad(-90.0,90.0,-90.0,90.0); + cpgbox("BC",0.,0,"BC",0.,0); + cpgsfs(2); + cpgcirc(0.0,0.0,90.0); + cpgpt1(0.0,0.0,2); + + // Plot orbit + if (satno>0 && plot_curve==1) { + // Initialize + imode=init_sgdp4(&orb); + if (imode==SGDP4_ERROR) + break; + + cpgsci(15); + // for (mjd=d.mjd0,i=0;mjd0) { + if (vtca*v<0.0) { + mjdtca=mjd; + } + } + + x=(90-alt)*sin(azi*D2R); + y=-(90-alt)*cos(azi*D2R); + if (i==0 || alt<0.0) + cpgmove(x,y); + else + cpgdraw(x,y); + vtca=v; + } + mjd2nfd(mjdtca,nfd); + cpgsci(1); + cpgsls(1); + + // Plot points + for (i=0;i0 && plot_curve==1 && residuals==0) { + + // Plot tle + format_tle(orb,line1,line2); + cpgmtxt("T",2.0,0.0,0.0,line1); + cpgmtxt("T",1.0,0.0,0.0,line2); + + // Initialize + imode=init_sgdp4(&orb); + if (imode==SGDP4_ERROR) + break; + + // Loop over sites for plotting model + for (j=0;j0.0) { + cpgsls(1); + cpgsci(color); + } else { + cpgsls(2); + cpgsci(14); + } + + if (i==0) + cpgmove(t,f); + else + cpgdraw(t,f); + } + cpgsci(1); + cpgsls(1); + } + } + + // Plot selected points + for (i=0;i=1 && c-'0'<7) { + if (ia[c-49]==0) + ia[c-49]=1; + else if (ia[c-49]==1) + ia[c-49]=0; + redraw=1; + } + + // Change + if (c=='c') { + printf("(1) Inclination, (2) Ascending Node, (3) Eccentricity,\n(4) Arg. of Perigee, (5) Mean Anomaly, (6) Mean Motion,\n(7) B* drag, (8) Epoch, (9) Satellite ID\n\nWhich parameter to change: "); + status=scanf("%i",&i); + if (i>=0 && i<=9) { + printf("\nNew value: "); + fgets(string,64,stdin); + status=scanf("%s",string); + // if (i==0) strcpy(d.satname,string); + if (i==1) orb.eqinc=RAD(atof(string)); + if (i==2) orb.ascn=RAD(atof(string)); + if (i==3) orb.ecc=atof(string); + if (i==4) orb.argp=RAD(atof(string)); + if (i==5) orb.mnan=RAD(atof(string)); + if (i==6) orb.rev=atof(string); + if (i==7) orb.bstar=atof(string); + if (i==8) { + orb.ep_year=2000+(int) floor(atof(string)/1000.0); + orb.ep_day=atof(string)-1000*floor(atof(string)/1000.0); + } + if (i==9) orb.satno=atoi(string); + if (i==0) { + printf("%f\n",d.ffit); + d.ffit=atof(string); + d.fitfreq=0; + } + redraw=1; + } + printf("\n================================================================================\n"); + } + + // Move + if (c=='m') { + printf("Provide frequency offset (kHz): "); + status=scanf("%lf",&v); + + // Loop over points + for (i=0;i0) { + redraw=1; + plot_curve=1; + } + printf("\n================================================================================\n"); + } + + // Parameter search + // if (c=='S') { + // search(); + // redraw=1; + // } + + // Write TLE + if (c=='w') { + printf("TLE filename to write: "); + status=scanf("%s",filename); + print_tle(orb,filename); + printf("\n================================================================================\n"); + } + + // Reread tle + if (c=='R') { + // Read TLE + fp=fopen(catalog,"rb"); + if (fp==NULL) + fatal_error("File open failed for reading %s\n",catalog); + read_twoline(fp,satno,&orb); + print_orb(&orb); + printf("\n================================================================================\n"); + fclose(fp); + d.ffit=d.f0; + redraw=1; + } + + // Diagonal select + if (c=='k') { + printf("Selecting diagonal\n"); + diagonal_select(xmin,ymin,xmax,ymax,2); + redraw=1; + } + + // Toggle residuals + if (c=='j') { + if (residuals==0) + residuals=1; + else if (residuals==1) + residuals=0; + redraw=1; + } + + // Delete + if (c=='X') { + if (click==0) { + deselect_nearest(x,y,xmin,ymin,xmax,ymax); + + click=0; + redraw=1; + } + } + + // Zoom + if (c=='z') { + click=1; + mode=2; + } + // Delete box + if (c=='d') { + click=2; + mode=2; + } + + // Execute zoom, or box delete + if (c=='A') { + if (click==0) { + click=1; + } else if (click==1 && mode==2) { + xmin=(x0x) ? x0 : x; + ymin=(y0y) ? y0 : y; + + click=0; + mode=0; + redraw=1; + } else if (click==2 && mode==2) { + xminsel=(x0x) ? x0 : x; + yminsel=(y0y) ? y0 : y; + deselect_inside(xminsel,yminsel,xmaxsel,ymaxsel); + click=0; + mode=0; + redraw=1; + } else if (click==3 && mode==2) { + xminsel=(x0x) ? x0 : x; + yminsel=(y0y) ? y0 : y; + printf("%f %f %f %f\n",xminsel,xmaxsel,yminsel,ymaxsel); + click=0; + mode=0; + redraw=1; + } else { + click=0; + mode=0; + redraw=1; + } + } + + // Highlight + if (c=='s') { + highlight(xmin,ymin,xmax,ymax,2); + for (i=0,nobs=0;i3) + elset=0; + print_orb(&orb); + printf("\n================================================================================\n"); + click=0; + redraw=1; + continue; + } + /* + // Default tle + if (c=='t') { + orb.satno=99999; + orb.eqinc=0.5*M_PI; + orb.ascn=0.0; + orb.ecc=0.0; + orb.argp=0.0; + orb.mnan=0.0; + orb.rev=14.0; + orb.bstar=0.0; + orb.ep_day=mjd2doy(0.5*(d.mjdmin+d.mjdmax),&orb.ep_year); + satno=99999; + print_orb(&orb); + printf("\n================================================================================\n"); + click=0; + redraw=1; + continue; + } + */ + // Unzoom + if (c=='r') { + xmin=d.mjdmin-d.mjd0; + xmax=d.mjdmax-d.mjd0; + ymin=-12.0*d.f0/C; + ymax=12.0*d.f0/C; + mode=0; + click=0; + redraw=1; + continue; + } + + + + // Help + if (c=='h') { + printf("Usage:\n================================================================================\nq Quit\np Toggle curve plotting\n1 Toggle fitting parameter (Inclination)\n2 Toggle fitting parameter (RA of ascending node)\n3 Toggle fitting parameter (Eccentricity)\n4 Toggle fitting parameter (Argyment of perigee)\n5 Toggle fitting parameter (Mean anomaly)\n6 Toggle fitting parameter (Mean motion)\nc Change parameter\nm Move highlighted points in frequency\nl Select points on flux limit\nf Fit highlighted points\ng Get TLE from catalog\ni Identify satellite from catalog\nw Write present TLE\nR Reread TLE from catalog\nX Delete nearest point (right mouse button)\nz Start box to zoom\nd Start box to delete points\nA Zoom/delete points (left mouse button)\nh Highlight points in present window\nx Deselect all except highlighted\nI Invert selection\nD Delete highlighted points\ns Save highlighted points into file\nU Deselect all points\nu Deselect highlighted points\nt Load template tle\nr Reset zoom\nh This help\n================================================================================\n\n"); + } + + + // Save + x0=x; + y0=y; + } + cpgclos(); + + + // Free + free(d.p); + + fclose(stderr); + + return 0; +} + +// Read a line of maximum length int lim from file FILE into string s +int fgetline(FILE *file,char *s,int lim) +{ + int c,i=0; + + while (--lim > 0 && (c=fgetc(file)) != EOF && c != '\n') + s[i++] = c; + if (c == '\t') + c=' '; + if (c == '\n') + s[i++] = c; + s[i] = '\0'; + return i; +} + +// Decode line +struct point decode_line(char *line) +{ + int year,month,iday,hour,min,sec,fsec; + double day; + struct point p; + int site_id,rsite_id,nread,status; + + // Get timestamp + if (line[5]=='.') { + site_id=-1; + rsite_id=-1; + nread=sscanf(line,"%lf %lf %f %d %d",&p.mjd,&p.freq,&p.flux,&site_id,&rsite_id); + } else { + strncpy(p.timestamp,line,19); + p.timestamp[19]='\0'; + + // Get MJD, freq, flux + site_id=-1; + status=sscanf(line,"%2d/%2d/%2d %2d:%2d:%2d.%1d %lf %f %d",&year,&month,&iday,&hour,&min,&sec,&fsec,&p.freq,&p.flux,&site_id); + + day=iday+hour/24.0+min/1440.0+(sec+0.1*fsec)/86400.0; + p.mjd=date2mjd(2000+year,month,day); + } + + // Decode site + p.s=get_site(site_id); + p.site_id=site_id; + if (rsite_id!=-1) { + p.r=get_site(rsite_id); + p.rsite_id=rsite_id; + } else { + p.rsite_id=0; + } + + // Change to kHz + p.freq*=1E-3; + + // Set flag + p.flag=1; + + return p; +} + +// Read data +struct data read_data(char *filename) +{ + int i=0; + char line[LIM]; + FILE *file; + struct data d; + double c; + + // Open file + file=fopen(filename,"r"); + if (file==NULL) { + fprintf(stderr,"Failed to open %s\n",filename); + exit(1); + } + + // Count lines + while (fgetline(file,line,LIM)>0) + i++; + d.n=i; + + // Allocate + d.p=(struct point *) malloc(sizeof(struct point)*d.n); + + // Rewind file + rewind(file); + + // Read data + i=0; + while (fgetline(file,line,LIM)>0) + d.p[i++]=decode_line(line); + + // Close file + fclose(file); + + d.mjdmin=d.mjdmax=d.p[0].mjd; + d.freqmin=d.freqmax=d.p[0].freq; + d.fluxmin=d.fluxmax=d.p[0].flux; + for (i=1;i0) { + if (d.p[i].mjdd.mjdmax) d.mjdmax=d.p[i].mjd; + if (d.p[i].freqd.freqmax) d.freqmax=d.p[i].freq; + if (d.p[i].fluxd.fluxmax) d.fluxmax=d.p[i].flux; + } + } + c=0.1*(d.mjdmax-d.mjdmin); + d.mjdmin-=c; + d.mjdmax+=c; + c=0.1*(d.freqmax-d.freqmin); + d.freqmin-=c; + d.freqmax+=c; + c=0.02*(d.fluxmax-d.fluxmin); + d.fluxmin-=c; + d.fluxmax+=c; + + // Center time and frequency + d.mjd0=floor(0.5*(d.mjdmax+d.mjdmin)); + d.f0=floor(0.5*(d.freqmax+d.freqmin)); + + // Compute times + for (i=0;ix=gc*cos(s.lat*D2R)*cos(theta*D2R)*XKMPER; + pos->y=gc*cos(s.lat*D2R)*sin(theta*D2R)*XKMPER; + pos->z=gs*sin(s.lat*D2R)*XKMPER; + vel->x=-gc*cos(s.lat*D2R)*sin(theta*D2R)*XKMPER*dtheta; + vel->y=gc*cos(s.lat*D2R)*cos(theta*D2R)*XKMPER*dtheta; + vel->z=0.0; + + return; +} + +// Greenwich Mean Sidereal Time +double gmst(double mjd) +{ + double t,gmst; + + t=(mjd-51544.5)/36525.0; + + gmst=modulo(280.46061837+360.98564736629*(mjd-51544.5)+t*t*(0.000387933-t/38710000),360.0); + + return gmst; +} + +// Greenwich Mean Sidereal Time +double dgmst(double mjd) +{ + double t,dgmst; + + t=(mjd-51544.5)/36525.0; + + dgmst=360.98564736629+t*(0.000387933-t/38710000); + + return dgmst; +} + +// Return x modulo y [0,y) +double modulo(double x,double y) +{ + x=fmod(x,y); + if (x<0.0) x+=y; + + return x; +} + +// Deselect inside box +void deselect_inside(float x0,float y0,float x,float y) +{ + int i; + float xmin,xmax,ymin,ymax; + + xmin=(x0x) ? x0 : x; + ymin=(y0y) ? y0 : y; + for (i=0;ixmin && d.p[i].tymin && d.p[i].fx) ? x0 : x; + ymin=(y0y) ? y0 : y; + for (i=0;ixmin && d.p[i].tymin && d.p[i].fxmax || d.p[i].fymax) + d.p[i].flag=0; + + return; +} + +// Deselect nearest point +void deselect_nearest(float x,float y,float xmin,float ymin,float xmax,float ymax) +{ + int i,imin,flag; + float r,rmin; + float dx,dy; + + for (i=0,flag=0;i0) { + dx=(x-d.p[i].t)/(xmax-xmin); + dy=(y-d.p[i].f)/(ymax-ymin);; + r=sqrt(dx*dx+dy*dy); + if (flag==0) { + imin=i; + rmin=r; + flag=1; + } + if (rxmin && d.p[i].tymin && d.p[i].f=1.0) + a[2]=0.999; + if (a[5]<0.05) + a[5]=0.05; + + // Set parameters + orb.eqinc=RAD(a[0]); + orb.ascn=RAD(modulo(a[1],360.0)); + orb.ecc=a[2]; + orb.argp=RAD(modulo(a[3],360.0)); + orb.mnan=RAD(modulo(a[4],360.0)); + orb.rev=a[5]; + + // Initialize + imode=init_sgdp4(&orb); + if (imode==SGDP4_ERROR) + printf("Error\n"); + + // Loop over highlighted points + for (i=0,sum1=0.0,sum2=0.0;i0) + rms=sqrt(rms/(float) n); + + return rms; +} + +// Compute Date from Julian Day +void mjd2nfd(double mjd,char *nfd) +{ + double f,jd,dday; + int z,alpha,a,b,c,d,e; + int year,month,day,hour,min; + float sec,x; + + jd=mjd+2400000.5; + jd+=0.5; + + z=floor(jd); + f=fmod(jd,1.); + + if (z<2299161) + a=z; + else { + alpha=floor((z-1867216.25)/36524.25); + a=z+1+alpha-floor(alpha/4.); + } + b=a+1524; + c=floor((b-122.1)/365.25); + d=floor(365.25*c); + e=floor((b-d)/30.6001); + + dday=b-d-floor(30.6001*e)+f; + if (e<14) + month=e-1; + else + month=e-13; + + if (month>2) + year=c-4716; + else + year=c-4715; + + day=(int) floor(dday); + x=24.0*(dday-day); + x=3600.*fabs(x); + sec=fmod(x,60.); + x=(x-sec)/60.; + min=fmod(x,60.); + x=(x-min)/60.; + hour=x; + sec=floor(1000.0*sec)/1000.0; + + sprintf(nfd,"%04d-%02d-%02dT%02d:%02d:%02.0f",year,month,day,hour,min,sec); + + return; +} + +// Compute Date from Julian Day +void mjd2date(double mjd,int *year,int *month,double *day) +{ + double f,jd; + int z,alpha,a,b,c,d,e; + + jd=mjd+2400000.5; + jd+=0.5; + + z=floor(jd); + f=fmod(jd,1.); + + if (z<2299161) + a=z; + else { + alpha=floor((z-1867216.25)/36524.25); + a=z+1+alpha-floor(alpha/4.); + } + b=a+1524; + c=floor((b-122.1)/365.25); + d=floor(365.25*c); + e=floor((b-d)/30.6001); + + *day=b-d-floor(30.6001*e)+f; + if (e<14) + *month=e-1; + else + *month=e-13; + + if (*month>2) + *year=c-4716; + else + *year=c-4715; + + return; +} + +// Print TLE +void print_tle(orbit_t orb,char *filename) +{ + int i,n; + FILE *file; + double mjdmin,mjdmax; + int year,month; + double day; + char line1[70],line2[70]; + + // Count number of points + for (i=0,n=0;imjdmax) mjdmax=d.p[i].mjd; + n++; + } + } + + // Write TLE + file=fopen(filename,"w"); + format_tle(orb,line1,line2); + fprintf(file,"%s\n%s\n",line1,line2); + + mjd2date(mjdmin,&year,&month,&day); + fprintf(file,"# %4d%02d%05.2lf-",year,month,day); + mjd2date(mjdmax,&year,&month,&day); + fprintf(file,"%4d%02d%05.2lf, %d measurements, %.3lf kHz rms\n",year,month,day,n,compute_rms()); + fclose(file); + + return; +} + +// Parameter search +void search(void) +{ + int i,j,n; + double a[6],da[6]; + FILE *file; + double xmin,xmax,ymin,ymax; + int nx,ny,status; + + // Get input + printf("Provide mean motion estimate: "); + status=scanf("%lf",&a[5]); + printf("Provide inclination estimate: "); + status=scanf("%lf",&a[0]); + printf("Provide mean anomaly range, steps [min max nstep]: "); + status=scanf("%lf %lf %d",&xmin,&xmax,&nx); + printf("Provide ascending node range, steps [min max nstep]: "); + status=scanf("%lf %lf %d",&ymin,&ymax,&ny); + + // Count highlighted points + for (i=0,n=0;i0) + versafit(n,6,a,da,chisq,0.0,1e-5,"n"); + + // Return parameters + orb.eqinc=RAD(a[0]); + orb.ascn=RAD(modulo(a[1],360.0)); + orb.ecc=a[2]; + orb.argp=RAD(modulo(a[3],360.0)); + orb.mnan=RAD(modulo(a[4],360.0)); + orb.rev=a[5]; + + // print_tle(orb); + rms=compute_rms(); + + return rms; +} + +// MJD to DOY +double mjd2doy(double mjd,int *yr) +{ + int year,month,k=2; + double day,doy; + + mjd2date(mjd,&year,&month,&day); + + if (year%4==0 && year%400!=0) + k=1; + + doy=floor(275.0*month/9.0)-k*floor((month+9.0)/12.0)+day-30; + + *yr=year; + + return doy; +} + +// DOY to MJD +double doy2mjd(int year,double doy) +{ + int month,k=2; + double day; + + if (year%4==0 && year%400!=0) + k=1; + + month=floor(9.0*(k+doy)/275.0+0.98); + + if (doy<32) + month=1; + + day=doy-floor(275.0*month/9.0)+k*floor((month+9.0)/12.0)+30.0; + + return date2mjd(year,month,day); +} diff --git a/simplex.c b/simplex.c new file mode 100644 index 0000000..68ac169 --- /dev/null +++ b/simplex.c @@ -0,0 +1,28 @@ +// Creates Simplex +#include +#include + +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;jj) p[i][j]=a[j]-da[j]; + } + } + + return p; +} + diff --git a/versafit.c b/versafit.c new file mode 100644 index 0000000..8b49a44 --- /dev/null +++ b/versafit.c @@ -0,0 +1,168 @@ +// Versatile Fitting Routine +#include +#include +#include +#include + +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