Program sluzy do szukania minimum funkcji dopasowania w przestrzeni parametrow opisujacy orbity planet. Wykorzystuje metode gradientowa opracowana przez autora tych slow. W programie nalezy podac parametry startowe oraz nazwe pliku z danymi pomiarowymi predkosci radialnej gwiazdy, oraz plikow wyjsciowych. Mozliwe, ze zamieszcze opis samej metody, ale nie wiem dokladnie kiedy. Prgram startuje z parametrow zadanych w programie i poruszajac sie w kierunku wyznaczanym orzez tak zwany "gradient wzgledny" szuka minimum funkcji dopasowania. Minimum globalne tej funkcji oznacza najlepsze parametry orbitalne. Jednak funkcja dopasowania jest skomplikowana, wiec wynik koncowy zalezy od miejsca w przestrzeni parametrow, z ktorego wystartujemy.
/*GRADIENT2.c*/
#include<stdio.h>
#include<math.h>
#define DZIEN 86400 /*liczba sekund w 1 dniu*/
#define HP 0.01 /*poczatkowy 'krok gradientowy'*/
#define ALFA 1.001 /*'wspolczynnik miniecia'*/
#define MP 2 /*liczba planet*/
#define LP 11 /*liczba parametrow*/
#define EPSILON 0.1 /*konczaca wartosc gradientu*/
#define EPSILON2 1E-20 /* Fp-Fn */
#define N 38 /*liczba pomiarow*/
#define B 0.0001 /*dokladnosc rozwiazania r-nia Keplera*/
#define PI 3.1415926536 /*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);
void g(float s[],float t[],float v[],float d[], int Npom, float b, int Lplan, float grad[], float pi);
main()
{
int i,j,numerKroku;
float t[N],v[N],d[N],s[LP],s0[LP],gradP[LP],gradN[LP],Fmin,smin[LP],Fn,Fp,modul,modul0,h[LP];
float V0,K[MP],P[MP],e[MP],w[MP],tau[MP];
FILE *plik, *wynik, *gradienty;
char *nazwa, *nazwa_w, *nazwa_g;
/*wczytanie danych*/
nazwa="HD114729.dat";
plik=fopen(nazwa,"r");
for (i=0; i<=N-1; ++i)
{
fscanf(plik,"%f%f%f",&t[i],&v[i],&d[i]);
t[i]=t[i]*DZIEN;
}
fclose(plik);
/*parametry poczatkowe*/
V0=-5.5;
K[0]=15;
P[0]=1080*DZIEN;
e[0]=0.19;
w[0]=81*PI/180;
tau[0]=539*DZIEN;
K[1]=6;
P[1]=10.5262*DZIEN;
e[1]=0.016;
w[1]=53*PI/180;
tau[1]=101.4*DZIEN;
/* glowna czesc programu*/
s0[0]=V0;
for (j=0; j<=MP-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 (j=0; j<=LP-1; ++j)
{
s[j]=s0[j];
}
nazwa_w="HD114729_GRADIENT_1_1.dat";
nazwa_g="HD114729_GRADIENT_2_1.dat";
wynik=fopen(nazwa_w,"w");
gradienty=fopen(nazwa_g,"w");
numerKroku=0;
Fn=f(s,t,v,d,N,LP,B,MP,PI);
V0=s[0];
for (j=0; j<=MP-1; ++j)
{
K[j]=s[5*j+1];
P[j]=s[5*j+2]/DZIEN;
e[j]=s[5*j+3];
w[j]=s[5*j+4]*180/PI;
tau[j]=s[5*j+5]/DZIEN;
}
for (j=0; j<=LP-1; ++j)
{
smin[j]=s[j];
}
Fmin=Fn;
g(s,t,v,d,N,B,MP,gradP,PI);
modul=0;
for (i=0; i<=LP-1; ++i)
{
modul=modul+atan(s0[i]*gradP[i])*atan(s0[i]*gradP[i]);
}
modul=sqrt(modul);
modul0=modul;
for (i=0; i<=LP-1; ++i)
{
h[i]=HP*Fn*modul0/sqrt((s0[i]*gradP[i])*(s0[i]*gradP[i]));
} /* poczatkowy "krok gradientowy" */
fprintf(wynik,"%10d%20.15f%12.5f%12.5f",numerKroku,Fn,modul,V0);
for (j=0; j<=MP-1; ++j)
{
fprintf(wynik,"%12.5f%12.5f%12.6f%12.5f%13.5f",K[j],P[j],e[j],w[j],tau[j]);
}
fprintf(wynik,"\n");
fprintf(gradienty,"%10d",numerKroku);
for (i=0; i<=LP-1; ++i)
{
fprintf(gradienty,"%12.5f",atan(s0[i]*gradP[i]));
}
fprintf(gradienty,"\n");
for (i=0; i<=LP-1; ++i)
{
s[i]=s[i]*(1-h[i]*atan(s0[i]*gradP[i])/modul);
}
Fp=Fn;
Fn=f(s,t,v,d,N,LP,B,MP,PI);
numerKroku=1;
V0=s[0];
for (j=0; j<=MP-1; ++j)
{
K[j]=s[5*j+1];
P[j]=s[5*j+2]/DZIEN;
e[j]=s[5*j+3];
w[j]=s[5*j+4]*180/PI;
tau[j]=s[5*j+5]/DZIEN;
}
g(s,t,v,d,N,B,MP,gradN,PI);
for (i=0; i<=LP-1; ++i)
{
if (gradP[i]*gradN[i]<0)
{
h[i]=h[i]/ALFA;
}
}
for (i=0; i<=LP-1; ++i)
{
gradP[i]=gradN[i];
}
modul=0;
for (i=0; i<=LP-1; ++i)
{
modul=modul+atan(s0[i]*gradP[i])*atan(s0[i]*gradP[i]);
}
modul=sqrt(modul);
fprintf(wynik,"%10d%20.15f%12.5f%12.5f",numerKroku,Fn,modul,V0);
for (j=0; j<=MP-1; ++j)
{
fprintf(wynik,"%12.5f%12.5f%12.6f%12.5f%13.5f",K[j],P[j],e[j],w[j],tau[j]);
}
fprintf(wynik,"\n");
fprintf(gradienty,"%10d",numerKroku);
for (i=0; i<=LP-1; ++i)
{
fprintf(gradienty,"%12.5f",atan(s0[i]*gradP[i]));
}
fprintf(gradienty,"\n");
if (Fn<Fmin)
{
Fmin=Fn;
for (j=0; j<=LP-1; ++j)
{
smin[j]=s[j];
}
}
while (numerKroku<1000)
{
numerKroku=numerKroku+1;
Fp=Fn;
for (i=0; i<=LP-1; ++i)
{
s[i]=s[i]*(1-h[i]*atan(s0[i]*gradP[i])/modul);
}
Fn=f(s,t,v,d,N,LP,B,MP,PI);
V0=s[0];
for (j=0; j<=MP-1; ++j)
{
K[j]=s[5*j+1];
P[j]=s[5*j+2]/DZIEN;
e[j]=s[5*j+3];
w[j]=s[5*j+4]*180/PI;
tau[j]=s[5*j+5]/DZIEN;
}
g(s,t,v,d,N,B,MP,gradN,PI);
for (i=0; i<=LP-1; ++i)
{
if (gradP[i]*gradN[i]<0)
{
h[i]=h[i]/ALFA;
}
}
for (i=0; i<=LP-1; ++i)
{
gradP[i]=gradN[i];
}
modul=0;
for (i=0; i<=LP-1; ++i)
{
modul=modul+atan(s0[i]*gradP[i])*atan(s0[i]*gradP[i]);
}
modul=sqrt(modul);
fprintf(wynik,"%10d%20.15f%12.5f%12.5f",numerKroku,Fn,modul,V0);
for (j=0; j<=MP-1; ++j)
{
fprintf(wynik,"%12.5f%12.5f%12.6f%12.5f%13.5f",K[j],P[j],e[j],w[j],tau[j]);
}
fprintf(wynik,"\n");
fprintf(gradienty,"%10d",numerKroku);
for (i=0; i<=LP-1; ++i)
{
fprintf(gradienty,"%12.5f",atan(s0[i]*gradP[i]));
}
fprintf(gradienty,"\n");
if (Fn<Fmin)
{
Fmin=Fn;
for (j=0; j<=LP-1; ++j)
{
smin[j]=s[j];
}
}
}
fclose(wynik);
fclose(gradienty);
/*----------- wyniki ------------------*/
V0=smin[0];
for (j=0; j<=MP-1; ++j)
{
K[j]=smin[5*j+1];
P[j]=smin[5*j+2]/DZIEN;
e[j]=smin[5*j+3];
w[j]=smin[5*j+4]*180/PI;
tau[j]=smin[5*j+5]/DZIEN;
}
return 0;
}
/*koniec funkcji main*/
/*funkcje*/
/*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;
}
/*funkcja liczaca gradient funkcji chi^2 w przestrzeni parametrow s*/
void g(float s[],float t[],float v[],float d[], int Npom, float b, int Lplan, float grad[], float pi)
{
float V0,K[Lplan],P[Lplan],e[Lplan],w[Lplan],tau[Lplan],beta2,beta;
float pe[Lplan],niu[Lplan],EE[Lplan],dniudP[Lplan],dniudE[Lplan],dniudtau[Lplan],lp;
int i,j;
lp=5*Lplan+1;
for (j=0; j<=lp-1; ++j)
{
grad[j]=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));
dniudE[j]=pe[j]/((cos(EE[j]/2))*(cos(EE[j]/2))+pe[j]*pe[j]*(sin(EE[j]/2))*(sin(EE[j]/2)));
dniudP[j]=-dniudE[j]*2*pi*(t[i]-tau[j])/(P[j]*P[j]*(1-e[j]*cos(EE[j])));
dniudtau[j]=-dniudE[j]*2*pi/(P[j]*(1-e[j]*cos(EE[j])));
}
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]));
}
beta=2*beta2/(d[i]*d[i]);
grad[0]=grad[0]-beta; /* dchi/dV0 */
for (j=0; j<=Lplan-1; ++j)
{
grad[5*j+1]=grad[5*j+1]-beta*(cos(w[j]+niu[j])+e[j]*cos(w[j])); /* dchi/dK(j) */
grad[5*j+2]=grad[5*j+2]+beta*K[j]*sin(w[j]+niu[j])*dniudP[j]; /* dchi/dP(j) */
grad[5*j+3]=grad[5*j+3]-beta*K[j]*cos(w[j]); /* dchi/de(j) */
grad[5*j+4]=grad[5*j+4]+beta*K[j]*(sin(w[j]+niu[j])+e[j]*sin(w[j])); /* dchi/dw(j) */
grad[5*j+5]=grad[5*j+5]+beta*K[j]*sin(w[j]+niu[j])*dniudtau[j]; /* dchi/dtau(j) */
}
}
}
--
CezaryMigaszewski - 26 Feb 2005
Daje mnie
Naruszenie ochrony pamiêci :(.
Mozna dac inny warynek konczacy petle, np. okreslona liczba krokow (juz zrobione).
Kod zostanie jeszcze poprawiony w przyszlosci (niedalekiej)
Brakuje plik wejść
HD114729.dat - Też mogłybyś zrobić test żeby uważać użytkownik jeśli plik nie istnieje - jest to więcej ,,user-friendly" niż błęd
Naruszenie ochrony pamięci 
np:
nazwa="HD114729.dat";
plik=fopen(nazwa,"r");
if(dane==NULL)
{
printf("Niestety, nie widze plik HD114729.dat . Prosze go przygotowac.\n");
exit (0);
};
Zrobiłem soft link
ln -s HD114729_modelowa.dat HD114729.dat z plik od
SymulowaneDane i teraz robi coś...