1
0
Fork 0

Misc changes and new additions

pull/2/head
Cees Bassa 2013-10-01 18:53:57 +02:00
parent c90740ab77
commit 377242d2ff
15 changed files with 716 additions and 47 deletions

View File

@ -8,7 +8,7 @@
#include <gsl/gsl_multifit.h>
#include <getopt.h>
#define NMAX 2048
#define NMAX 8192
#define LIM 128
#define D2R M_PI/180.0
#define R2D 180.0/M_PI

View File

@ -58,7 +58,7 @@ int main(int argc,char *argv[])
struct catalog cat,ast;
float sx,sy,q;
char *env,starfile[128];
float r,rmin=1.0,rmax=10.0,mmin=1.0,mmax=8.0,mag=8.0;
float r,rmin=1.0,rmax=10.0,mmin=2.0,mmax=8.0,mag=8.0;
// Environment variables
env=getenv("ST_DATADIR");
@ -68,11 +68,11 @@ int main(int argc,char *argv[])
img=read_fits(argv[1],0);
// Hard coded
img.ra0=337.29;
img.de0=-5.10;
img.ra0=atof(argv[2]);
img.de0=atof(argv[3]);
sx=-10.0;
sy=10.0;
q=0.0;
q=atof(argv[4]);
img.x0=0.5*(float) img.naxis1;
img.y0=0.5*(float) img.naxis2;

View File

@ -14819,4 +14819,5 @@
39231 13037H
39232 13043A
99212 13706B
99068 13770A
99999 99000A

View File

@ -49,7 +49,8 @@
6242 JM 42.9453 -2.8284 623 Jon Mikel
6241 JM 42.9565 -2.8203 619 Jon Mikel
4160 BD 51.2793 5.4768 35 Bram Dorreman
0000 MB 44.6500 -63.6000 0 Michael Boschat
0000 JH 28.7564 -17.8889 2308 Jan Hattenbach
#0000 MB 44.6500 -63.6000 0 Michael Boschat
#0000 IS 54.2825 26.7466 0 Ivan Serhey
#0000 RL 38.9483 -104.5617 2073 Ron Lee
#0000 CW 59.3104 18.1003 0 Christoffer Wallstenius

View File

@ -32,7 +32,7 @@ struct image read_fits(char *filename,int pnum);
void forward(double ra0,double de0,double ra,double de,double *x,double *y)
{
int i;
char pcode[4]="TAN";
char pcode[4]="STG";
double phi,theta;
struct celprm cel;
struct prjprm prj;
@ -72,8 +72,8 @@ int main(int argc,char *argv[])
struct jpeg_image jpg,out;
double rx,ry,rx0,ry0;
double x,y,d;
double drx=-5.0,dry=5.0;
double ra0=297.695833,de0=8.868333;
double drx=-10.0,dry=10.0;
double ra0=303.91,de0=47.43;
// Read image
img=read_fits(argv[1],0);
@ -82,8 +82,8 @@ int main(int argc,char *argv[])
// Offset
forward(img.ra0,img.de0,ra0,de0,&rx0,&ry0);
out.nx=12000;
out.ny=8000;
out.nx=6000;
out.ny=4000;
out.nz=3;
out.z=(float *) malloc(sizeof(float)*out.nx*out.ny*out.nz);

View File

@ -10,8 +10,6 @@
struct image {
int nx,ny,nz;
float *z;
float xmin,xmax,ymin,ymax;
double tr[6];
double mjd;
char nfd[32],observer[32];
int cospar;
@ -23,18 +21,52 @@ double date2mjd(int year,int month,double day);
double nfd2mjd(char *date);
void mjd2nfd(double mjd,char *nfd);
struct image rebin(struct image raw,int nbin)
{
int i,j,k;
int ii,jj,kk;
struct image img;
img.nx=raw.nx/nbin;
img.ny=raw.ny/nbin;
img.nz=1;
img.z=(float *) malloc(sizeof(float)*img.nx*img.ny*img.nz);
for (i=0;i<img.nx;i++) {
for (j=0;j<img.ny;j++) {
k=i+img.nx*j;
img.z[k]=0.0;
for (ii=0;ii<nbin;ii++) {
for (jj=0;jj<nbin;jj++) {
kk=ii+nbin*i+raw.nx*(jj+nbin*j);
img.z[k]+=raw.z[kk];
}
}
}
}
img.mjd=raw.mjd;
img.cospar=raw.cospar;
img.exptime=raw.exptime;
strcpy(img.nfd,raw.nfd);
strcpy(img.observer,raw.observer);
return img;
}
int main(int argc,char *argv[])
{
int arg;
struct image img;
char infile[64],outfile[64],nfd[32]="2000-01-01T00:00:00";
struct image img,raw;
char infile[64],outfile[64]="",nfd[32]="2000-01-01T00:00:00";
double mjd=51544.0,delay=0.0,tz=0.0;
int cospar=0;
char observer[32]="Cees Bassa";
float exptime=10.06;
int flag=0,nbin=1;
// Decode options
while ((arg=getopt(argc,argv,"i:t:o:d:Z:c:T:O:"))!=-1) {
while ((arg=getopt(argc,argv,"i:t:o:d:Z:c:T:O:b:"))!=-1) {
switch(arg) {
case 'i':
@ -45,12 +77,17 @@ int main(int argc,char *argv[])
delay=(double) atof(optarg);
break;
case 'b':
nbin=atoi(optarg);
break;
case 'Z':
tz=(double) atof(optarg);
break;
case 'o':
strcpy(outfile,optarg);
flag=1;
break;
case 't':
@ -75,8 +112,14 @@ int main(int argc,char *argv[])
}
}
if (infile!=NULL)
img=read_jpg(infile);
if (infile!=NULL) {
if (nbin==1) {
img=read_jpg(infile);
} else {
raw=read_jpg(infile);
img=rebin(raw,nbin);
}
}
if (nfd!=NULL) {
// Compute time
@ -89,6 +132,9 @@ int main(int argc,char *argv[])
img.mjd=mjd;
}
if (flag==0)
sprintf(outfile,"%s.fits",img.nfd);
// Set properties
img.cospar=cospar;
img.exptime=exptime;
@ -97,7 +143,10 @@ int main(int argc,char *argv[])
if (outfile!=NULL)
write_fits(img,outfile);
// Free
free(img.z);
if (nbin!=1)
free(raw.z);
return 0;
}
@ -160,20 +209,6 @@ struct image read_jpg(char *filename)
}
}
// Limits
img.xmin=0.0;
img.xmax=(float) img.nx;
img.ymin=0.0;
img.ymax=(float) img.ny;
// Translation matrix
img.tr[1]=(img.xmax-img.xmin)/(float) img.nx;
img.tr[2]=0.0;
img.tr[0]=img.xmin-0.5*img.tr[1];
img.tr[4]=0.0;
img.tr[5]=(img.ymax-img.ymin)/(float) img.ny;
img.tr[3]=img.ymin-0.5*img.tr[5];
// Free allocated memory
free(row_pointer[0]);
free(raw_image);

View File

@ -11,12 +11,17 @@ CC = gcc
F77 = gfortran
all:
make addwcs angular calibrate dec2sex faketle fitsheader fitskey imgstat jpg2fits jpgstack measure pgm2fits plotfits pstrack rde2iod reduce residuals runsched satfit satid satmap satorbit sex2dec skymap tle2ole tleinfo uk2iod viewer wcsfit
make addwcs angular calibrate dec2sex faketle fitsheader fitskey imgstat jpg2fits jpgstack measure pgm2fits plotfits pstrack rde2iod reduce residuals runsched satfit satid satmap satorbit sex2dec skymap tle2ole tleinfo uk2iod viewer wcsfit deproject slewto waitfor
slewto: slewto.o
$(CC) -o slewto slewto.o -lm
waitfor: waitfor.o
$(CC) -o waitfor waitfor.o -lm
deproject: deproject.o
$(F77) -o deproject deproject.o $(LFLAGS)
jpgstack: jpgstack.o
$(CC) -o jpgstack jpgstack.o -ljpeg

View File

@ -90,6 +90,8 @@ void plot_objects(char *filename)
cpgsci(0);
if (strstr(catalog,"classfd")!=NULL)
cpgsci(4);
if (strstr(catalog,"inttles")!=NULL)
cpgsci(3);
cpgmove(x0,y0);
cpgdraw(x1,y1);

View File

@ -151,16 +151,21 @@ int main(int argc,char *argv[])
char *datafile,catalog[LIM];
char *env;
float dt;
int verbose=0;
env=getenv("ST_TLEDIR");
sprintf(catalog,"%s/classfd.tle",env);
// Decode options
while ((arg=getopt(argc,argv,"d:c:hs:"))!=-1) {
while ((arg=getopt(argc,argv,"d:c:hs:v"))!=-1) {
switch(arg) {
case 'd':
datafile=optarg;
break;
case 'v':
verbose=1;
break;
case 'c':
strcpy(catalog,optarg);
break;
@ -184,6 +189,12 @@ int main(int argc,char *argv[])
// Read data
d=read_data(datafile);
if (verbose==1) {
for (i=0;i<d.n;i++) {
printf("%14.8lf %10.4f %10.4f %10.4f %10.6f %10.6f\n",d.p[i].mjd,d.p[i].obspos.x,d.p[i].obspos.y,d.p[i].obspos.z,d.p[i].ra,d.p[i].de);
}
}
if (split==1) {
split_file(d,dt);
} else {

View File

@ -1111,8 +1111,8 @@ double chisq(double a[])
} else if (a[0]>180.0) {
a[0]=180.0;
}
if (a[5]>17.0)
a[5]=17.0;
if (a[5]>20.0)
a[5]=20.0;
if (a[5]<0.1)
a[5]=0.1;

View File

@ -335,6 +335,8 @@ int main(int argc,char *argv[])
env=getenv("ST_TLEDIR");
sprintf(filename,"%s/classfd.tle",env);
plot_satellites(filename,img,0,img.mjd,img.exptime,4);
sprintf(filename,"%s/inttles.tle",env);
plot_satellites(filename,img,0,img.mjd,img.exptime,3);
sprintf(filename,"%s/catalog.tle",env);
plot_satellites(filename,img,0,img.mjd,img.exptime,0);

View File

@ -7,7 +7,7 @@
#include "cpgplot.h"
#include "sgdp4h.h"
#define LIM 80
#define LIM 90
#define NMAX 1024
#define MMAX 28368
#define D2R M_PI/180.0
@ -29,8 +29,8 @@ struct map {
int length;
char orientation[LIM];
char nfd[LIM],tlefile[LIM],observer[32];
char datadir[LIM],tledir[LIM];
int site_id;
char datadir[LIM],tledir[LIM],notamfile[LIM],xyzfile[LIM];
int site_id,notamflag,xyzflag;
float w;
} m;
struct globe {
@ -77,6 +77,8 @@ void initialize_setup(void)
m.mjd=nfd2mjd(m.nfd);
m.w=1.2;
m.h0=gmst(m.mjd);
m.notamflag=0;
m.xyzflag=0;
// Default settings
strcpy(m.observer,"Unknown");
@ -293,7 +295,7 @@ void plot_track(void)
cpgpt1(x,y,17);
}
cpgmove(x,y);
} else {
} else if (s.r>XKMPER) {
if (sqrt(x*x+y*y)<XKMPER && z<0.0)
cpgmove(x,y);
else
@ -318,6 +320,115 @@ void plot_track(void)
return;
}
// plot satellite track
void plot_xyz(void)
{
int i=0,nstep=500;
orbit_t orb;
xyz_t pos,vel;
double jd,dt,h,mjd,mjd0;
FILE *fp=NULL;
float x,y,z,r,v;
long imode;
int isci;
float isch;
char norad[7],line[LIM];
struct sat s;
double rsun,rearth;
double dx,dy,dz,sra,sde;
xyz_t sunpos,satpos;
cpgqci(&isci);
cpgqch(&isch);
cpgsci(8);
fp=fopen(m.xyzfile,"rb");
if (fp==NULL) {
fatal_error("File open failed for reading \"%s\"",m.xyzfile);
}
h=gmst(m.mjd);
while (fgetline(fp,line,LIM)>0) {
// Get satellite position
sscanf(line,"%lf %lf %lf %lf",&mjd,&satpos.x,&satpos.y,&satpos.z);
printf("%lf %f %f %f\n",mjd,satpos.x,satpos.y,satpos.z);
// Get positions
sunpos_xyz(m.mjd,&sunpos,&sra,&sde);
// Sun position from satellite
dx=-satpos.x+sunpos.x;
dy=-satpos.y+sunpos.y;
dz=-satpos.z+sunpos.z;
// Distances
rsun=sqrt(dx*dx+dy*dy+dz*dz);
rearth=sqrt(satpos.x*satpos.x+satpos.y*satpos.y+satpos.z*satpos.z);
// Angles
s.psun=asin(696.0e3/rsun)*R2D;
s.pearth=asin(6378.135/rearth)*R2D;
s.p=acos((-dx*satpos.x-dy*satpos.y-dz*satpos.z)/(rsun*rearth))*R2D;
// s.p=acos(((sunpos.x+satpos.x)*satpos.x+(sunpos.y+satpos.y)*satpos.y+(sunpos.z+satpos.z)*satpos.z)/(rsun*rearth))*R2D;
s.p-=s.pearth;
// Celestial position
s.r=sqrt(satpos.x*satpos.x+satpos.y*satpos.y+satpos.z*satpos.z);
s.ra=atan2(satpos.y,satpos.x)*R2D;
s.de=asin(satpos.z/s.r)*R2D;
// Latitude and longitude
s.lng=s.ra-gmst(m.mjd);;
s.lat=s.de;
s.x=satpos.x;
s.y=satpos.y;
s.z=satpos.z;
x=satpos.x;
y=satpos.y;
z=satpos.z;
rotate(0,-90.0,&x,&y,&z);
rotate(1,90.0,&x,&y,&z);
rotate(1,m.l0+h,&x,&y,&z);
rotate(0,m.b0,&x,&y,&z);
// Visibility
if (s.p<-s.psun)
cpgsci(14);
else if (s.p>-s.psun && s.p<s.psun)
cpgsci(15);
else if (s.p>s.psun)
cpgsci(7);
// Plot
if (i==0) {
plot_footprint(s);
if (!(sqrt(x*x+y*y)<XKMPER && z<0.0)) {
sprintf(norad," xyz");
cpgsch(0.6);
cpgtext(x,y,norad);
cpgsch(isch);
cpgpt1(x,y,17);
}
cpgmove(x,y);
i++;
} else if (s.r>XKMPER) {
if (sqrt(x*x+y*y)<XKMPER && z<0.0)
cpgmove(x,y);
else
cpgdraw(x,y);
}
}
cpgsls(1);
cpgsci(isci);
cpgsch(isch);
return;
}
void read_globe(void)
{
int i,status;
@ -585,6 +696,54 @@ void plot_terminator(void)
return;
}
// 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;
}
void plot_notam(char *filename)
{
int i,flag=0;
float x,y,z;
float l,b;
char line[LIM];
FILE *file;
file=fopen(filename,"r");
while (fgetline(file,line,LIM)>0) {
sscanf(line,"%f %f",&b,&l);
z=cos(l*D2R)*cos(b*D2R)*XKMPER;
x=sin(l*D2R)*cos(b*D2R)*XKMPER;
y=sin(b*D2R)*XKMPER;
rotate(1,m.l0,&x,&y,&z);
rotate(0,m.b0,&x,&y,&z);
if (z>0.0) {
if (flag==0) {
cpgmove(x,y);
flag=1;
} else {
cpgdraw(x,y);
}
} else {
flag=0;
}
}
return;
}
void plot_map(void)
{
int redraw=1,status;
@ -653,8 +812,18 @@ void plot_map(void)
cpgsci(1);
cpgbox("BC",0.,0,"BC",0.,0);
// Plot notam
if (m.notamflag==1) {
cpgsci(3);
plot_notam(m.notamfile);
cpgsci(1);
}
// Plot track
plot_track();
if (m.xyzflag==0)
plot_track();
else
plot_xyz();
}
// Reset redraw
@ -689,7 +858,7 @@ void plot_map(void)
m.lng=atan2(y,x)*R2D;
m.lat=asin(z/XKMPER)*R2D;
printf("%8.3f %8.3f\n",m.lng,m.lat);
m.l0=m.lng;
m.b0=m.lat;
redraw=1;
@ -765,7 +934,7 @@ int main(int argc,char *argv[])
initialize_setup();
// Decode options
while ((arg=getopt(argc,argv,"t:c:i:s:l:h"))!=-1) {
while ((arg=getopt(argc,argv,"t:c:i:s:l:hN:p:"))!=-1) {
switch (arg) {
case 't':
@ -789,6 +958,16 @@ int main(int argc,char *argv[])
m.length=atoi(optarg);
break;
case 'N':
strcpy(m.notamfile,optarg);
m.notamflag=1;
break;
case 'p':
strcpy(m.xyzfile,optarg);
m.xyzflag=1;
break;
case 'h':
usage();
return 0;

View File

@ -23,7 +23,7 @@ long Isatsel=0;
extern double SGDP4_jd0;
struct map {
double alpha0,delta0,ra0,de0,azi0,alt0;
double alpha0,delta0,ra0,de0,azi0,alt0,q;
double fov,mjd,gmst,w,wl,wb;
float length;
float minmag,maxmag,minrad,maxrad;
@ -1876,6 +1876,19 @@ long identify_satellite(char *filename,int satno,double mjd,float rx,float ry)
return Isatmin;
}
double parallactic_angle(double mjd,double ra,double de)
{
double h,q;
h=gmst(mjd)+m.lng-ra;
q=atan2(sin(h*D2R),(tan(m.lat*D2R)*cos(de*D2R)-sin(de*D2R)*cos(h*D2R)))*R2D;
if (fabs(m.lat-de)<0.001)
q=0.0;
return q;
}
void skymap_plotsun(void)
{
double rx,ry;
@ -1895,6 +1908,8 @@ void skymap_plotsun(void)
return;
}
// plot skymap
int plot_skymap(void)
{
@ -1918,6 +1933,9 @@ int plot_skymap(void)
horizontal2equatorial(m.mjd,m.azi0,m.alt0,&m.ra0,&m.de0);
else if (strcmp(m.orientation,"equatorial")==0)
equatorial2horizontal(m.mjd,m.ra0,m.de0,&m.azi0,&m.alt0);
// Parallactic angle
m.q=parallactic_angle(m.mjd,m.ra0,m.de0);
// Get sun position
sunpos_xyz(m.mjd,&sunpos,&m.sra,&m.sde);
@ -2003,7 +2021,7 @@ int plot_skymap(void)
// Bottom string
dec2sex(m.ra0/15.0,sra,0,5);
dec2sex(m.de0,sde,0,4);
sprintf(text,"R: %s; D: %s; A: %.1f; E: %.1f; S: %.1fx%.1f deg; L: %d; O: %s; m < %.1f; f: %.0f mm; l: %.0f s",sra,sde,modulo(m.azi0-180.0,360.0),m.alt0,3.0*m.w,2.0*m.w,m.level,m.orientation,m.maxmag,focallength[fov],m.length);
sprintf(text,"R:%s; D:%s; A: %.1f; E: %.1f; q: %.2f; S: %.1fx%.1f deg; L: %d; O: %s; m < %.1f; f: %.0f mm; l: %.0f s",sra,sde,modulo(m.azi0-180.0,360.0),m.alt0,m.q,3.0*m.w,2.0*m.w,m.level,m.orientation,m.maxmag,focallength[fov],m.length);
cpgmtxt("B",1.0,0.0,0.0,text);
cpgsch(1.0);
@ -2419,7 +2437,7 @@ void dec2sex(double x,char *s,int f,int len)
int i;
double sec,deg,min;
char sign;
char *form[]={":: ",",, ","hms"," "};
char *form[]={"::",",,","hms"," "};
sign=(x<0 ? '-' : ' ');
x=3600.*fabs(x);

397
slewto.c 100644
View File

@ -0,0 +1,397 @@
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <getopt.h>
#include <time.h>
#include <netdb.h>
#include <sys/socket.h>
#include <netinet/in.h>
#include <netinet/in_systm.h>
#include <netinet/ip.h>
#include <netinet/tcp.h>
#include <arpa/inet.h>
#define LIM 128
#define D2R M_PI/180.0
#define R2D 180.0/M_PI
#define PORT 7264
#define IP "127.0.0.1"
struct map {
double alpha0,delta0,ra0,de0,azi0,alt0,q;
char orientation[LIM];
char nfd[LIM];
char datadir[LIM],observer[32];
double lat,lng;
double mjd;
float alt,timezone;
int site_id;
} m;
void usage(void)
{
return;
}
// Send new position to telescope
void send_position(char *sra,char *sde)
{
int skt;
struct hostent *he;
struct sockaddr_in addr;
char packet[2048];
FILE *file;
float ra,de;
// Old packet style
// sprintf(packet,"<newNumberVector device='Celestron GPS' name='EQUATORIAL_EOD_COORD_REQUEST'><oneNumber name='RA'>%s</oneNumber><oneNumber name='DEC'>%s</oneNumber></newNumberVector>",sra,sde);
// New packet style (as of 2013-08-20)
sprintf(packet,"<newNumberVector device='Celestron GPS' name='EQUATORIAL_EOD_COORD'><oneNumber name='RA'>%s</oneNumber><oneNumber name='DEC'>%s</oneNumber></newNumberVector>",sra,sde);
// Send TCP packet
skt=socket(AF_INET,SOCK_STREAM,0);
addr.sin_family=AF_INET;
addr.sin_port=htons(PORT);
he=gethostbyname(IP);
bcopy(he->h_addr,(struct in_addr *) &addr.sin_addr,he->h_length);
if(connect(skt,(struct sockaddr *) &addr,sizeof(addr))<0) {
fprintf(stderr,"Connection refused by remote host.\n");
return;
}
write(skt,packet,strlen(packet));
close(skt);
// Set restart
// file=fopen("/media/video/satobs/control/state.txt","w");
// if (file!=NULL) {
// fprintf(file,"restart");
// fclose(file);
// }
// Set position
// file=fopen("/media/video/satobs/control/position.txt","w");
// if (file!=NULL) {
// fprintf(file,"%s %s\n",sra,sde);
// fclose(file);
// }
return;
}
// Get observing site
void get_site(int site_id)
{
int i=0;
char line[LIM];
FILE *file;
int id;
double lat,lng;
float alt;
char abbrev[3],observer[64],filename[LIM];
sprintf(filename,"%s/data/sites.txt",m.datadir);
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
sscanf(line,"%4d %2s %lf %lf %f",
&id,abbrev,&lat,&lng,&alt);
strcpy(observer,line+38);
// Change to km
alt/=1000.0;
if (id==site_id) {
m.lat=lat;
m.lng=lng;
m.alt=alt;
m.site_id=id;
strcpy(m.observer,observer);
}
}
fclose(file);
return;
}
void initialize(void)
{
int i;
char *env;
// Default Map parameters
m.azi0=0;
m.alt0=90.0;
m.lat=0.0;
m.lng=0.0;
m.alt=0.0;
// Default settings
m.site_id=0;
// Get environment variables
env=getenv("ST_DATADIR");
if (env!=NULL) {
strcpy(m.datadir,env);
} else {
printf("ST_DATADIR environment variable not found.\n");
}
env=getenv("ST_COSPAR");
if (env!=NULL) {
get_site(atoi(env));
} else {
printf("ST_COSPAR environment variable not found.\n");
}
return;
}
// Return x modulo y [0,y)
double modulo(double x,double y)
{
x=fmod(x,y);
if (x<0.0) x+=y;
return x;
}
// 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;
}
// Convert Sexagesimal into Decimal
double sex2dec(char *s)
{
double x;
float deg,min,sec;
char t[LIM];
strcpy(t,s);
deg=fabs(atof(strtok(t," :")));
min=fabs(atof(strtok(NULL," :")));
sec=fabs(atof(strtok(NULL," :")));
x=(double) deg+(double) min/60.+(double) sec/3600.;
if (s[0]=='-') x= -x;
return x;
}
// Compute Julian Day from Date
double date2mjd(int year,int month,double day)
{
int a,b;
double jd;
if (month<3) {
year--;
month+=12;
}
a=floor(year/100.);
b=2.-a+floor(a/4.);
if (year<1582) b=0;
if (year==1582 && month<10) b=0;
if (year==1852 && month==10 && day<=4) b=0;
jd=floor(365.25*(year+4716))+floor(30.6001*(month+1))+day+b-1524.5;
return jd-2400000.5;
}
// Present nfd
void nfd_now(char *s)
{
time_t rawtime;
struct tm *ptm;
// Get UTC time
time(&rawtime);
ptm=gmtime(&rawtime);
sprintf(s,"%04d-%02d-%02dT%02d:%02d:%02d",ptm->tm_year+1900,ptm->tm_mon+1,ptm->tm_mday,ptm->tm_hour,ptm->tm_min,ptm->tm_sec);
return;
}
// nfd2mjd
double nfd2mjd(char *date)
{
int year,month,day,hour,min,sec;
double mjd,dday;
sscanf(date,"%04d-%02d-%02dT%02d:%02d:%02d",&year,&month,&day,&hour,&min,&sec);
dday=day+hour/24.0+min/1440.0+sec/86400.0;
mjd=date2mjd(year,month,dday);
return mjd;
}
// Convert horizontal into equatorial coordinates
void horizontal2equatorial(double mjd,double azi,double alt,double *ra,double *de)
{
double h;
h=atan2(sin(azi*D2R),cos(azi*D2R)*sin(m.lat*D2R)+tan(alt*D2R)*cos(m.lat*D2R))*R2D;
*ra=modulo(gmst(mjd)+m.lng-h,360.0);
*de=asin(sin(m.lat*D2R)*sin(alt*D2R)-cos(m.lat*D2R)*cos(alt*D2R)*cos(azi*D2R))*R2D;
if (*ra<0.0)
*ra+=360.0;
return;
}
// Convert equatorial into horizontal coordinates
void equatorial2horizontal(double mjd,double ra,double de,double *azi,double *alt)
{
double h;
h=gmst(mjd)+m.lng-ra;
*azi=modulo(atan2(sin(h*D2R),cos(h*D2R)*sin(m.lat*D2R)-tan(de*D2R)*cos(m.lat*D2R))*R2D,360.0);
*alt=asin(sin(m.lat*D2R)*sin(de*D2R)+cos(m.lat*D2R)*cos(de*D2R)*cos(h*D2R))*R2D;
return;
}
double parallactic_angle(double mjd,double ra,double de)
{
double h,q;
h=gmst(mjd)+m.lng-ra;
q=atan2(sin(h*D2R),(tan(m.lat*D2R)*cos(de*D2R)-sin(de*D2R)*cos(h*D2R)))*R2D;
return q;
}
// Convert Decimal into Sexagesimal
void dec2sex(double x,char *s,int f,int len)
{
int i;
double sec,deg,min;
char sign;
char *form[]={"::",",,","hms"," "};
sign=(x<0 ? '-' : ' ');
x=3600.*fabs(x);
sec=fmod(x,60.);
x=(x-sec)/60.;
min=fmod(x,60.);
x=(x-min)/60.;
// deg=fmod(x,60.);
deg=x;
if (len==7) sprintf(s,"%c%02i%c%02i%c%07.4f%c",sign,(int) deg,form[f][0],(int) min,form[f][1],sec,form[f][2]);
if (len==6) sprintf(s,"%c%02i%c%02i%c%06.3f%c",sign,(int) deg,form[f][0],(int) min,form[f][1],sec,form[f][2]);
if (len==5) sprintf(s,"%c%02i%c%02i%c%05.2f%c",sign,(int) deg,form[f][0],(int) min,form[f][1],sec,form[f][2]);
if (len==4) sprintf(s,"%c%02i%c%02i%c%04.1f%c",sign,(int) deg,form[f][0],(int) min,form[f][1],sec,form[f][2]);
if (len==2) sprintf(s,"%c%02i%c%02i%c%02i%c",sign,(int) deg,form[f][0],(int) min,form[f][1],(int) floor(sec),form[f][2]);
return;
}
int main(int argc,char *argv[])
{
int arg=0;
char sra[16],sde[16];
// Initialize
initialize();
// Get time
nfd_now(m.nfd);
m.mjd=nfd2mjd(m.nfd);
// Decode options
while ((arg=getopt(argc,argv,"t:R:D:A:E:h"))!=-1) {
switch(arg) {
case 't':
strcpy(m.nfd,optarg);
m.mjd=nfd2mjd(m.nfd);
break;
case 'R':
m.ra0=15.0*sex2dec(optarg);
strcpy(m.orientation,"equatorial");
break;
case 'D':
m.de0=sex2dec(optarg);
strcpy(m.orientation,"equatorial");
break;
case 'A':
m.azi0=modulo(atof(optarg)+180.0,360.0);
strcpy(m.orientation,"horizontal");
break;
case 'E':
m.alt0=atof(optarg);
strcpy(m.orientation,"horizontal");
break;
case 'h':
usage();
return 0;
break;
default:
usage();
return 0;
}
}
// Compute RA and Dec
if (strcmp(m.orientation,"horizontal")==0)
horizontal2equatorial(m.mjd,m.azi0,m.alt0,&m.ra0,&m.de0);
else if (strcmp(m.orientation,"equatorial")==0)
equatorial2horizontal(m.mjd,m.ra0,m.de0,&m.azi0,&m.alt0);
// Parallactic angle
m.q=parallactic_angle(m.mjd,m.ra0,m.de0);
// Get sexagesimal
dec2sex(m.ra0/15.0,sra,0,5);
dec2sex(m.de0,sde,0,4);
// Print to screen
printf("%s R:%s D:%s A: %4.1f E: %4.1f q: %5.2f\n",m.nfd,sra,sde,modulo(m.azi0-180.0,360.0),m.alt0,m.q);
// Send position
send_position(sra,sde);
return 0;
}

18
waitfor.c 100644
View File

@ -0,0 +1,18 @@
#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <math.h>
#include <sys/time.h>
int main(int argc,char *argv[])
{
struct timeval tv;
for (;;) {
gettimeofday(&tv,0);
if (tv.tv_usec>999000)
break;
}
return 0;
}