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);
    };