You are here:
Foswiki
>
Main Web
>
TWikiUsers
>
CezaryMigaszewski
>
MonteCarlo
(revision 2) (raw view)
Edit
Attach
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. -------------------- <verbatim> /*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; } </verbatim> -- Main.CezaryMigaszewski - 26 Feb 2005
Edit
|
Attach
|
P
rint version
|
H
istory
:
r3
<
r2
<
r1
|
B
acklinks
|
V
iew topic
|
Edit WikiText
|
More topic actions...
Topic revision: r2 - 05 Mar 2005,
CezaryMigaszewski
Main
Log In
or
Register
Toolbox
Create New Topic
Index
Search
Changes
Notifications
RSS Feed
Statistics
Preferences
Users
Groups
Webs
Cosmo
Main
Sandbox
System
English
Français
Polski
Copyright © CC-BY-SA by the contributing authors. All material on this collaboration platform is copyrighted under CC-BY-SA by the contributing authors unless otherwise noted.
Ideas, requests, problems regarding Foswiki?
Send feedback