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 ((modul>EPSILON) && (((Fp-Fn)>EPSILON2) || ((Fp-Fn)<-EPSILON2))) { 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
This topic: Main
>
TWikiUsers
>
CezaryMigaszewski
>
MetodaGradientowa
Topic revision: revision 2 (raw view)
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