You are here:
Foswiki
>
Main Web
>
TWikiUsers
>
CezaryMigaszewski
>
FunkcjaDopasowania
(revision 1) (raw view)
Edit
Attach
Omowienie programu wkrotce --------------- <verbatim> /*chi2.c*/ #include<stdio.h> #include<math.h> #include<stdlib.h> #define DZIEN 86400 /*liczba sekund w 1 dniu*/ #define N 2 /*liczba planet*/ #define LP 11 /*liczba parametrow*/ #define M 38 /*liczba pomiarow*/ #define B 0.0001 /*dokladnosc rozwiazania r-nia Keplera*/ #define KROK 0.001 /*podzial siatki argumentow*/ #define ZAKRES 0.05 /*zakres zmian parametrow*/ #define PI 3.141592653589793115997963468544185161590576172 /*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); main() { int i,j; float Fn; float V0,K[N],P[N],e[N],w[N],tau[N]; float s[LP],s0[LP],t[M],v[M],d[M]; FILE *plik, *plik1, *plik2, *plik3, *plik4, *plik5, *plik6, *plik7, *plik8; FILE *plik9, *plik10, *plik11, *wynik; char *nazwa; float param[LP]={1,1,DZIEN,1,PI/180,DZIEN,1,DZIEN,1,PI/180,DZIEN}; /*wczytanie danych*/ nazwa="HD114729.dat"; plik=fopen(nazwa,"r"); for (i=0; i<=M-1; ++i) { fscanf(plik,"%f%f%f",&t[i],&v[i],&d[i]); t[i]=t[i]*DZIEN; } fclose(plik); /*parametry poczatkowe*/ V0=-5.32061; K[0]=15.29573; P[0]=1088.65051*DZIEN; e[0]=0.205788; w[0]=75.49266*PI/180; tau[0]=509.81573*DZIEN; K[1]=5.69828; P[1]=10.52720*DZIEN; e[1]=0.016552; w[1]=38.30322*PI/180; tau[1]=100.91705*DZIEN; nazwa="HD114729_chi2_V0_1.dat"; plik1=fopen(nazwa,"w"); nazwa="HD114729_chi2_K1_1.dat"; plik2=fopen(nazwa,"w"); nazwa="HD114729_chi2_P1_1.dat"; plik3=fopen(nazwa,"w"); nazwa="HD114729_chi2_e1_1.dat"; plik4=fopen(nazwa,"w"); nazwa="HD114729_chi2_w1_1.dat"; plik5=fopen(nazwa,"w"); nazwa="HD114729_chi2_tau1_1.dat"; plik6=fopen(nazwa,"w"); nazwa="HD114729_chi2_K2_1.dat"; plik7=fopen(nazwa,"w"); nazwa="HD114729_chi2_P2_1.dat"; plik8=fopen(nazwa,"w"); nazwa="HD114729_chi2_e2_1.dat"; plik9=fopen(nazwa,"w"); nazwa="HD114729_chi2_w2_1.dat"; plik10=fopen(nazwa,"w"); nazwa="HD114729_chi2_tau2_1.dat"; plik11=fopen(nazwa,"w"); /* glowna czesc programu*/ s0[0]=V0; for (j=0; j<=N-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 (i=0; i<=LP-1; ++i) { for (j=0; j<=LP-1; ++j) { s[j]=s0[j]; } s[i]=(1-ZAKRES)*s0[i]; if (i==0) wynik=plik1; if (i==1) wynik=plik2; if (i==2) wynik=plik3; if (i==3) wynik=plik4; if (i==4) wynik=plik5; if (i==5) wynik=plik6; if (i==6) wynik=plik7; if (i==7) wynik=plik8; if (i==8) wynik=plik9; if (i==9) wynik=plik10; if (i==10) wynik=plik11; if (s0[i]<0) { while (s[i]>((1+ZAKRES)*s0[i])) { Fn=f(s,t,v,d,M,LP,B,N,PI); fprintf(wynik,"%10.5f%10.5f\n",s[i]/param[i],Fn); s[i]=s[i]+KROK*s0[i]; } } else { while (s[i]<((1+ZAKRES)*s0[i])) { Fn=f(s,t,v,d,M,LP,B,N,PI); fprintf(wynik,"%10.5f%10.5f\n",s[i]/param[i],Fn); s[i]=s[i]+KROK*s0[i]; } } } fclose(plik1); fclose(plik2); fclose(plik3); fclose(plik4); fclose(plik5); fclose(plik6); fclose(plik7); fclose(plik8); fclose(plik9); fclose(plik10); fclose(plik11); return 0; } /*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; } </verbatim> -- Main.CezaryMigaszewski - 26 Feb 2005
Edit
|
Attach
|
P
rint version
|
H
istory
:
r2
<
r1
|
B
acklinks
|
V
iew topic
|
Edit WikiText
|
More topic actions...
Topic revision: r1 - 26 Feb 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