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
This topic: Main
>
TWikiUsers
>
CezaryMigaszewski
>
FunkcjaDopasowania
Topic revision: revision 1 (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