DODATEK 1
PROGRAM
WYZNACZAJĄCY MIEJSCE WIDOME
KSIĘŻYCA
(moon.c)
#include<conio.h>
#include<stdio.h>
#include<math.h>
void skyjpl(double
JD,double r2[9]);
void skymoon(double
r3[7]) /* miejsce widome Ksizyca
*/
{int i,ii,j;
double
alfa=0,delta=0,
JD=0,JD1=0,JD2=0,q[3],T,dt,T1,TDB_UT,
X,Y,Z,fi,S,S1,C,ff,kkk,X1,Y1,Z1,longitude,height,
VX,VY,VZ,VX1,VY1,VZ1,
EB1[3],P[3],Q[3],p[3],p1[3],E[3],e[3],V[3],p2[3],p4[3],
SM[3],SE[3],EM[3],
PP,EE,QQ,VV,pe,pq,qe,p1V,b,r2[9],
dpsi=0,depsil=0,epsil0=0,epsil=0,N[3][3],R[3][3],R1[3][3],r,
eta,zeta,teta,prs[3][3],RR;
height=r3[3];
longitude=r3[4];
fi=r3[5];
JD=JD1=r3[0];
T=(JD-2451545.0)/36525.0;
/****************************** NUTACJA
******************************/
r2[7]=3;r2[8]=14;
skyjpl(JD1,r2);
dpsi=r2[0];
depsil=r2[1];
/* dpsi=-12.602;
depsil=-6.710;
*/
epsil0= 84381.448
-(46.8150*T) -(0.00059*T*T) +(0.001813*T*T*T);
/*const in arcsec*/
epsil0=epsil0/3600.;
epsil0=epsil0*M_PI; epsil0=epsil0/180.;
epsil= epsil0 +
depsil; /* ok */
/********_nutation_matrix_*********/
N[0][0]=cos(dpsi);
N[1][0]=sin(dpsi)*cos(epsil);
N[2][0]=sin(dpsi)*sin(epsil);
N[0][1]=-sin(dpsi)*cos(epsil0);
N[1][1]=cos(dpsi)*cos(epsil)*cos(epsil0)+sin(epsil)*sin(epsil0);
N[2][1]=cos(dpsi)*sin(epsil)*cos(epsil0)-cos(epsil)*sin(epsil0);
N[0][2]=-sin(dpsi)*sin(epsil0);
N[1][2]=cos(dpsi)*cos(epsil)*sin(epsil0)-sin(epsil)*cos(epsil0);
N[2][2]=cos(dpsi)*sin(epsil)*sin(epsil0)+cos(epsil)*cos(epsil0);
/*********************** PRECESJA
************************************/
/* const.in arcsec */
eta =
2306.2182 * T + 0.30188 * T * T + 0.017998 * T * T * T ;
zeta =
2306.2181 * T + 1.09468 * T * T + 0.018203 * T * T * T ;
teta =
2004.3109 * T - 0.42665 * T * T - 0.041833 * T * T * T ;
/* transf.to rad */
eta =
eta * 4.84813681109536e-6 ;
zeta = zeta * 4.84813681109536e-6 ;
teta = teta * 4.84813681109536e-6 ;
prs[0][0]=
cos(eta)*cos(teta)*cos(zeta)-sin(eta)*sin(zeta);
prs[1][0]=
cos(eta)*cos(teta)*sin(zeta)+sin(eta)*cos(zeta);
prs[2][0]= sin(teta)*cos(eta);
prs[0][1]=-sin(eta)*cos(teta)*cos(zeta)-cos(eta)*sin(zeta);
prs[1][1]=-sin(eta)*cos(teta)*sin(zeta)+cos(eta)*cos(zeta);
prs[2][1]=-sin(teta)*sin(eta);
prs[0][2]=-cos(zeta)*sin(teta);
prs[1][2]=-sin(teta)*sin(zeta);
prs[2][2]= cos(teta);
/**********************
MACIERZ R=P.N i R1
***************************/
RR=0;
for(i=0;i<3;i++)
for(ii=0;ii<3;ii++)
{
for(j=0;j<3;j++)
{
RR=RR+((N[i][j])*(prs[j][ii]));
}
R[i][ii]=RR;
RR=0;
}
RR=R[0][0]*(R[1][1]*R[2][2]-R[2][1]*R[1][2])
-R[1][0]*(R[0][1]*R[2][2]-R[2][1]*R[0][2])
+R[2][0]*(R[0][1]*R[1][2]-R[1][1]*R[0][2]);
R1[0][0]=
(R[1][1]*R[2][2]-R[2][1]*R[1][2])/RR;
R1[0][1]=-(R[0][1]*R[2][2]-R[2][1]*R[0][2])/RR;
R1[0][2]=
(R[0][1]*R[1][2]-R[1][1]*R[0][2])/RR;
R1[1][0]=-(R[1][0]*R[2][2]-R[2][0]*R[1][2])/RR;
R1[1][1]= (R[0][0]*R[2][2]-R[2][0]*R[0][2])/RR;
R1[1][2]=-(R[0][0]*R[1][2]-R[1][0]*R[0][2])/RR;
R1[2][0]=
(R[1][0]*R[2][1]-R[2][0]*R[1][1])/RR;
R1[2][1]=-(R[0][0]*R[2][1]-R[2][0]*R[0][1])/RR;
R1[2][2]=
(R[0][0]*R[1][1]-R[1][0]*R[0][1])/RR;
/**************************
SIDERAL TIME *****************************/
JD2=floor(JD)+0.5;
T=(JD2-2451545.)/36525.;
S1=24110.54841 +
8640184.812866*T + 0.093104*T*T - (6.2E-6)*T*T*T;
S1=S1/3600.;
S1=S1+dpsi*cos(epsil)*180/(M_PI*15);
kkk=366.2422/365.2422;
T1=((JD-2378496.5)/36525.);
dt=JD-JD2;
S1=S1+dt*kkk*24-kkk*(
63.1839/3600.) +longitude/15.;
/*printf("\nS1=
%2.9lf\n",S1);*/
S1=24*(S1/24-floor(S1/24));
/*printf("\n\tS1
= %2.1lf h %2.1lf m %2.4lf
s\n\n",S1,(S1-floor(S1))*60,((S1-floor(S1))*60-floor((S1-floor(S1))*60))*60);
*/
/*********************************************************************/
/****** GEOCENTRYCZNA PREDKOSC I POLOZENIE
OBSERWATORA ******/
/*********************************************************************/
S1=S1*15*M_PI/180;
fi=fi*M_PI/180; ff=0.003352813177897;
C=1/(sqrt( cos(fi)*cos(fi) +
(1-ff)*(1-ff)*sin(fi)*sin(fi) ));
S=(1-ff)*(1-ff)*C;
X=((6378140.*C+height)*cos(fi)*cos(S1)) ;
Y=((6378140.*C+height)*cos(fi)*sin(S1)) ;
Z=((6378140.*S+height)*sin(fi)) ;
X=X / (1.49597870E+11);
Y=Y / (1.49597870E+11);
Z=Z / (1.49597870E+11);
X1= R1[0][0]*X +
R1[0][1]*Y + R1[0][2]*Z;
Y1= R1[1][0]*X +
R1[1][1]*Y + R1[1][2]*Z;
Z1= R1[2][0]*X + R1[2][1]*Y
+ R1[2][2]*Z;
VX=7.2921151467E-5
*(6378140.*C+height)*cos(fi)*sin(S1)*86400. / (1.49597870E+11);
VY=7.2921151467E-5
*(6378140.*C+height)*cos(fi)*cos(S1)*86400. / (1.49597870E+11);
VZ=0;
VX1= R1[0][0]*VX +
R1[0][1]*VY + R1[0][2]*VZ;
VY1= R1[1][0]*VX +
R1[1][1]*VY + R1[1][2]*VZ;
VZ1= R1[2][0]*VX +
R1[2][1]*VY + R1[2][2]*VZ;
/************************
CZAS ABERRACJI ***************************/
/****************
itteracyjne przyblizenie polozenia ***************/
i=0;
r2[7]=12;r2[8]=3;
skyjpl(JD,r2);
SE[0]=r2[0]+X1;
SE[1]=r2[1]+Y1;
SE[2]=r2[2]+Z1;
while(i<=4){ i++;
r2[7]=12;r2[8]=10;
skyjpl(JD1,r2);
SM[0]=r2[0];
SM[1]=r2[1];
SM[2]=r2[2];
EM[0]=SE[0]-SM[0];
EM[1]=SE[1]-SM[1];
EM[2]=SE[2]-SM[2];
r=sqrt( EM[0]*EM[0]
+ EM[1]*EM[1] + EM[2]*EM[2] );
JD1=JD-r/173.144633;
}
/*********** TOPOCENTRYCZNY WEKTOR POLOZENIA KSIEZYCA
P[i] ********/
r2[7]=11;r2[8]=10;
skyjpl(JD1,r2);
Q[0]=r2[0];
Q[1]=r2[1];
Q[2]=r2[2];
r2[7]=11;r2[8]=3;
skyjpl(JD,r2);
E[0]=r2[0]+X1;
E[1]=r2[1]+Y1;
E[2]=r2[2]+Z1;
P[0]=Q[0]-E[0];
P[1]=Q[1]-E[1];
P[2]=Q[2]-E[2];
r3[6]=(sqrt(P[0]*P[0]+P[1]*P[1]+P[2]*P[2]))*149597870;
/***
ODLEGLOSC !!!!! *****/
T=(JD-2451545.)/36525.;
/*********** e - unitvector Bary. Earth Pos. ***************/
EE=sqrt(E[0]*E[0]+E[1]*E[1]+E[2]*E[2]);
for(i=0;i<=2;i++)
e[i]=E[i]/EE;
/*********** p - unitvector Earth-Planet
***************/
PP=sqrt(P[0]*P[0]+P[1]*P[1]+P[2]*P[2]);
for(i=0;i<=2;i++)
p[i]=P[i]/PP;
/************ q - unitvector Helioc.
Planet Pos. ********/
QQ=sqrt(Q[0]*Q[0]+Q[1]*Q[1]+Q[2]*Q[2]);
for(i=0;i<=2;i++)
q[i]=Q[i]/QQ;
/******************** p.q q.p e.p
*******************************/
pe=p[0]*e[0]+p[1]*e[1]+p[2]*e[2];
pq=p[0]*q[0]+p[1]*q[1]+p[2]*q[2];
qe=q[0]*e[0]+q[1]*e[1]+q[2]*e[2];
/********** p1
***** LIGHT DEFLECTION
****************************/
for(i=0;i<=2;i++)
{p1[i]=p[i]+(2*9.87*EE*(e[i]*pq-pe*q[i]))/(1000000000.*(1+qe));}
/********** p2
***** ABERRACJA *************************************/
r2[7]=12;r2[8]=3;
skyjpl(JD,r2);
EB1[0]=r2[3]-VX1;
EB1[1]=r2[4]-VY1;
EB1[2]=r2[5]-VZ1;
for(i=0;i<=2;i++)
{V[i]=EB1[i]*0.0057755;}
VV=sqrt(V[0]*V[0]+V[1]*V[1]+V[2]*V[2]);
p1V=p1[0]*V[0]+p1[1]*V[1]+p1[2]*V[2];
b=pow(1-VV*VV,-0.5);
for(i=0;i<=2;i++){
/* plase God help me */
p2[i]=((p1[i]/b+V[i]*(1+p1V/(1/b+1))))/(1+p1V);
}
/**************** R
- precession & nutation matrix
***************/
p4[0]=R[0][0]*p2[0]+R[0][1]*p2[1]+R[0][2]*p2[2];
p4[1]=R[1][0]*p2[0]+R[1][1]*p2[1]+R[1][2]*p2[2];
p4[2]=R[2][0]*p2[0]+R[2][1]*p2[1]+R[2][2]*p2[2];
alfa=atan2(p4[1],p4[0]);
if(alfa<0)alfa=alfa+2*M_PI;
r3[1]=alfa;
delta=asin(p4[2]);
r3[2]=delta;
return;
}
DODATEK 2
PROGRAM
SŁUŻĄCY DO UTWORZENIA
KATALOGU ZODIAKALNEGO
(hipzod.c)
#include
<stdio.h>
#include
<conio.h>
#include
<string.h>
#include
<stdlib.h>
#include
<math.h>
void main(void)
{
FILE *u2,*u3;
char n2[15],n3[15];
char hdhipp[7],mag[6],paralaks[8];
int
nrHIPP;
int
alfag,alfam,deltast,deltam;
double
alfas,deltas;
double
al,del,paralaksa,dececl;
char bufor[500];
int uu,dd,ZNAK,licznik;
double cc,bb,ee/*sprawdzenie*/;
clrscr();
mag[5]=paralaks[7]='\0';
printf("\n Program HIPP
2000m\n");
/***********************************************************************/
strcpy(n2,"d:\\wd\\trnship1.dat");
if((u2=fopen(n2,"r"))==NULL)
{
printf("\n\nfile not
found\n");
exit(0);
}
/************************************************************/
strcpy(n3,"d:\\hip2000m.dat");
if((u3=fopen(n3,"w"))==NULL)
{
printf("\n\nfile not
found\n");
exit(0);
}
/******************************************************/
licznik=0;
while(fgets(bufor,500,u2)!=NULL)
{
sscanf(bufor,"%d %d %d %lf %d %d %lf"
,
&nrHIPP,&alfag,&alfam,&alfas,&deltast,&deltam,&deltas);
al=((double)alfag+(double)alfam/60.0+(double)alfas/3600.0)*15.0;
del=(fabs(deltast)+deltam/60.0+deltas/3600.0);
if(deltast>0)del=-del;
al=al*M_PI/180;
dececl=-23.45*sin(al);
if(del>(dececl-9)&&del<(dececl+9))
fputs(bufor,u3);
licznik++;
/*
if(licznik==1000) break; */
}
fclose (u2);
fclose (u3);
}
DODATEK 3
PROGRAM
TWORZĄCY PODKATALOG
#include<graphics.h> /* biblioteki */
#include<conio.h>
#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#include<string.h>
void skymoon(double
r3[6]);
void sstar(double
JD);
void efemain(void);
void efem(double
data,double szer,double dlug,double wys)
/***************************DEKLARACJA************************************/
{int ah,am,dd,dm;
double NR1;
double
alfa=0,delta=0,alfamoon,r3[6],alfamoon1,deltamoon1,deltamoon2[2]
,deltamoon,mialfa=0,midelta=0,JD,r5[4],
PI=0,vel=0,
mag,as,ds;
char bufor[81],plik[15];
FILE *f, *f1;
r5[0]=JD=data; r5[1]=r3[3]=wys; r5[2]=r3[4]=dlug; r5[3]=r3[5]=szer;
strcpy(plik,"h:\\HIP\\hip2000b.dat");
if((f=fopen(plik,"r"))==NULL){printf("\n\t no way
");exit(0);}
f1=fopen("stars1.dat","w");
r3[0]=JD=data+0.75;
skymoon(r3);
alfamoon=r3[1];
deltamoon=r3[2];
r3[0]=data+1.25;
skymoon(r3);
alfamoon1=r3[1];
deltamoon1=r3[2];
deltamoon2[0]=deltamoon;
deltamoon2[1]=deltamoon;
if(deltamoon1>=deltamoon)deltamoon2[0]=deltamoon1;
if(deltamoon1<deltamoon) deltamoon2[1]=deltamoon1;
printf("\nskrajne polozenia ksiezyca\nalfa %lf delta
%lf",alfamoon,deltamoon2[0]);
printf("\nalfa %lf delta
%lf\n",alfamoon1,deltamoon2[1]);
printf("\nProgram tworzy
podkatalog Hipparcos'a \n");
alfa=0;
NR1=0.;
while(!feof(f))
{
fgets(bufor,81,f);
sscanf(bufor,"%lf %d %d %lf %d %d %lf %lf %lf %lf %lf %lf",
&NR1,&ah,&am,&as,&dd,&dm,&ds,&mialfa,&midelta
,&mag ,&PI,&vel);
alfa= 15*ah *M_PI/180 +15 *am*M_PI/10800.;
if(alfa<(alfamoon1+0.1))
{
delta=fabs(dd) *M_PI/180 + dm*M_PI/10800.;
if(dd<0)delta=-delta;
if(
alfa>alfamoon-0.006
&&
alfa<alfamoon1+0.006
&& delta<deltamoon2[0]+0.006
&& delta>deltamoon2[1]-0.006)
{
fputs(bufor,f1);
printf(".");
}
}
}
fclose(f);fclose(f1);
printf(" END\n");
printf("\nliczenie miejsc widomych");
sstar(JD); efemain();
}
DODATEK 4
PROGRAM
LICZĄCY MIEJSCA WIDOME
GWIAZD (star.c)
#include<graphics.h> /* biblioteki */
#include<conio.h>
#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#include<string.h>
void skyjpl(double
JD,double r2[9]);
void sstar(double
JD) /* miejsca widome */
/***************************DEKLARACJA************************************/
{int i,ii,j,ah,am,dd,dm,NR;
double alfa=0,delta=0,r2[9],NR1=0
,mialfa=0,midelta=0,
PI=0,vel=0,q[3],m[3],T,
SB[3],EB[3],EB1[3],P[3],p[3],p1[3],E[3],e[3],V[3],p2[3],
p4[3],PP,EE,VV,pe,p1V,b,
mag,as,ds,dpsi=0,depsil=0,epsil0=0,epsil=0,N[3][3],R[3][3],ahh,amm,dmm,
ddd,eta,zeta,teta,prs[3][3],RR;
char bufor[81],plik[15];
FILE *f,*f1;
/*****************************************************************/
/*
strcpy(plik,"c:\\c\\turboc\\newfk5.jof"); */
/*
strcpy(plik,"h:\\praca\\newfk5.jof"); */
strcpy(plik,"h:\\skyefem\\stars1.dat");
if((f=fopen(plik,"r"))==NULL){printf("\n\t no way
/sstar line64/");exit(0);}
f1=fopen("stars2.dat","w");
while(fgets(bufor,100,f)!=NULL) /*POPRAWKA WARUNKU*/
{
sscanf(bufor,"%lf %d %d %lf %d %d %lf %lf %lf %lf %lf %lf",
&NR1,&ah,&am,&as,&dd,&dm,&ds,&mialfa,&midelta
,&mag ,&PI,&vel);
alfa=
15*ah *M_PI/180 +15
*am*M_PI/10800. +as*15 *M_PI/648000. ;
delta= fabs(dd) *M_PI/180 +
dm*M_PI/10800. +ds*
M_PI/648000. ;
if(dd<0)delta=-delta;
mialfa=
mialfa *15*M_PI/64800000.;
midelta= midelta *M_PI/6480000.;
PI=PI*M_PI/648000.;
vel=vel/0.47404707978472;
/**********************************************************************/
/********************** liczenie miejsca
widomego *********************/
/**********************************************************************/
q[0]=cos(alfa)*cos(delta); /* Rectangular coordinates */
q[1]=sin(alfa)*cos(delta);
q[2]=sin(delta);
/* space
motion vector */
m[0]=-mialfa*cos(delta)*sin(alfa)-midelta*sin(delta)*cos(alfa)+PI*vel*cos(delta)*cos(alfa);
m[1]=
mialfa*cos(delta)*cos(alfa)-midelta*sin(delta)*sin(alfa)+PI*vel*cos(delta)*sin(alfa);
m[2]= midelta*cos(delta)+PI*vel*sin(delta);
T=(JD-2451545.)/36525.;/*printf("\n\t T
= %2.8lf",T);printf("\n");
*/
/*.....................................................................*/
/*
soma(JD,E,SB,EB1);
for(i=0;i<=2;i++)SB[i]=SB[i]*(-1);
for(i=0;i<=2;i++)EB[i]=E[i]+SB[i];
*/
r2[7]=11; r2[8]=12;
skyjpl(JD,r2);
SB[0]=-r2[0];
SB[1]=-r2[1];
SB[2]=-r2[2];
r2[7]=11; r2[8]=3;
skyjpl(JD,r2);
E[0]=r2[0];
E[1]=r2[1];
E[2]=r2[2];
r2[7]=12; r2[8]=3;
skyjpl(JD,r2);
EB1[0]=r2[3];
EB1[1]=r2[4];
EB1[2]=r2[5];
for(i=0;i<=2;i++)EB[i]=E[i]+SB[i];
/*.....................................................................*/
/* geocentric vector of star */
for(i=0;i<=2;i++){P[i]=q[i]+T*m[i]-PI*EB[i];
/*printf("\n\tP %d = %2.9lf
",i,P[i]);*/}
PP=sqrt(P[0]*P[0]+P[1]*P[1]+P[2]*P[2]);
for(i=0;i<=2;i++){p[i]=P[i]/PP;
/*
printf("\n\tp %d = %2.9lf ",i,p[i]);*/ }
/******************** e
**************************************/
EE=sqrt(E[0]*E[0]+E[1]*E[1]+E[2]*E[2]);
for(i=0;i<=2;i++){e[i]=E[i]/EE;}
/********************* p1
*************************************/
/* geocentric direction of star corrected for
light deflection */
pe=p[0]*e[0]+p[1]*e[1]+p[2]*e[2];
for(i=0;i<=2;i++)
{p1[i]=p[i]+(2*9.87*EE*(e[i]-pe*p[i]))/(1000000000.*(1+pe));
}
/********************* p2
*************************************/
/* proper direction in the geocentric natural
frame */
for(i=0;i<=2;i++){V[i]=EB1[i]*0.0057755;/*printf("\n\tV
%d = %2.9lf ",i,V[i]);*/}
VV=sqrt(V[0]*V[0]+V[1]*V[1]+V[2]*V[2]);
p1V=p1[0]*V[0]+p1[1]*V[1]+p1[2]*V[2];
b=pow(1-VV*VV,-0.5);
for(i=0;i<=2;i++){
/* plase God help me */
p2[i]=((p1[i]/b+V[i]*(1+p1V/(1/b+1))))/(1+p1V);
/*
printf("\n\tp2
%d = %2.9lf ",i,p2[i]);*/}
/*********************_NUTATION_*****************************************/
r2[7]=3; r2[8]=14;
skyjpl(JD,r2);
dpsi=r2[0];
depsil=r2[1];
epsil0=84381.448-(46.8150*T)-(0.00059*T*T)+(0.001813*T*T*T); /*const in arcsec*/
epsil0=epsil0/3600;epsil0=epsil0*M_PI;epsil0=epsil0/180;
epsil=
epsil0+depsil; /* ok */
/********_nutation_matrix_*********/
N[0][0]=cos(dpsi);
N[1][0]=sin(dpsi)*cos(epsil);
N[2][0]=sin(dpsi)*sin(epsil);
N[0][1]=-sin(dpsi)*cos(epsil0);
N[1][1]=cos(dpsi)*cos(epsil)*cos(epsil0)+sin(epsil)*sin(epsil0);
N[2][1]=cos(dpsi)*sin(epsil)*cos(epsil0)-cos(epsil)*sin(epsil0);
N[0][2]=-sin(dpsi)*sin(epsil0);
N[1][2]=cos(dpsi)*cos(epsil)*sin(epsil0)-sin(epsil)*cos(epsil0);
N[2][2]=cos(dpsi)*sin(epsil)*sin(epsil0)+cos(epsil)*cos(epsil0);
/************************_end_of_nutation_****************************/
/***********************_precession_**********************************/
/* const.in arcsec */
eta =
2306.2182 * T + 0.30188 * T * T + 0.017998 * T * T * T ;
zeta =
2306.2181 * T + 1.09468 * T * T + 0.018203 * T * T * T ;
teta =
2004.3109 * T - 0.42665 * T * T - 0.041833 * T * T * T ;
/* transf.to rad */
eta =
eta * 0.00000484813681109536 ;
zeta = zeta * 0.00000484813681109536
;
teta = teta * 0.00000484813681109536
;
prs[0][0]=
cos(eta)*cos(teta)*cos(zeta)-sin(eta)*sin(zeta);
prs[1][0]=
cos(eta)*cos(teta)*sin(zeta)+sin(eta)*cos(zeta);
prs[2][0]= sin(teta)*cos(eta);
prs[0][1]=-sin(eta)*cos(teta)*cos(zeta)-cos(eta)*sin(zeta);
prs[1][1]=-sin(eta)*cos(teta)*sin(zeta)+cos(eta)*cos(zeta);
prs[2][1]=-sin(teta)*sin(eta);
prs[0][2]=-cos(zeta)*sin(teta);
prs[1][2]=-sin(teta)*sin(zeta);
prs[2][2]= cos(teta);
/***************************END***********************************/
/*******************_macierz_R=P.N_*******************************/
RR=0;
for(i=0;i<3;i++)
for(ii=0;ii<3;ii++)
{
for(j=0;j<3;j++)
{
RR=RR+((N[i][j])*(prs[j][ii]));
}
R[i][ii]=RR;
RR=0;
}
p4[0]=R[0][0]*p2[0]+R[0][1]*p2[1]+R[0][2]*p2[2];
p4[1]=R[1][0]*p2[0]+R[1][1]*p2[1]+R[1][2]*p2[2];
p4[2]=R[2][0]*p2[0]+R[2][1]*p2[1]+R[2][2]*p2[2];
/***************************END*************************************/
if( fabs(alfa-atan2(p4[1],p4[0]) )>1 )
alfa=atan2(p4[1],p4[0])+2*M_PI
;
else alfa=atan2(p4[1],p4[0]);
delta=asin(p4[2]);
ahh=amm=ddd=dmm=0;
ah=(alfa*180)/(3.141592654*15);
ahh=(alfa*180)/(3.141592654*15);
am=(ahh-ah)*60;
amm=(ahh-ah)*60;
as=(amm-am)*60;
dd=(delta*180)/3.141592654;
ddd=(delta*180)/3.141592654;
dm=ddd*60-dd*60;
dmm=(ddd-dd)*60;
ds=(dmm-dm)*60;
if(dd<0){dm=-dm;ds=-ds;}
printf(".");
fprintf(f1," %6.0lf %2.2lf
%d %d %2.4lf %d %d %2.3lf %1.8lf %1.8lf\n",
NR1,mag,ah,am,as,dd,dm,ds,alfa,delta);
} fclose(f);
fclose(f1);
}
DODATEK 5
PROGRAM
LICZĄCY FAZE I
ELONGACJĘ
(phase.c)
#include<graphics.h> /* biblioteki */
#include<conio.h>
#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#include<string.h>
void skyjpl(double
JD,double r2[9]);
void skyphase(double
JD,double r4[3])
{int i,ii,j;
double
eta,zeta,teta,JD1,dpsi,depsil,epsil0,epsil,prs[3][3],T,N[3][3],R[3][3],
R1[3][3],RR,P[3],Q[3],E[3],r2[9],PP,QQ,EE,b,C,S,X,Y,Z,longitude,JD2,S1,kkk
,T1,TDB_UT,dt,fi=52.3986,ff,height,latitude;
JD1=JD;
T=(JD-2451545.0)/36525.0;
/****************************** NUTACJA
******************************/
r2[7]=3;r2[8]=14;
skyjpl(JD1,r2);
dpsi=r2[0];
depsil=r2[1];
epsil0= 84381.448
-(46.8150*T) -(0.00059*T*T) +(0.001813*T*T*T);
/*const in arcsec*/
epsil0=epsil0/3600.;
epsil0=epsil0*M_PI; epsil0=epsil0/180.;
epsil= epsil0 +
depsil; /* ok */
/********_nutation_matrix_*********/
N[0][0]=cos(dpsi);
N[1][0]=sin(dpsi)*cos(epsil);
N[2][0]=sin(dpsi)*sin(epsil);
N[0][1]=-sin(dpsi)*cos(epsil0);
N[1][1]=cos(dpsi)*cos(epsil)*cos(epsil0)+sin(epsil)*sin(epsil0);
N[2][1]=cos(dpsi)*sin(epsil)*cos(epsil0)-cos(epsil)*sin(epsil0);
N[0][2]=-sin(dpsi)*sin(epsil0);
N[1][2]=cos(dpsi)*cos(epsil)*sin(epsil0)-sin(epsil)*cos(epsil0);
N[2][2]=cos(dpsi)*sin(epsil)*sin(epsil0)+cos(epsil)*cos(epsil0);
/*********************** PRECESJA
******************************/
/* const.in arcsec */
eta =
2306.2182 * T + 0.30188 * T * T + 0.017998 * T * T * T ;
zeta =
2306.2181 * T + 1.09468 * T * T + 0.018203 * T * T * T ;
teta =
2004.3109 * T - 0.42665 * T * T - 0.041833 * T * T * T ;
/* transf.to rad */
eta =
eta * 4.84813681109536e-6 ;
zeta = zeta * 4.84813681109536e-6 ;
teta = teta * 4.84813681109536e-6 ;
prs[0][0]=
cos(eta)*cos(teta)*cos(zeta)-sin(eta)*sin(zeta);
prs[1][0]=
cos(eta)*cos(teta)*sin(zeta)+sin(eta)*cos(zeta);
prs[2][0]= sin(teta)*cos(eta);
prs[0][1]=-sin(eta)*cos(teta)*cos(zeta)-cos(eta)*sin(zeta);
prs[1][1]=-sin(eta)*cos(teta)*sin(zeta)+cos(eta)*cos(zeta);
prs[2][1]=-sin(teta)*sin(eta);
prs[0][2]=-cos(zeta)*sin(teta);
prs[1][2]=-sin(teta)*sin(zeta);
prs[2][2]= cos(teta);
/**********************
MACIERZ R=P.N i R1
***************************/
RR=0;
for(i=0;i<3;i++)
for(ii=0;ii<3;ii++)
{
for(j=0;j<3;j++)
{
RR=RR+((N[i][j])*(prs[j][ii]));
}
R[i][ii]=RR;
RR=0;
}
RR=R[0][0]*(R[1][1]*R[2][2]-R[2][1]*R[1][2])
-R[1][0]*(R[0][1]*R[2][2]-R[2][1]*R[0][2])
+R[2][0]*(R[0][1]*R[1][2]-R[1][1]*R[0][2]);
R1[0][0]=
(R[1][1]*R[2][2]-R[2][1]*R[1][2])/RR;
R1[0][1]=-(R[0][1]*R[2][2]-R[2][1]*R[0][2])/RR;
R1[0][2]=
(R[0][1]*R[1][2]-R[1][1]*R[0][2])/RR;
R1[1][0]=-(R[1][0]*R[2][2]-R[2][0]*R[1][2])/RR;
R1[1][1]=
(R[0][0]*R[2][2]-R[2][0]*R[0][2])/RR;
R1[1][2]=-(R[0][0]*R[1][2]-R[1][0]*R[0][2])/RR;
R1[2][0]=
(R[1][0]*R[2][1]-R[2][0]*R[1][1])/RR;
R1[2][1]=-(R[0][0]*R[2][1]-R[2][0]*R[0][1])/RR;
R1[2][2]=
(R[0][0]*R[1][1]-R[1][0]*R[0][1])/RR;
/**************************
SIDERAL TIME **************************/
JD2=floor(JD)+0.5;
T=(JD2-2451545.)/36525.;
S1=24110.54841 +
8640184.812866*T + 0.093104*T*T - (6.2E-6)*T*T*T;
S1=S1/3600.;
S1=S1+dpsi*cos(epsil)*180/(M_PI*15);
kkk=366.2422/365.2422;
T1=((JD-2378496.5)/36525.);
dt=JD-JD2;
S1=S1+dt*kkk*24-kkk*(
64/3600.) +longitude/15.;
/*printf("\nS1=
%2.9lf\n",S1);*/
S1=24*(S1/24-floor(S1/24));
/*printf("\n\tS1
= %2.1lf h %2.1lf m %2.4lf
s\n\n",S1,(S1-floor(S1))*60,((S1-floor(S1))*60-floor((S1-floor(S1))*60))*60);
*/
/*******************************************************************/
/****** GEOCENTRYCZNA PREDKOSC I POLOZENIE
OBSERWATORA ****/
/*******************************************************************/
S1=S1*15*M_PI/180;
fi=fi*M_PI/180;
ff=0.003352813177897;
C=1/(sqrt( cos(fi)*cos(fi) +
(1-ff)*(1-ff)*sin(fi)*sin(fi) ));
S=(1-ff)*(1-ff)*C;
/* S1=4.44038606503932;*/
X=((6378140.*C+height)*cos(fi)*cos(S1)) ;
Y=((6378140.*C+height)*cos(fi)*sin(S1)) ;
Z=((6378140.*S+height)*sin(fi)) ;
X=X / (1.49597870E+11);
Y=Y / (1.49597870E+11);
Z=Z / (1.49597870E+11);
r2[7]=10;r2[8]=11;
skyjpl(JD,r2);
Q[0]=r2[0]; Q[1]=r2[1];
Q[2]=r2[2];
r2[7]=10;r2[8]=3;
skyjpl(JD,r2);
E[0]=r2[0]+X;
E[1]=r2[1]+Y;
E[2]=r2[2]+Z;
r2[7]=11;r2[8]=3;
skyjpl(JD,r2);
P[0]=r2[0]; /* -0.060418315; */
P[1]=r2[1]; /*
0.900587724; */
P[2]=r2[2]; /*
0.390463880; */
PP=sqrt( P[0]*P[0] +
P[1]*P[1] + P[2]*P[2] );
QQ=sqrt( Q[0]*Q[0] +
Q[1]*Q[1] + Q[2]*Q[2] );
EE=sqrt( E[0]*E[0] +
E[1]*E[1] + E[2]*E[2] );
/*
printf("\n\nzimia
slonce=%lf",PP);
printf("\n\nksiezyc
slonce=%lf",QQ);
printf("\n\nksiezyc
ziemia=%2.9lf\n\n",EE);
*/
r4[2]=b=acos((QQ*QQ+EE*EE-PP*PP)/(2*QQ*EE));
b=(1+cos(b))/2;
r4[0]=b;
r4[1]=EE*149597870.;
return;
}
DODATEK 6
PROGRAM
EFEMERYDALNY
(efemain.c)
#include<graphics.h> /* biblioteki */
#include<conio.h>
#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#include<string.h>
void skymoon(double
r3[7]);
void skyphase(double
JD, double r4[3]);
void efemain(void)
/***************************
DEKLARACJA ********************************/
{int i,ii,nr,j,ah,am,dd,dm,sd,dnn=0,elon;
int
licznik,dzien,miesiac,rok,godzina1,minuta1;
double godzina,minuta,sekunda;
double
alfa=0,delta=0,alfamoon,r3[7],alfamoon1,deltamoon1,deltamoon2[2],r4[3]
,deltamoon,da,dr,dr1,db,NR1=0,distance,faza,promien,promien1,
oldstat[2],newstat[2],medstat[2],
q[3],m[3],T,x,y,JD,longitude,altitude,height,data,
SB[3],EB[3],EB1[3],P[3],p[3],p1[3],E[3],e[3],V[3],p2[3],p3[3],p4[3],PP,EE,VV,pe,p1V,b,
mag,as,ds,N[3][3],R[3][3],ahh,amm,dmm,
ddd,eta,zeta,teta,prs[3][3],RR,RRR[3][3],gwiazdy1[200][6],PA,
gwiazdy3[200][6],gwiazdy2[200][2];
char bufor[81],plik[15];
FILE *f, *f1, *f2, *f3;
/*************************************************************/
strcpy(plik,"h:\\skyefem\\stars2.dat");
if((f=fopen(plik,"r"))==NULL){printf("\n\t no way { EFEM
line 41}");exit(0);}
/******************_wczytanie_danych_gwiazd_****************************/
nr=0;
while( fgets(bufor,81,f)!=NULL)
{
nr++;
sscanf(bufor,"%lf %lf %d %d %lf %d %d %lf %lf %lf",
&NR1,&mag,&ah,&am,&as,&sd,&dm,&ds,&alfa,&delta);
gwiazdy2[nr][0]=NR1;
gwiazdy1[nr][0]=mag;
gwiazdy1[nr][2]=delta;
gwiazdy1[nr][1]=alfa;
gwiazdy3[nr][0]=ah;
gwiazdy3[nr][1]=am;
gwiazdy3[nr][2]=as;
gwiazdy3[nr][3]=sd;
gwiazdy3[nr][4]=dm;
gwiazdy3[nr][5]=ds;
printf("\n%lf %lf %lf %lf", NR1,mag,alfa,delta);
}
fclose(f);
printf("\n\nliczba=%d\n\n",nr); getch();
f2=fopen("dmy.dat","r");
fscanf(f2," %d %d %d
%lf",&dzien,&miesiac,&rok,&data);
fclose(f2);
f3=fopen("place.dat","r");
fscanf(f3," %lf %lf
%lf",&height,&longitude,&altitude);
fclose(f3);
/************************ stan poczatkowy
****************************/
/* zerowanie statusu */
for(i=1;i<=nr;i++){ gwiazdy2[i][1]=0; }
/* wsp. pocz. moon */
r3[0]=JD=data + 0.75;
r3[3]=height;
r3[4]=longitude;
r3[5]=altitude;
skymoon(r3);
alfamoon=r3[1];
deltamoon=r3[2];
distance=r3[6];
/*
promien=atan(1738/distance); */
promien=asin(1738/distance);
printf("\na=%lf\nd=%lf\n\n",alfamoon,deltamoon);
getch();
/* okreslenie statusu pocz. */
licznik=1;
while(licznik<=nr)
{ /*
da=(alfamoon-gwiazdy1[licznik][1])*cos(gwiazdy1[licznik][2]);
dd=(deltamoon-gwiazdy1[licznik][2]);
dr=sqrt( da*da + dd*dd ); */
da=M_PI/2-deltamoon;
db=M_PI/2-gwiazdy1[licznik][2];
dr=acos(cos(da)*cos(db)+sin(da)*sin(db)*cos(alfamoon-gwiazdy1[licznik][1]));
if(dr<=promien /*0.004654211338652 0.00436332313*/ )gwiazdy2[licznik][1]=1;
licznik++;
}
for(i=1;i<=nr;i++){printf("\n%d
nr%6.0lf status%1.0lf",i, gwiazdy2[i][0], gwiazdy2[i][1]); getch();}
f1=fopen("efem.dat","w");
/************************ sprawdzenie czy
jest zakrycie **************/
while(r3[0]<(JD+0.25))
{
r3[0]=r3[0]+0.001 /*0.00001157407407407*/ ;
printf(".");
if(
(r3[0]-floor(r3[0]))>0.5 && dnn==0 ){ dzien++; dnn=1; }
skymoon(r3);
alfamoon=r3[1];
deltamoon=r3[2];
distance=r3[6];
/* promien=atan(1738/distance); */
promien=asin(1738/distance);
licznik=1;
while(licznik<=nr)
{
da=M_PI/2-deltamoon;
db=M_PI/2-gwiazdy1[licznik][2];
dr=acos(cos(da)*cos(db)+sin(da)*sin(db)*cos(alfamoon-gwiazdy1[licznik][1]));
if((dr<=promien
&& gwiazdy2[licznik][1]==0)||(dr>promien && gwiazdy2[licznik][1]==1) )
{
printf("*");
oldstat[0]=r3[0]-0.001;
oldstat[1]=gwiazdy2[licznik][1];
newstat[0]=r3[0];
if(gwiazdy2[licznik][1]==0) newstat[1]=1;
if(gwiazdy2[licznik][1]==1) newstat[1]=0;
for(i=0;i<=8;i++)
{
r3[0]=medstat[0]=(oldstat[0]+newstat[0])/2;
skymoon(r3);
alfamoon=r3[1];
deltamoon=r3[2];
distance=r3[6];
/* promien=atan(1738/distance);*/
promien=asin(1738/distance);
da=M_PI/2-deltamoon;
db=M_PI/2-gwiazdy1[licznik][2];
dr1=acos(cos(da)*cos(db)+sin(da)*sin(db)*cos(alfamoon-gwiazdy1[licznik][1]));
medstat[1]=0;
if(dr1<=promien )medstat[1]=1;
if(medstat[1]==oldstat[1])
{
newstat[1]=medstat[1];
newstat[0]=medstat[0];
}
else {
oldstat[1]=medstat[1];
oldstat[0]=oldstat[0]; }
} /* for */
skyphase(r3[0],r4);
elon=180-(r4[2]*180)/M_PI;
faza=r4[0];
PA=atan2( (alfamoon-gwiazdy1[licznik][1]),(deltamoon-gwiazdy1[licznik][2])
);
PA=PA*180/M_PI+180;
godzina1=godzina=(r3[0]-JD)*24+18;
if(godzina>24){godzina1=godzina=godzina-24;}
minuta1=minuta=(godzina-godzina1)*60;
sekunda=(minuta-minuta1)*60;
godzina=floor(godzina);
minuta =floor(minuta );
if(dr<=promien &&
gwiazdy2[licznik][1]==0)
{
printf("\n%2d %2d %4d %2.0lf %2.0lf
%2.0lf %6.0lf %lf zak %2.0lf"
,dzien,miesiac,rok,godzina,minuta,sekunda,gwiazdy2[licznik][0],r3[0],PA);
fprintf(f1,"\n%2d %2d %4d %2.0lf %2.0lf %2.0lf D %6.0lf %2.1lf
%lf %1.2lf %d %2.0lf %2.0lf %2.0lf %2.3lf %2.0lf %2.0lf %2.3lf"
,dzien,miesiac,rok,godzina,minuta,sekunda,gwiazdy2[licznik][0],gwiazdy1[licznik][0],r3[0],faza,elon,PA,
gwiazdy3[licznik][0],gwiazdy3[licznik][1],gwiazdy3[licznik][2],gwiazdy3[licznik][3],gwiazdy3[licznik][4],gwiazdy3[licznik][5]);
gwiazdy2[licznik][1]=1;
}
if(dr>promien &&
gwiazdy2[licznik][1]==1)
{
printf("\n%2d %2d %4d %2.0lf %2.0lf %2.0lf %6.0lf %lf odk
%2.0lf"
,dzien,miesiac,rok,godzina,minuta,sekunda,gwiazdy2[licznik][0],r3[0],PA);
fprintf(f1,"\n%2d %2d %4d %2.0lf %2.0lf %2.0lf R %6.0lf %2.1lf
%lf %1.2lf %d %2.0lf %2.0lf %2.0lf %2.3lf %2.0lf %2.0lf %2.3lf"
,dzien,miesiac,rok,godzina,minuta,sekunda,gwiazdy2[licznik][0],gwiazdy1[licznik][0],r3[0],faza,elon,PA,
gwiazdy3[licznik][0],gwiazdy3[licznik][1],gwiazdy3[licznik][2],gwiazdy3[licznik][3],gwiazdy3[licznik][4],gwiazdy3[licznik][5]);
gwiazdy2[licznik][1]=0;
}
} /* if */
licznik++;
}
}
fclose(f1);
printf("\n END ");
}
/**********THE END**********/
DODATEK 7
EFEMERYDA
(23.III.99r.)
23 4 1999
18 33 8 D 45043 8.5 2451292.273008 0.61
102 34 9 10 24.567 15 35 21.093
23 4 1999
18 46 6 R 44942 8.2 2451292.282010 0.61
102 248 9 9
10.830 15 14 38.257
23 4 1999
18 53 18 R 45043 8.5
2451292.287014 0.61 102
2 9 10 24.567 15 35 21.093
23 4 1999
19 3 23 R 44985 8.4 2451292.294021 0.61
102 303 9 9
44.652 15 26 56.312
23 4 1999
19 39 24 D 45170 6.5
2451292.319029 0.61 102
148 9 12 11.513 15
0 7.341
23 4 1999
20 34 8 R 45170 6.5 2451292.357031 0.61
103 248 9 12 11.513 15 0 7.341
23 4 1999
21 5 48 D 45410 5.4 2451292.379033 0.62
103 107 9 15 10.683 14 56 38.70
23 4 1999
21 27 25 D 45462 7.8
2451292.394041 0.62 103
82 9 15 49.708 15
0 41.112
23 4 1999
21 31 45 D 45474 8.6
2451292.397049 0.62 103
106 9 15 56.244 14 53 14.080
23 4 1999
21 50 28 D 45495 8.6
2451292.410051 0.62 103
39 9 16 8.397 15
7 24.895
23 4 1999
22 10 38 R 45495 8.6
2451292.424053 0.62 103
2 9 16 8.397 15
7 24.895
23 4 1999
22 12 5 R 45410 5.4 2451292.425055 0.62
104 289 9 15 10.683 14 56 38.709
23 4 1999
22 25 3 R 45462 7.8 2451292.434063 0.62
104 315 9 15 49.708 15 0 41.112
23 4 1999
22 36 35 R 45474 8.6
2451292.442070 0.62 104
289 9 15 56.244 14 53 14.080