sublimiter 2011. jan. 07. 21:09 | válasz | #1344
1.215947e-04 2.432175e-04 2.000232e+00 Lol, a kapott szog ketszerese a newtoninak, pedig csak a fenysebesseg kiszamolasahoz hasznaltam az egyik Schwarzschild egyenletet, mit nemreg linkeltem. A progi a fenti egyenlettel szamolja a feny iranyvaltozasat, tehat fenytorest szammol. dpx2 = dpx *c2*c2/(c1*c1); Ennyi. Vagyis nem egeszen. Tomeggel rendelkezo testre fenyszeru belso mozgasokat kell szamolni. Ez hazi feladat. Az eredmeny tobb, mint meglepo. Nem en fogom leirni a megoldast. #include <stdio.h> #include <stdlib.h> #include <math.h> typedef long double skalar; struct vektor { skalar x,y,z; vektor() {x=0;y=0;z=0;}; vektor(int x_,int y_,int z_) {x=x_;y=y_;z=z_;}; vektor(skalar x_,skalar y_,skalar z_) {x=x_;y=y_;z=z_;}; vektor operator + (vektor v) {vektor v_;v_.x=x+v.x;v_.y=y+v.y;v_.z=z+v.z;return v_;}; vektor operator - (vektor v) {vektor v_;v_.x=x-v.x;v_.y=y-v.y;v_.z=z-v.z;return v_;}; vektor operator * (skalar s) {vektor v_;v_.x=x*s;v_.y=y*s;v_.z=z*s;return v_;}; vektor operator / (skalar s) {vektor v_;v_.x=x/s;v_.y=y/s;v_.z=z/s;return v_;}; }; skalar skalar_szorzat(vektor v1,vektor v2) { return (v1.x*v2.x + v1.y*v2.y + v1.z*v2.z);} skalar hossz(vektor v1) {return sqrtl(skalar_szorzat(v1,v1));} vektor normalt(vektor v1) { return (v1/hossz(v1));} vektor vektor_szorzat(vektor va,vektor vb) { vektor v_; v_.x=(va.y*vb.z)-(vb.y*va.z); v_.y=(va.z*vb.x)-(vb.z*va.x); v_.z=(va.x*vb.y)-(vb.x*va.y); return v_; }; skalar g,c,m,rs,dr,rx,ry,dt; vektor p,dp,ddp; int lepes; void newton_gravity() { skalar r; while(p.x<rx) { p = p + dp*dt ; ddp = normalt(p); r = hossz(p); ddp = normalt(p)*(-g*m/(r*r)); dp = dp + ddp*dt; p = p + ddp*dt*dt/2; } } void refract_gravity() { skalar r,dpx,dpy,dpx2,dpy2,c1,c2; vektor nx,ny; while(p.x<rx) { r = hossz(p); c1 = c*(1-g*m/(2*c*c*r)) / pow(1+g*m/(2*c*c*r),3.0); p = p + dp*dt ; r = hossz(p); c2 = c*(1-g*m/(2*c*c*r)) / pow(1+g*m/(2*c*c*r),3.0);//uj fenysebesseg ny = normalt(p);//beesesi meroleges nx.x = ny.y; nx.y = -ny.x; nx.z = ny.z; dpx = skalar_szorzat(nx,dp); dpy = skalar_szorzat(ny,dp); dpx2 = dpx *c2*c2/(c1*c1); dpy2 = sqrtl(c2*c2 - dpx2*dpx2); if(dpy<0) dpy2 = -dpy2; dp = (nx*dpx2 + ny*dpy2); } } void alapallapot() { g = (skalar)6.67428e-11; c = (skalar)2.99792458e8; m = (skalar)1.9891e30; rs = (skalar)2*m*g/(c*c); ry = 1392000e3;//sun R rx = 150e9; lepes = 2000000; dt = (rx*2/c)/lepes; p = vektor(-rx,ry,(skalar)0.0); dp = vektor((skalar)c,(skalar)0.0,(skalar)0.0); } int main() { alapallapot(); newton_gravity(); skalar fi = -atanl((p.y-ry)/rx)*180/M_PI; printf("%Le \n",fi); alapallapot(); refract_gravity(); skalar fi2 = -atanl((p.y-ry)/rx)*180/M_PI; printf("%Le \n",fi2); printf("%Le \n",fi2/fi); return 0; }