*
oh ja indem ich folgendes tue
for (i=start; i<particles;i++)
{
//r(t+delta/2)
F2_r[i][0] = (pold[i][0]+ delta/2 * F1_p[i][0]) + gam * (rold[i][1]+ delta/2 * F1_r[i][1]) ;
F2_r[i][1] = pold[i][1] + delta/2 * F1_p[i][1];
F2_r[i][2] = pold[i][2] + delta/2 * F1_p[i][2];
//temporaere Teilchenkraft
temp_F[i][0] = F[i][0]+ delta/2 * tempKraft[i][0];
temp_F[i][1] = F[i][1]+ delta/2 * tempKraft[i][1];
temp_F[i][2] = F[i][2]+ delta/2 * tempKraft[i][2];
//temporaere Impulse
temp_p[i][0] = pold[i][0] + delta/2 * F1_p[i][0];
temp_p[i][1] = pold[i][1] + delta/2 * F1_p[i][1];
temp_p[i][2] = pold[i][2] + delta/2 * F1_p[i][2];
temp_omz = omzold + delta/2 * F1_omz;
//p(t+delta/2)
F2_p[i][0] = temp_F[i][0] - gam * temp_p[i][1]
- zet(temp_F,temp_p) * temp_p[i][0]
- nustir/temp * Fdotp(temp_F,temp_p) * temp_omz * temp_p[i][1];
F2_p[i][1] = temp_F[i][1] - zet(temp_F,temp_p) * temp_p[i][1]
+ nustir/temp * Fdotp(temp_F,temp_p) * temp_omz * temp_p[i][0];
F2_p[i][2] = temp_F[i][2] - zet(temp_F,temp_p) * temp_p[i][2] ;
//omz(t+delta/2)
F2_omz = omz_dt(temp_F,temp_p);
}
... mache ich tatsächlich einen riesen Fehler ... die Größen zet, Fdotp und omz_dt hängen natürlich von allen temp_F und temp_p ab. Ich benötige zuerst das gesamte Array von temp_F und temp_p um F2_p berechnen zu können. Alles mit 0 zu initialiesieren hat deshalb geholfen, weil die 0 eine gutartige Zahl ist, die die Beiträge, die ich falsch berechne einfach löscht. Ohne Initialisierung mit 0 stand in diesen Arrays zu Beginn nur Datenmüll, was das sofortige Auftreten von nan erklärt.
Das ist alles ein Copy and paste Fehler denn bei der Berechnung von F1_p müsste es wie oben gut gehen, da ich dort nur rold, pold und F benötige, die als vollständig gefüllte Arrays vorliegen.
Aber dass der gesamte Code murks ist würd ich jetzt nicht behaupten. Meine Änderung sollte die folgende sein:
//Kraft am Ort F1_r ausrechnen
Gesamtkraft(F1_r,tempKraft); //tempKraft liegt als vollständiger Satz vor :)
//Schleife über alle Teilchen
for (i=start; i<particles;i++)
{
//alle F1_p, und F1_r, sowie pold liegen vollständig vor
F2_r[i][0] = (pold[i][0]+ delta/2 * F1_p[i][0]) + gam * (rold[i][1]+ delta/2 * F1_r[i][1]) ;
F2_r[i][1] = pold[i][1] + delta/2 * F1_p[i][1];
F2_r[i][2] = pold[i][2] + delta/2 * F1_p[i][2];
//auch alle tempKraft habe ich
//temporaere Teilchenkraft
temp_F[i][0] = F[i][0]+ delta/2 * tempKraft[i][0];
temp_F[i][1] = F[i][1]+ delta/2 * tempKraft[i][1];
temp_F[i][2] = F[i][2]+ delta/2 * tempKraft[i][2];
//temporaere Impulse
temp_p[i][0] = pold[i][0] + delta/2 * F1_p[i][0];
temp_p[i][1] = pold[i][1] + delta/2 * F1_p[i][1];
temp_p[i][2] = pold[i][2] + delta/2 * F1_p[i][2];
//F1_omz gibt es aus dem vorherigen Schritt, ist auch kein Array
temp_omz = omzold + delta/2 * F1_omz;
} //hier endet die Schleife über alle Teilchen
//hier beginnt wieder Schleife über alle Teilchen
// ich besitze temp_F und temp_p vollständig, da schleife oben wieder geschlossen
//ich arbeite also mit den physikalisch richtigen Werten
for (i=start; i<particles;i++)
{
F2_p[i][0] = temp_F[i][0] - gam * temp_p[i][1]
- zet(temp_F,temp_p) * temp_p[i][0]
- nustir/temp * Fdotp(temp_F,temp_p) * temp_omz * temp_p[i][1];
F2_p[i][1] = temp_F[i][1] - zet(temp_F,temp_p) * temp_p[i][1]
+ nustir/temp * Fdotp(temp_F,temp_p) * temp_omz * temp_p[i][0];
F2_p[i][2] = temp_F[i][2] - zet(temp_F,temp_p) * temp_p[i][2] ;
//omz(t+delta/2)
F2_omz = omz_dt(temp_F,temp_p);
}
Danke für den Tip ... ich wär da so schnell nicht drauf gekommen ... mal schauen ob's dann geht ...
jesus (schämt sich ein wenig ... )
edit: temp_p[0] und so brauche ich eigentlich nicht, ich setze auf das Teilchen 0 meinen Koordinatenursprung. In meinem mitbewegte Koordinatensystem dürfen keine Kräfte auf Teilchen 0 wirken (sonst wär das Koordinatensystem nicht fest auf 0 ) deshalb durchlaufe ich alle Schleifen mit dem Startwert 1 (#define start 1). Ich hab das so gewählt, weil ich später auch mal die Dynamik des Systems gegenüber einem festen Beobachtungsort sehen wollte, dann muss ich nur start auf 0 setzen und schon wird alles für Teilchen 0 mitberechnet. solnage ich das nciht tue bleiben alle Parameter für Teilchen 0 wie sie anfangs initialisiert wurden