diff --git a/addwcs.c b/addwcs.c index 9e5b8e2..8de503e 100644 --- a/addwcs.c +++ b/addwcs.c @@ -8,7 +8,7 @@ #include #include -#define NMAX 2048 +#define NMAX 8192 #define LIM 128 #define D2R M_PI/180.0 #define R2D 180.0/M_PI diff --git a/calibrate.c b/calibrate.c index b3647e9..4144b31 100644 --- a/calibrate.c +++ b/calibrate.c @@ -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; diff --git a/data/desig.txt b/data/desig.txt index 3a6d1ac..11080f3 100644 --- a/data/desig.txt +++ b/data/desig.txt @@ -14819,4 +14819,5 @@ 39231 13037H 39232 13043A 99212 13706B +99068 13770A 99999 99000A diff --git a/data/sites.txt b/data/sites.txt index cc45b58..6c4fa8e 100644 --- a/data/sites.txt +++ b/data/sites.txt @@ -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 diff --git a/deproject.c b/deproject.c index 986daed..bf45dc0 100644 --- a/deproject.c +++ b/deproject.c @@ -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); diff --git a/jpg2fits.c b/jpg2fits.c index b7058ac..6a65eb0 100644 --- a/jpg2fits.c +++ b/jpg2fits.c @@ -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;i180.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; diff --git a/satid.c b/satid.c index e2aff17..a1b188b 100644 --- a/satid.c +++ b/satid.c @@ -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); diff --git a/satorbit.c b/satorbit.c index 04862ff..8f0762d 100644 --- a/satorbit.c +++ b/satorbit.c @@ -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)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.ps.psun) + cpgsci(7); + + // Plot + if (i==0) { + plot_footprint(s); + if (!(sqrt(x*x+y*y)XKMPER) { + if (sqrt(x*x+y*y) 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; diff --git a/skymap.c b/skymap.c index 3ba09dd..0d2a599 100644 --- a/skymap.c +++ b/skymap.c @@ -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); diff --git a/slewto.c b/slewto.c new file mode 100644 index 0000000..19f6a24 --- /dev/null +++ b/slewto.c @@ -0,0 +1,397 @@ +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +#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,"%s%s",sra,sde); + + // New packet style (as of 2013-08-20) + sprintf(packet,"%s%s",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; +} diff --git a/waitfor.c b/waitfor.c new file mode 100644 index 0000000..b8b61cc --- /dev/null +++ b/waitfor.c @@ -0,0 +1,18 @@ +#include +#include +#include +#include +#include + +int main(int argc,char *argv[]) +{ + struct timeval tv; + + for (;;) { + gettimeofday(&tv,0); + if (tv.tv_usec>999000) + break; + } + + return 0; +}