Omowienie programu wkrotce
/*chi2.c*/
#include<stdio.h>
#include<math.h>
#include<stdlib.h>
#define DZIEN 86400 /*liczba sekund w 1 dniu*/
#define N 2 /*liczba planet*/
#define LP 11 /*liczba parametrow*/
#define M 38 /*liczba pomiarow*/
#define B 0.0001 /*dokladnosc rozwiazania r-nia Keplera*/
#define KROK 0.001 /*podzial siatki argumentow*/
#define ZAKRES 0.05 /*zakres zmian parametrow*/
#define PI 3.141592653589793115997963468544185161590576172 /*liczba pi*/
float NR(float P, float tau, float t, float e, float b, float pi);
float f(float s[], float t[], float v[], float d[], int Npom, int Lpar, float b, int Lplan, float pi);
main()
{
int i,j;
float Fn;
float V0,K[N],P[N],e[N],w[N],tau[N];
float s[LP],s0[LP],t[M],v[M],d[M];
FILE *plik, *plik1, *plik2, *plik3, *plik4, *plik5, *plik6, *plik7, *plik8;
FILE *plik9, *plik10, *plik11, *wynik;
char *nazwa;
float param[LP]={1,1,DZIEN,1,PI/180,DZIEN,1,DZIEN,1,PI/180,DZIEN};
/*wczytanie danych*/
nazwa="HD114729.dat";
plik=fopen(nazwa,"r");
for (i=0; i<=M-1; ++i)
{
fscanf(plik,"%f%f%f",&t[i],&v[i],&d[i]);
t[i]=t[i]*DZIEN;
}
fclose(plik);
/*parametry poczatkowe*/
V0=-5.32061;
K[0]=15.29573;
P[0]=1088.65051*DZIEN;
e[0]=0.205788;
w[0]=75.49266*PI/180;
tau[0]=509.81573*DZIEN;
K[1]=5.69828;
P[1]=10.52720*DZIEN;
e[1]=0.016552;
w[1]=38.30322*PI/180;
tau[1]=100.91705*DZIEN;
nazwa="HD114729_chi2_V0_1.dat";
plik1=fopen(nazwa,"w");
nazwa="HD114729_chi2_K1_1.dat";
plik2=fopen(nazwa,"w");
nazwa="HD114729_chi2_P1_1.dat";
plik3=fopen(nazwa,"w");
nazwa="HD114729_chi2_e1_1.dat";
plik4=fopen(nazwa,"w");
nazwa="HD114729_chi2_w1_1.dat";
plik5=fopen(nazwa,"w");
nazwa="HD114729_chi2_tau1_1.dat";
plik6=fopen(nazwa,"w");
nazwa="HD114729_chi2_K2_1.dat";
plik7=fopen(nazwa,"w");
nazwa="HD114729_chi2_P2_1.dat";
plik8=fopen(nazwa,"w");
nazwa="HD114729_chi2_e2_1.dat";
plik9=fopen(nazwa,"w");
nazwa="HD114729_chi2_w2_1.dat";
plik10=fopen(nazwa,"w");
nazwa="HD114729_chi2_tau2_1.dat";
plik11=fopen(nazwa,"w");
/* glowna czesc programu*/
s0[0]=V0;
for (j=0; j<=N-1; ++j)
{
s0[5*j+1]=K[j];
s0[5*j+2]=P[j];
s0[5*j+3]=e[j];
s0[5*j+4]=w[j];
s0[5*j+5]=tau[j];
}
for (i=0; i<=LP-1; ++i)
{
for (j=0; j<=LP-1; ++j)
{
s[j]=s0[j];
}
s[i]=(1-ZAKRES)*s0[i];
if (i==0) wynik=plik1;
if (i==1) wynik=plik2;
if (i==2) wynik=plik3;
if (i==3) wynik=plik4;
if (i==4) wynik=plik5;
if (i==5) wynik=plik6;
if (i==6) wynik=plik7;
if (i==7) wynik=plik8;
if (i==8) wynik=plik9;
if (i==9) wynik=plik10;
if (i==10) wynik=plik11;
if (s0[i]<0)
{
while (s[i]>((1+ZAKRES)*s0[i]))
{
Fn=f(s,t,v,d,M,LP,B,N,PI);
fprintf(wynik,"%10.5f%10.5f\n",s[i]/param[i],Fn);
s[i]=s[i]+KROK*s0[i];
}
}
else
{
while (s[i]<((1+ZAKRES)*s0[i]))
{
Fn=f(s,t,v,d,M,LP,B,N,PI);
fprintf(wynik,"%10.5f%10.5f\n",s[i]/param[i],Fn);
s[i]=s[i]+KROK*s0[i];
}
}
}
fclose(plik1);
fclose(plik2);
fclose(plik3);
fclose(plik4);
fclose(plik5);
fclose(plik6);
fclose(plik7);
fclose(plik8);
fclose(plik9);
fclose(plik10);
fclose(plik11);
return 0;
}
/*funkcja obliczajaca E w czasie t metoda Newtona-Raphsona IV rzedu z dokladnoscia b*/
float NR(float P, float tau, float t, float e, float b, float pi)
{
float MM,EE,f0,f1,f2,f3,d1,d2,d3;
MM=2*pi*(t-tau)/P;
EE=MM;
f0=EE-e*sin(EE)-MM;
f1=1-e*cos(EE);
f2=e*sin(EE);
f3=e*cos(EE);
d1=-f0/f1;
d2=-f0/(f1+d1*f2/2);
d3=-f0/(f1+d2*f2/2+d2*d2*f3/6);
while ((d3>b) || (d3<-b))
{
EE=EE+d3;
f0=EE-e*sin(EE)-MM;
f1=1-e*cos(EE);
f2=e*sin(EE);
f3=e*cos(EE);
d1=-f0/f1;
d2=-f0/(f1+d1*f2/2);
d3=-f0/(f1+d2*f2/2+d2*d2*f3/6);
}
return EE;
}
/*funkcja obliczajaca funkcje dopasowania chi^2*/
float f(float s[], float t[], float v[], float d[], int Npom, int Lpar, float b, int Lplan, float pi)
{
float chi,V0,K[Lplan],P[Lplan],e[Lplan],w[Lplan],tau[Lplan],beta2,pe[Lplan],niu[Lplan],EE[Lplan];
int i,j;
chi=0;
for (i=0; i<=Npom-1; ++i)
{
V0=s[0];
for (j=0; j<=Lplan-1; ++j)
{
K[j]=s[5*j+1];
P[j]=s[5*j+2];
e[j]=s[5*j+3];
w[j]=s[5*j+4];
tau[j]=s[5*j+5];
}
for (j=0; j<=Lplan-1; ++j)
{
pe[j]=sqrt((1+e[j])/(1-e[j]));
EE[j]=NR(P[j],tau[j],t[i],e[j],b,PI);
niu[j]=2*atan(pe[j]*tan(EE[j]/2));
}
beta2=v[i]-V0;
for (j=0; j<=Lplan-1; ++j)
{
beta2=beta2-K[j]*(cos(w[j]+niu[j])+e[j]*cos(w[j]));
}
beta2=(beta2/d[i])*(beta2/d[i]);
chi=chi+beta2;
}
chi=chi/(Npom-Lpar-1);
return chi;
}
--
CezaryMigaszewski - 26 Feb 2005