Metoda Monte Carlo jest losowa metoda numeryczna. Program znajduje w sposob losowy najlepsze parametry planet. Metoda ta nie jest zbyt wydajna. Polega na tym, ze losujemy zespol parametrow i sprawdzamy jaka wartosc ma funkcja dopasowania, po czym znowu losujemy parametry i sprawdzamy. Jezeli dopasowanie jest lepsze zachowujemy je w pamieci. Czynimy tak do czasu znalezienia odpowiednio dobrego dopasowania lub po wykonaniu odpowiedniej liczby krokow. W programie nalezy podac zakrs parametrow (granice strzalu) oraz nazwe pliku wyjsciowego.
/*MonteCarlo*/
#include<stdio.h>
#include<math.h>
#include<stdlib.h>
#define LKMAX 100000000
#define LP 11
#define M 142
#define N 2
#define KMAX 150
#define DZIEN 86400
#define PMAX 1000
#define VMAX 100
#define PI 3.1415926
#define B 0.0001
float kepler(float okres,float tau, float t, float e, float b);
main()
{
FILE *dane;
FILE *wynik;
char *plik_dane;
char *plik_wynik;
float t[M],v[M],d[M];
int i,j,lk;
float pom,v0,chi,chi2,sum,EE,niu;
float K[N],P[N],e[N],w[N],tau[N];
int pomocnicza;
plik_dane="hr7.dat";
dane=fopen(plik_dane,"r");
for (j=0; j<=M-1; ++j)
{
fscanf(dane,"%f%f%f\n",&t[j],&v[j],&d[j]);
t[j]=t[j]*DZIEN;
}
fclose(dane);
/* losowanie parametrow poczatkowych */
pom=rand();
v0=VMAX*(pom/RAND_MAX-0.5);
for (i=0; i<N; ++i)
{
pom=rand();
K[i]=KMAX*pom/RAND_MAX;
pom=rand();
e[i]=pom/RAND_MAX;
pom=rand();
w[i]=2*PI*pom/RAND_MAX;
pom=rand();
P[i]=PMAX*DZIEN*pom/RAND_MAX;
pom=rand();
tau[i]=P[i]*pom/RAND_MAX+t[0];
}
/* obliczamy chi^2/(M-LP-1) */
chi=0.0;
for (j=0; j<=M-1; ++j)
{
sum=0.0;
for (i=0; i<N; ++i)
{
EE=kepler(P[i],tau[i],t[j],e[i],B);
niu=2*atan(sqrt((1+e[i])/(1-e[i]))*tan(EE/2));;
if (niu<0) niu=niu+2*PI;
sum=sum+K[i]*(cos(w[i]+niu)+e[i]*cos(w[i]));
}
chi=chi+((v[j]-v0-sum)/d[j])*((v[j]-v0-sum)/d[j]);
}
chi=chi/(M-LP-1);
/* zapisujemy parametry oraz chi do pliku */
plik_wynik="fit9.dat";
wynik=fopen(plik_wynik,"a");
for (i=0; i<N; ++i) {
fprintf(wynik,"%10.4f%10.4f%10.4f%15.4f%15.4f",K[i],e[i],w[i]*180/PI,P[i]/DZIEN,tau[i]/DZIEN);
}
fprintf(wynik,"%10.4f%10.4f\n",v0,chi);
/* losujemy parametry i sprawdzamy, czy daja lepsze dopasowanie. Jesli tak, to zapisaujemy dane do pliku */
pomocnicza=0;
for (lk=1; lk<=LKMAX; ++lk)
{
pom=rand();
v0=VMAX*(pom/RAND_MAX-0.5);
for (i=0; i<N; ++i)
{
pom=rand();
K[i]=KMAX*pom/RAND_MAX;
pom=rand();
e[i]=pom/RAND_MAX;
pom=rand();
w[i]=2*PI*pom/RAND_MAX;
pom=rand();
P[i]=PMAX*DZIEN*pom/RAND_MAX;
pom=rand();
tau[i]=P[i]*pom/RAND_MAX+t[0];
}
chi2=0.0;
for (j=0; j<=M-1; ++j)
{
sum=0.0;
for (i=0; i<N; ++i)
{
EE=kepler(P[i],tau[i],t[j],e[i],B);
niu=2*atan(EE/2);
if (niu<0) niu=niu+2*PI;
sum=sum+K[i]*(cos(w[i]+niu)+e[i]*cos(w[i]));
}
chi2=chi2+((v[j]-v0-sum)/d[j])*((v[j]-v0-sum)/d[j]);
}
chi2=chi2/(M-LP-1);
if (chi2<chi)
{
chi=chi2;
for (i=0; i<N; ++i) {
fprintf(wynik,"%10.4f%10.4f%10.4f%15.4f%15.4f",K[i],e[i],w[i]*180/PI,P[i]/DZIEN,tau[i]/DZIEN); }
fprintf(wynik,"%10.4f%10.4f\n",v0,chi);
}
pomocnicza=pomocnicza+1;
if (pomocnicza==100000)
{
fclose(wynik);
wynik=fopen(plik_wynik,"a");
pomocnicza=0;
}
}
fclose(wynik);
return 0;
}
/* funkcja rozwiazujaca rownanie Keplera */
float kepler(float okres,float tau, float t, float e, float b)
{
float MM,EE,f0,f1,f2,f3,d1,d2,d3,pi;
pi=3.1415926;
MM=2*pi*(t-tau)/okres;
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;
}
--
CezaryMigaszewski - 26 Feb 2005
Brakuje plik wejść
hr7.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:
plik_dane="hr7.dat";
dane=fopen(plik_dane,"r");
if(dane==NULL)
{
printf("Niestety, nie widze plik hr7.dat . Prosze go przygotowac.\n");
exit (0);
};