Warning: Undefined array key "HTTP_REFERER" in /var/www/vhosts/web61970.ssd-space.de/phi.pf-control.de/uni/download.php on line 13
#include #include #include #include #include #include #include #include #include using namespace std; /*Dieses Programm simuliert Teilchenkollisionen über Runge Kutta 4. Ordnung. Dabei werden mehrere Teilchen mit definierter Masse zunächst in einem Gitter angeordnet und in den Grundzustand fallen gelassen. Die so entstandene Datei kann wieder eingelesen werden, um die Teilchen weiter abklingen zu lassen oder um die Teilchen 2 Mal einzufügen und aufeinander zu schießen. Stoßparameter, Anfangsimpulse und Temperatur sind dabei anpassbar. Ausgegeben wird der Endzustand der Teilchen bzw. Wahlweise eine Animation der Teilchenbewegung mit aktueller Uhrzeit als Datum. Benötigt wird gnuplot und microsoft visual studio c++ 2010 express o.Ä. */ //Neue Blockgröße int anzahl=4*4*4*2; //2 4*4*4-Blöcke aufeinander int kantenlang=4; //Kantenlänge der Blöcke für neue Blöcke. pow(4*4*4,1/3.0) ist hier zu ungenau... int tstart=0; //nur proforma, braucht nicht geändert werden. Erweiterbar: falls zeitabhängigkeit in Rungekutta erwünscht! long schritte=10000;//=>0.005ps double tende=0.05,//ns inabst1=4.3, //Rasterabstand Angström Block1 (bei neuen Dateien) temp=300; //Temperatur (K) über Monte-Carlo (bis 10000K). Restliche Parameter //restliche Parameter: Stoßparameter, Anfangsimpuls, Datei einlesen oder neue erzeugen siehe in letzter Routine "main"! double const masse=0.0028;//28u, alle Teilchen gleiche Masse, 1u =1.036426866d-4 eV ps^2/A^2 //Parameter für Potential const double a=5, b=0.27, sigma1=1.3, sigma2=3.0, gamma=10, //Reibung für Grundzustand kb=0.000086173324, //Boltzmann pi=3.14159265359, c=299792458; //Zwischenvariabeln für Drehimpuls/Drehtransfos. vectormitte(3,0),drehimp(3,0),vektor(3,0),winkelg(3,0),l(anzahl*6,0),r0(anzahl*6,0),r01(anzahl*6,0),schwerp(3,0),winkelg1(3,0); vector>traeg(3,vector(3,0)); vector>temp1(3,vector(4,0)),temp2(3,vector(4,0)); double rw,phiw,thetaw; //Berechnet Drehimpuls und Winkelgeschwindigkeit (über Gauß). //Dadurch lässt sich der Grundzustand so drehen, dass die Abweichung leichter ablesbar wird. void Drehimpuls(vector &y){ int i,j; for(i=0;i &y){ double energie; int i=0; for(i=0;i>drehz(3,vector(3,0)); drehz[0][0]=cos(-phiw); drehz[0][1]=sin(-phiw); drehz[0][2]=0; drehz[1][0]=-sin(-phiw); drehz[1][1]=cos(-phiw); drehz[1][2]=0; drehz[2][0]=0; drehz[2][1]=0; drehz[2][2]=1; vector>drehy(3,vector(3,0)); drehy[0][0]=cos(-thetaw); drehy[0][1]=0; drehy[0][2]=-sin(-thetaw); drehy[1][0]=0; drehy[1][1]=1; drehy[1][2]=0; drehy[2][0]=sin(-thetaw); drehy[2][1]=0; drehy[2][2]=cos(-thetaw); for(int i=0;i<3;i++){ for(int j=0;j<3;j++){ winkelg1[i]+=drehz[i][j]*winkelg[j]; }} for(int i=0;i<3;i++) cout<<"winkelg1"<(3,0); for(int i=0;i<3;i++){ for(int j=0;j<3;j++){ winkelg[i]+=drehy[i][j]*winkelg1[j]; } } for(int i=0;i<3;i++) cout<<"Winkelgeschwindigkeit transformiert"<(anzahl*6,0); for(int i=0;i(anzahl*6,0); for(int i=0;i(anzahl*6,0); for(int i=0;i &y){ //vertmax: Funktion für Maxima der Verteilung (über Mathematica) double vertmax = 181.69723638904276- 49036.1596250/(temp*temp) + 12007.293973803038/temp - 0.2326018089098843 *temp + 0.00020143492718019792 *temp*temp - 9.889627385975435*pow(10.0,-8)*pow(temp,3)+ 2.8465004602038517*pow(10.0,-11) *pow(temp,4) - 4.898019128141887*pow(10.0,-15) *pow(temp,5) + 4.950284303906514*pow(10.0,-19) *pow(temp,6) - 2.7071514724946854*pow(10.0,-23) *pow(temp,7) + 6.176463667836648*pow(10.0,-28) *pow(temp,8); double betragmax = 0.0017642832638790238+ 0.0012577681693344063*sqrt(temp)+0.002; //Verteilung hätte hier 10^-5 double betrag,test,vert,theta,phi; srand(time(NULL)); int j=0; for(int i=0;ipow(2*kb*temp*masse*pi,-3/2) *4*pi*betrag*betrag*exp(-test))){ break; } } if(j==2000000){cout<<"Mehr als 2d6 Überschreitungen!";} y[6*i+3] = betrag*sin(theta)*cos(phi); y[6*i+4] = betrag*sin(theta)*sin(phi); y[6*i+5] = betrag*cos(theta); } } //Runge Kutta Algorithmus 4. Dimension übernommen aus numerical recipies, jedoch angepasst auf double und vector. void rk4(vector y, vector dydx, double n, double x, double h, vector &yout, void (*derivs)(double, vector, vector&)){ //derives sit Gleichungssystem int i; double xh,hh,h6; vector dym(n,1); vector dyt(n,1); vector yt(n,1); hh=h*0.5; h6=h/6.0; xh=x+hh; for (i=0;i vstart, double nvar, double x1, double x2, double nstep,vector > &yy,vector &xx, void (*derivs)(double, vector,vector&)) { int i,k; double x,h,betr; vector v(nvar,1); //v ist im jeweiligen Schritt der aktuelle Satz an Orts/Impulskoordinaten vector vout(nvar,1); vector dv(nvar,1); for (i=0;i(3,0); for(int j=0;j(3,0); for(int j=0;jabs(gesenergie*1.1)) cout<<"\nEnergieerhaltung verletzt: Schrittweite verringern: "< y, vector &dxdy){ double betr=0; double temp=0; for(int i=0;i T fromString(const std::string& s) { std::istringstream stream (s); T t; stream >> t; return t; } //aktuelle Uhrzeit für Dateiausgabe von Animationen char* currentDateTime() { time_t now = time(0); struct tm tstruct; char buf[80]; tstruct = *localtime(&now); strftime(buf, sizeof(buf), "%Y_%m_%d_%H_%M", &tstruct); return buf; } int main(){ cout<<"hello!"<t(schritte+1,0); //Zeitschritte, hier nicht benutzt, aber möglich für Zeitabhängigkeit ifstream readfile; vectory(anzahl*6,0); //die Teilchen! vectory1(anzahl/2*6,0); //doppeltes einlesen vom Grundzustand sodass Teilchen gegeneinander geschossen werden können //x1 y1 z1 px1 py1 pz1 , x0 y0 z0 px0 py0 pz0 //Dieser Bereich belegt ein Gitter der Kantenlänge "kantenlang" im Abstand von inabst1 mit Teilchen. /* for(int i=0;i(line); } readfile.close(); //Wie vorheriger Bereich, Unterschied: hier werden Daten direkt eingelesen //Reichen die bisherigen Schritte/Zeit nicht aus, kann hier alles erneut geladen und weiter gerechnet werden /*readfile.open("endzustand.temp"); std::string line; int z=0; for(int i=0;i(line); } readfile.close();*/ //Test für das erstellen von 2 Gittern, die man aufeinander schießen könnte //Achtung: keine Grundzustände hier! /* for(int i=0;i abstand if(i%6==0){y[i]+=15;}//x-achse => stoßparam } anzahl*=2; //tatsächliche Anzahl wieder herstellen vector > yy(anzahl*6,vector(schritte+1,0)); //Ergebnis: y[index][schritt] vectordxdy(anzahl*6,0); vectortout(anzahl*6,0); //Runge Kutta-Schnittstelle rkdumb(y,anzahl*6,tstart,tende,schritte,yy,t,myderiv); wofstream file; file.open("endzustand.temp"); //Ausgabe-Datei des letzten Schrittes! jede Zeile eine Zahl: x y z px py pz x ... for(int i=0;i