You are here:
Foswiki
>
Main Web
>
TWikiUsers
>
CezaryMigaszewski
>
MetodaGradientowa
(15 Mar 2005,
BoudRoukema
)
(raw view)
E
dit
A
ttach
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. ----------- <verbatim> /*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) */ } } } </verbatim> -- Main.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: <verbatim> nazwa="HD114729.dat"; plik=fopen(nazwa,"r"); if(dane==NULL) { printf("Niestety, nie widze plik HD114729.dat . Prosze go przygotowac.\n"); exit (0); }; </verbatim> Zrobiłem soft link _ln -s HD114729_modelowa.dat HD114729.dat_ z plik od SymulowaneDane i teraz robi coś...
E
dit
|
A
ttach
|
P
rint version
|
H
istory
: r5
<
r4
<
r3
<
r2
|
B
acklinks
|
V
iew topic
|
Edit
w
iki text
|
M
ore topic actions
Topic revision: r5 - 15 Mar 2005,
BoudRoukema
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