problem mit function

[ZuK]

Wenn du keine grösse wie zb. "%20.3f" angibst schreibt sprintf drüber hinaus.
Kurt
EDIT: habs gerade ausprobiert hilft nichts schreibt trotzden drüber hinaus
mit "%20.3g" funktionierts

[*char*jesus]

... nee das hilft nix ... ich kann machen was ich will sobald ich meine Variablen lokal definiere, macht das Programm nicht mehr mit ... dann steht in allen Ergebnissen nur noch nan ...

gibt es denn eine bessere Möglichkeit zu prüfen, ob ein double nan ist als in einen String zu konvertieren ?

jesus

[Don Carsto]

Hallo,
ich hab das Dingen nun auch mal kompiliert.
Unter Windows mit Dev-C++.
Läuft einwandfrei, auch wenn ich das ganze als Funktion
ändere. Allerding schmeißt mir der Compiler jede Menge
Warnungen wegen falschen Typumwandlungen raus.
Habe nun mal durch Type-Casting und abändern einiger Variablen
von long int in long double usw. den Code so verändert, daß
keine Warnungen mehr kommen.
Läuft immer noch. Mußte mal prüfen, ob er nun auch bei dir
läuft. Wenn ja, lags daran, ansonsten mußte weitersuchen.

Hier der geänderte Code:

//Teilchen mit DipolWW unter Scherung 
//Integrator: Runge-Kutta-Verfahren 4. Ordnung 

#include <stdio.h> 
#include <stdlib.h> 
#include <math.h> 
#include <time.h> 
#include <string.h> 

// Anzahl der zu simulierenden Teilchen -1 
#define particles 3         // 3 bedeutet 2 unabhaengige Teilchen !!! eines ist im Ursprung 
#define start 1         // Teilchen 0 mit einbeziehen ? 1: nein 0: ja 
#define schreibschritt 1000000    // welche Schritte sollen geschrieben werden ? >>stepcount % (steps/ schreibschritt) 
//wenn schreibschritt==steps, dann wird jeder schritt geschrieben 
#define dateiname "ausgabe.csv" //diese Datei wird geschrieben 
#define reload 1        //alle wieviel Sekunden soll die Ausgabe erfolgen ? 

double dipst   = 5.0; 
double temp    = 2.5; 
double gam     = 1.0; 
int nugauss = 0;     //1: thermostat ein, 0: thermostat aus 
double nustir  = 1.0;     //Staerke der Ankopplung des twirlers   

long double delta = 0.000001;     // Schrittweite 
long double steps = 1000000.0;       // Anzahl Schritte 
long double stepcount = 0.;     // Schritt      

//füllt uebergebenen String mit Zeitdaten 
void timetostr(char ttostr[], long int time) 
{ 
  double std, min, s; 
  char zeichen[5]; 
//std    ((time)/3600), 
//min    (fmod((time)/60,60.)), 
//s    (time ,60.0), 
    std = ((time)/3600); 
    min = fmod((time/60),60.); 
    s   = fmod(time ,60.0); 
        if (std < 10) sprintf(zeichen," 0%.0f:",floor(std)); else sprintf(zeichen," %.0f:",floor(std)); 
    strcpy(ttostr,zeichen); 
    if (min < 10) sprintf(zeichen,"0%.0f:",floor(min) ); else sprintf(zeichen,"%.0f:",floor(min)); 
    strcat(ttostr,zeichen); 
    if (s < 10) sprintf(zeichen,"0%.0f ", floor(s) ); else sprintf(zeichen,"%.0f ", floor(s) ); 
    strcat(ttostr,zeichen); 

    //sprintf(ttostr," %.0f:%.0f:%.0f ",floor(std),floor(min), floor(s) ); 
} 

//xKraft eines Teilchens am Ort (x,y,z) 
double Fx(double x, double y, double z) 
{ 
    double WandKraft, DipKraft; 
    double ergebnis=0; 

    WandKraft = - 46656/3125 *(     (pow(2-sqrt(pow(x,2)+pow(y,2)+pow(z,2)),5) 
                                    *(-x)*1/sqrt(pow(x,2)+pow(y,2)+pow(z,2))) 

                                +     (1-sqrt(pow(x,2)+pow(y,2)+pow(z,2))) 
                                    *5*pow(2-sqrt(pow(x,2)+pow(y,2)+pow(z,2)),4) 
                                    *(-x)*1/sqrt(pow(x,2)+pow(y,2)+pow(z,2)) 

                                ); 
    DipKraft = dipst * (    -( (6*pow(y,2)*x)*(0.01*dipst+pow(sqrt(pow(x,2)+pow(y,2)+pow(z,2)),3)) 
                                    /pow(pow(x,2)+pow(y,2)+pow(z,2),2)) 

                            -((3*pow(y,2)/(pow(x,2)+pow(y,2)+pow(z,2))) 
                                *3*x *sqrt(pow(x,2)+pow(y,2)+pow(z,2))) 

                        )/ pow(0.01*dipst+pow(sqrt(pow(x,2)+pow(y,2)+pow(z,2)),3),2) ; 

    ergebnis = WandKraft+ DipKraft; 

// printf("\n\t%1.5f\t%1.5f\t%1.5f\n",x,y,z); 

    return ergebnis; 
    } 

//yKraft eines Teilchens am Ort (x,y,z) 
double Fy(double x, double y, double z) 
{ 
    double WandKraft, DipKraft; 
    double ergebnis = 0.0; 

    WandKraft = - 46656/3125 *(     (pow(2-sqrt(pow(x,2)+pow(y,2)+pow(z,2)),5) 
                                    *(-y)*1/sqrt(pow(x,2)+pow(y,2)+pow(z,2))) 

                                +     (1-sqrt(pow(x,2)+pow(y,2)+pow(z,2))) 
                                    *5 *pow(2-sqrt(pow(x,2)+pow(y,2)+pow(z,2)),4) 
                                    *(-y)*1/sqrt(pow(x,2)+pow(y,2)+pow(z,2)) 

                                ); 
    DipKraft = dipst * (     ( (6*y*(pow(x,2)+pow(y,2)+pow(z,2)))- (6*pow(y,3)) 

                                *(0.01*dipst+pow(sqrt(pow(x,2)+pow(y,2)+pow(z,2)),3)) 
                                    /pow(pow(x,2)+pow(y,2)+pow(z,2),2) 
                                ) 

                            -( 
                                (3*pow(y,2)/(pow(x,2)+pow(y,2)+pow(z,2))) 
                                *3*y*sqrt(pow(x,2)+pow(y,2)+pow(z,2)) 
                                ) 
                        )/ pow(0.01*dipst+pow(sqrt(pow(x,2)+pow(y,2)+pow(z,2)),3),2) ; 

    ergebnis = WandKraft+ DipKraft; 
//    printf("\n\t%1.5f\t%1.5f\t%1.5f\n",x,y,z); 
//    printf("dipst : %1.5f , Dipkraft : %1.5f , Wandkraft : %1.5f\n",dipst, DipKraft, WandKraft); 

    return ergebnis; 
    } 

//zKraft eines Teilchens am Ort (x,y,z) 
double Fz(double x, double y, double z) 
{ 
    double WandKraft, DipKraft; 
    double ergebnis = 0.0; 

    WandKraft = - 46656/3125 *(     (pow(2-sqrt(pow(x,2)+pow(y,2)+pow(z,2)),5) 
                                    *(-z)*1/sqrt(pow(x,2)+pow(y,2)+pow(z,2))) 

                                +     (1-sqrt(pow(x,2)+pow(y,2)+pow(z,2))) 
                                    *5 *pow(2-sqrt(pow(x,2)+pow(y,2)+pow(z,2)),4) 
                                    *(-z)*1/sqrt(pow(x,2)+pow(y,2)+pow(z,2)) 

                                ); 
    DipKraft = dipst * (     -(- (6*pow(y,2)*z) 

                                *(0.01*dipst+pow(sqrt(pow(x,2)+pow(y,2)+pow(z,2)),3)) 
                                    /pow(pow(x,2)+pow(y,2)+pow(z,2),2) 
                                ) 

                            -( 
                                (3*pow(y,2)/(pow(x,2)+pow(y,2)+pow(z,2))) 
                                *3*z*sqrt(pow(x,2)+pow(y,2)+pow(z,2)) 
                                ) 
                        )/ pow(0.01*dipst+pow(sqrt(pow(x,2)+pow(y,2)+pow(z,2)),3),2) ; 

    ergebnis = WandKraft+DipKraft; 
//printf("\n\t%1.5f\t%1.5f\t%1.5f\n",x,y,z); 

    return ergebnis; 
    } 

//Gesamtkraft ausrechnen 
int Gesamtkraft(double r[][3] , double F[][3]) 
{ 
int i,j; 
double xwert =0.0, ywert=0.0, zwert=0.0; 

//Teilchen i 
for (i=start; i< particles;i=i+1) 

{ 

    xwert =0.0; 
    ywert =0.0; 
    zwert =0.0; 

    for (j=0; j< particles;j++) // mit allen anderen Teilchen 
        { 
            //printf("i: %d , %2.3f\n",i,xwert); 

            if (i!=j) { 
                        //    printf("\t%d -- %d\n",i,j); 
                        xwert = xwert + Fx( r[i][0]-r[j][0], r[i][1]-r[j][1], r[i][2]-r[j][2] ); 
                        ywert = ywert + Fy( r[i][0]-r[j][0], r[i][1]-r[j][1], r[i][2]-r[j][2] ); 
                        zwert = zwert + Fz( r[i][0]-r[j][0], r[i][1]-r[j][1], r[i][2]-r[j][2] ); 
                        //    printf("\t%2.5f\n",xwert); 

                        } 
                        else {} 

        } //Ende j-Schleife 

    F[i][0] = xwert; 
    //printf("F[%d][x] : %2.3f\n",i,F[i][0]); 
    F[i][1] = ywert; 
    F[i][2] = zwert; 

}//Ende i-Schleife 

//printf("%2.3f\n",F[1][1]); 
//printf("fertig\n"); 
return 0; 
} 

double Fdotp(double F[][3] , double p[][3]) 
{    
    double ergebnis = 0.0; 
    int i = 0; 

    //for-schleife : Summation ueber ergebnis 
    for (i=start; i< particles; i++) 
        {    
            ergebnis += (F[i][0]*p[i][0] + F[i][1]*p[i][1] + F[i][2]*p[i][2]); 

            } 
    //printf("%2.4f\n",ergebnis); 
    return ergebnis;        

    } 

double sumpypx(double p[][3]) 
{    
    double ergebnis = 0.0; 
    int i = 0; 

    //for-schleife : Summation über ergebnis 
    for (i=start; i< particles; i++) 
        {    
            ergebnis += (p[i][1]*p[i][0]); 

            } 

    return ergebnis;        

    } 

double psq(double p[][3]) 
{    
    double ergebnis = 0.0; 
    int i = 0; 

    //for-schleife : Summation über ergebnis 
    for (i=start; i< particles; i++) 
        {    
            ergebnis += pow(p[i][0],2) + pow(p[i][1],2) + pow(p[i][2],2) ; 

            } 

    return ergebnis;        

    } 

double omz_dt(double F[][3], double p[][3]) 
{    
    double ergebnis = 0.0; 
    int i = 0; 

    //for-schleife : Summation über ergebnis 
    for (i=start; i< particles; i++) 
        {    
            ergebnis += (F[i][0]*p[i][0] - F[i][1]*p[i][1]); 

            } 
    //ergebnis mit -nustir/temp    multiplizieren    
    ergebnis *= - nustir/temp; 

    return ergebnis;        

    } 

double zet(double F[][3], double p[][3]) 
{    
    double ergebnis = 0.0; 

    ergebnis = nugauss* (Fdotp(F,p) - gam* (sumpypx(p)/psq(p)) ); 

    return ergebnis;        

    } 

//funktioniert nicht als Funktion, aber als main schon ?!? 
int RunSimulation() 
//int main() 
{ 

    FILE *datei; 

    //Anfangsorte und Impulse [Teilchenzahl][3 Dimensionen] 
    double rold[particles][3]= {     {0.0   , 0.0  , 0.0}, 
                    {1.166 , 0.01 , 0.01}, 
                    {-1.166, 0.01 ,    0.01} 
                                  }; 
    double rnew[particles][3];                              

    double pold[particles][3];     // zur Zeit von r(t); Anfbed hier rein !    
    double pnew[particles][3];        

    double F[particles][3];        //Kraefte vor dem Sprung 

    //Twirler 
    double omzold = 0.01;        // Anfangsbedingung fuer omz 
    double omznew = 0.0;        // neues omz 

    //Runge-Kutta-Schritte 
    double F1_r[particles][3]; 
    double F2_r[particles][3]; 
    double F3_r[particles][3]; 
    double F4_r[particles][3]; 

    double F1_p[particles][3]; 
    double F2_p[particles][3]; 
    double F3_p[particles][3]; 
    double F4_p[particles][3]; 

    double F1_omz; 
    double F2_omz; 
    double F3_omz; 
    double F4_omz; 

    double tempKraft[particles][3]; //temporaere Teichenkraft 
    double temp_p[particles][3];    //temporaere Impulse (t+delta/2 *t) 
    double temp_F[particles][3];    // temporaere Kräfte (t+delta/2 *t) 
    double temp_omz; 

    //Variablendeklaration 
    int i,j; 
    double progress; 
    clock_t prgstart, stop, stop2; 
    time_t starttime, stoptime, stoptime2; 
    long double timediff; 
    double restzeit=0., restminuten=0.; 
    double restsekunden=0.; 
    char string1[15], string2[15], string3[15], test[20]; 

    printf("... here we go ...  \n"); 

    starttime = time(NULL);     stoptime2=starttime; 
    prgstart=clock();      stop2=prgstart; 

    //Zeitdiffenrz mit 0 initialisieren 
    timediff= 0.; 

    //popen("setterm -cursor off","w");    //blende Cursor aus 
    printf("\033[?25l");            //Variante 2 (besser für bash) 

    //Anfangswerte der Impulse pold setzen 
    for (i=0;i<3;i++) 
        { pold[1][i]=sqrt(temp); 
          pold[2][i]=sqrt(temp); 
            } 
    pold[2][2]=-sqrt(temp); 

    //Arrays initialisieren    
    for (i=0;i< particles; i++) 
    { for (j=0;j< 3; j++) 
        { 
        rnew[i][j] = 0; 
        pnew[i][j] = 0; 
        F[i][j]    = 0; 
        } 
    } 

    printf("Anfangswerte xold & pold - Zeile 1: Teilchen 2, Zeile 2: Teilchen 3\n\n"); 
    for    (i=0;i<3;i++) 
    { 
    printf("Koordinate %d \n %f & %f\n %f & %f\n\n ",i+1,rold[1][i],pold[1][i],rold[2][i],pold[2][i] ); 
    } 
//<<<    
    //Tabellenkopf ausgeben 
    printf("Zeit in s\t| CPU-Zeit in s\t| CPU-Auslastung gesamt\t| Fortschritt \t| Restzeit [m:s]\n"); 

    //datei zum Schreiben neu anlegen Achtung: bestehende Datei wird überschrieben !! 
    datei=fopen(dateiname,"w+"); 

    //Schleife über Zeitschritte 
    for (stepcount=0; stepcount<steps;stepcount++) 
    {    

    //Schritt 1: Gesamtkraft am Ort rold ausrechnen 
    Gesamtkraft(rold,F); 

    //Schleife über alle Teilchen    
    for (i=start; i<particles;i++) 
        { 
        //Schritt 2: Berechnung von F1 

//r(t+delta/2) 
/*x-Kompon.*/    F1_r[i][0] = pold[i][0] + gam * rold[i][1]; 
/*y-Kompon.*/    F1_r[i][1] = pold[i][1]; 
/*z-Kompon.*/    F1_r[i][2] = pold[i][2]; 
//p(t+delta/2) 
/*x-Kompon.*/    F1_p[i][0] = F[i][0] - gam * pold[i][1] - zet(F,pold) * pold[i][0] 
                            - nustir/temp * Fdotp(F,pold) * omzold * pold[i][1]; 
/*y-Kompon.*/    F1_p[i][1] = F[i][1]            - zet(F,pold) * pold[i][1] 
                            + nustir/temp * Fdotp(F,pold) * omzold * pold[i][0]; 
/*z-Kompon.*/    F1_p[i][2] = F[i][2]            - zet(F,pold) * pold[i][2] ; 

//omz(t+delta/2) 
        F1_omz = omz_dt(F,pold); 

        } 

    //Schritt 2: Berechnung von F2 

    //Kraft am Ort F1_r ausrechnen    
    Gesamtkraft(F1_r,tempKraft); 

        //Schleife über alle Teilchen    
    for (i=start; i<particles;i++) 
        { 
//r(t+delta/2) 
/*x-Kompon.*/    F2_r[i][0] = (pold[i][0]+ delta/2 * F1_p[i][0]) + gam * (rold[i][1]+ delta/2 * F1_r[i][1]) ; 
/*y-Kompon.*/    F2_r[i][1] = pold[i][1] + delta/2 * F1_p[i][1]; 
/*z-Kompon.*/    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) 

/*x-Kompon.*/    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]; 
/*y-Kompon.*/    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]; 
/*z-Kompon.*/    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); 

        }        
//printf("%f\n",temp_F[1][1]);                
    //Schritt 3: Berechnung von F3 

    //Kraft am Ort F2_r ausrechnen    
    Gesamtkraft(F2_r,tempKraft); 

        //Schleife über alle Teilchen    
    for (i=start; i<particles;i++) 
        { 
//r(t+delta/2) 
/*x-Kompon.*/    F3_r[i][0] = (pold[i][0]+ delta/2 * F2_p[i][0]) + gam * (rold[i][1]+ delta/2 * F2_r[i][1]) ; 
/*y-Kompon.*/    F3_r[i][1] = pold[i][1] + delta/2 * F2_p[i][1]; 
/*z-Kompon.*/    F3_r[i][2] = pold[i][2] + delta/2 * F2_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 * F2_p[i][0]; 
        temp_p[i][1] = pold[i][1] + delta/2 * F2_p[i][1]; 
        temp_p[i][2] = pold[i][2] + delta/2 * F2_p[i][2]; 

        temp_omz = omzold + delta/2 * F2_omz; 
//p(t+delta/2) 

/*x-Kompon.*/    F3_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]; 
/*y-Kompon.*/    F3_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]; 
/*z-Kompon.*/    F3_p[i][2] = temp_F[i][2]            - zet(temp_F,temp_p) * temp_p[i][2] ; 

//omz(t+delta/2) 
        F3_omz = omz_dt(temp_F,temp_p); 

        }        

    //Schritt 4: Berechnung von F4 

    //Kraft am Ort F3_r ausrechnen    
    Gesamtkraft(F3_r,tempKraft); 

        //Schleife über alle Teilchen    
    for (i=start; i<particles;i++) 
        { 
//r(t+delta/2) 
/*x-Kompon.*/    F4_r[i][0] = (pold[i][0]+ delta * F3_p[i][0]) + gam * (rold[i][1]+ delta * F3_r[i][1]) ; 
/*y-Kompon.*/    F4_r[i][1] = pold[i][1] + delta * F3_p[i][1]; 
/*z-Kompon.*/    F4_r[i][2] = pold[i][2] + delta * F3_p[i][2]; 

        //temporaere Teilchenkraft 
        temp_F[i][0] = F[i][0]+ delta * tempKraft[i][0]; 
        temp_F[i][1] = F[i][1]+ delta * tempKraft[i][1]; 
        temp_F[i][2] = F[i][2]+ delta * tempKraft[i][2]; 
        //temporaere Impulse 
        temp_p[i][0] = pold[i][0] + delta * F3_p[i][0]; 
        temp_p[i][1] = pold[i][1] + delta * F3_p[i][1]; 
        temp_p[i][2] = pold[i][2] + delta * F3_p[i][2]; 

        temp_omz = omzold + delta * F3_omz; 
//p(t+delta/2) 

/*x-Kompon.*/    F4_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]; 
/*y-Kompon.*/    F4_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]; 
/*z-Kompon.*/    F4_p[i][2] = temp_F[i][2]            - zet(temp_F,temp_p) * temp_p[i][2] ; 

//omz(t+delta/2) 
        F4_omz = omz_dt(temp_F,temp_p); 

        }    

    //Schritt 5: Berechnung von rnew,pnew und omznew 

        //Schleife über alle Teilchen    
    for (i=start; i<particles;i++) 
        { 

        rnew[i][0] = rold[i][0] + delta * (F1_r[i][0]/6 + F2_r[i][0]/3 + F3_r[i][0]/3 + F4_r[i][0]/6); 
        rnew[i][1] = rold[i][1] + delta * (F1_r[i][1]/6 + F2_r[i][1]/3 + F3_r[i][1]/3 + F4_r[i][1]/6); 
        rnew[i][2] = rold[i][2] + delta * (F1_r[i][2]/6 + F2_r[i][2]/3 + F3_r[i][2]/3 + F4_r[i][2]/6); 

        pnew[i][0] = pold[i][0] + delta * (F1_p[i][0]/6 + F2_p[i][0]/3 + F3_p[i][0]/3 + F4_p[i][0]/6); 
        pnew[i][1] = pold[i][1] + delta * (F1_p[i][1]/6 + F2_p[i][1]/3 + F3_p[i][1]/3 + F4_p[i][1]/6); 
        pnew[i][2] = pold[i][2] + delta * (F1_p[i][2]/6 + F2_p[i][2]/3 + F3_p[i][2]/3 + F4_p[i][2]/6); 

        omznew = omzold + delta*(F1_omz/6 + F2_omz/3 + F3_omz/3 + F4_omz/6); 
        } 

    //Zeit bestimmen 
    stoptime=time(NULL);    
    //timediff = (float)(stop-stop2) /CLOCKS_PER_SEC; 
    timediff = stoptime-stoptime2; 

    //wenn Zeitdifferenz zur letzten Ausgabe > reload 
    if (timediff >= reload) 
    {    
        progress = (double)stepcount*100/(double)steps ; 
         //Zeit der neuen Referenz 
        stop= clock(); 
        stoptime2=time(NULL); 

        restzeit = ((float)(stoptime2-starttime) * 100. / progress) - (float)(stoptime2-starttime) ; 

        restminuten = fmod(restzeit/60,60.) ; 
        restsekunden = fmod(restzeit ,60.0); 

         timetostr(string1,(stoptime-starttime)); 
        timetostr(string2,((long int)(stop-prgstart) /CLOCKS_PER_SEC)); 
        timetostr(string3,(long int)restzeit); 

//        printf("\n  String : %s\n",string1); 

        //Ausgabe der Zeit zum Start 
        printf("\r %s \t|  %s \t| %3.1f %% \t\t| %3.4f %% \t| %s        ", 
            string1, 
            //CPU-Zeit 
            string2, 
            //Auslastung 
            (float) (100*(float)(stop-prgstart) /CLOCKS_PER_SEC/(stoptime-starttime)), 
            progress, 
            string3 
            ); 

        fflush(stdout); 

    } 

//    //nur alle step/schreibschritt Schritte speichern 
    if (((int) stepcount % ((int)steps/schreibschritt))== 0) 
    { 

    //Schritt 6: abspeichern von rold, pold, omzold, Fold 
    for (i=start; i<particles;i++) 
        { 
        fprintf(datei,"%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t", 
                rold[i][0],rold[i][1],rold[i][2], 
                pold[i][0],pold[i][1],pold[i][2], 
                F[i][0],F[i][1],F[i][2], 
                omzold); 
        } 

    fprintf(datei,"\n");    

    //Umwandeln der Kraft F[1][0] in einen String 
    sprintf(test,"%lf",F[1][0]); 
    //wenn "nan" im String steht, macht das Programm Bockmist und soll sich beenden 
    if (strstr(test,"nan")!=0) { printf("\nError ! NAN !!!\n\a"); goto finalcase;}    

    } 

    //Schritt 7: Setzen der Variablen 
    for (i=start; i<particles;i++) 
        { 
            rold[i][0] = rnew[i][0]; 
            rold[i][1] = rnew[i][1]; 
            rold[i][2] = rnew[i][2]; 

            pold[i][0] = pnew[i][0]; 
            pold[i][1] = pnew[i][1]; 
            pold[i][2] = pnew[i][2]; 
            omzold = omznew; 
        }    

    //schritt beendet auf zur nächsten Iteration        

    } //Ende der Step-Schleife 

//Sprungmarke für den Fehlerfall    
finalcase  : 

    //datei wird geschlossen 
    fclose(datei); 
    printf("\n  ... ready\n\a");    

    //popen("setterm -cursor off","w");    //schalte Cursor an 
    printf("\033[?25h");            //Variante 2 
return 5; 
}    

//funktioniert so nicht 

int main() 
{ 
    int ergebnis; 

    ergebnis=RunSimulation(); 
    printf("\n%d\n",ergebnis); 
return 0;    

}

[*char*jesus]

ohne den code jetzt testen zu können ... ich hab gestern auch mal g++ zum compilieren genommen und auh so 5 oder 6 falsche casts gefunden ... das hatte ich geändert, ohne Erfolg beim meinem Problem zu erzielen ...

ich denke mal irgenwo hat ZuK recht und ich hab irgendwo so einen Array-Fehler drin. Ich vermute es ist eine der String Opereationen ... na ja und dass der code läuft mag unbestritten sein, aber hat bei dir die Ausgabedatei dieselben ergebnisse wenn du es einmal nicht als Funktion aufrufst und einmal doch ? Im letzteren Fall sind nahezu alle Ergebnisse sofort 'nan' bei mir ...

grüße jesus

[Konfusius]

wenn du denkst, daß es an den char-arrays string1,... liegt, dann vergrößere die einfach mal auf 256 zeichen oder so. tritt der fehler dann nicht mehr auf, dann hast du mit an sicherheit grenzender warscheinlichkeit den fehler gefunden.

[*char*jesus]

das hab ich auch schon ausprobiert und es bringt keine Verbesserung ... ich habn auch schon alle Stringoperationen aus meinem Programm auskommentiert (inklusive deklarationen dazu) ...

ich habe nachdem ich gestern ja festgestellt habe, dass es geht, wenn ich alle Variablen global definiere jetzt mal probiert, die Deklaration schrittweise lokal vorzunehmen und konnte das Problem auf die Arrays

double temp_F[particles][3];	// temporaere Kräfte (t+delta/2 *t)
	double tempKraft[particles][3]; //temporaere Teichenkraft 
	double temp_p[particles][3];	//temporaere Impulse (t+delta/2 *t)

eingrenzen. Das witzige ist, je nachdem in welcher Reihenfolge ich die Arrays einführe funktioniert es mal und mal nicht. Obiges Beispiel geht bei mir gerade ... hatz jemand ein Idee ?

Grüße jesus

[*char*jesus]

@Don Carsto: nein, dein Code macht bei mir den gleichen Fehler, wie ich vorhin schon vermutet hatte ... trotzdem danke

jesus

[*char*jesus]

... also ich hab jetzt n neue Idee:

ich initialisiere meine Array jetz wie folgt:

double temp_p[3][3] = { {0},{0} };	
	double temp_F[3][3] = { {0},{0} };

mein Code scheint zu laufen (nur dass ich jetzt mal schauen muss ob ich auch physikalisch rechne).

Ich hatte die Arrays nicht initialisiert, weil sowieso das erste, was mit ihnen passiert ein Schreibzugriff ist ... ich schreibe allerdings nicht auf die Komponente mit Index [0][i] ... ich glaub das könnte ein Problem sein ... weiss ihc aber noch nicht genau ... aber warum lief es dann vorher als main ?

jesus

[Konfusius]

ich habe den code bei mir mit vc++ 6 auch mal kompiliert. bei mir kommt immer NAN raus, egal ob main() als separater funktion oder nicht.

EDIT: vc++ 6 initialisiert variablen immer mit 0, gcc nicht. warscheinlich dividierst du mit uninitialisierten double-array-feldern.

[Konfusius]

//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) 

/*x-Kompon.*/    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];

ich denke ich weiß jetzt, wo der hund begraben liegt: du du benutzt uninitialisierte arrayfelder.

du denkst warscheinlich, du kannst zet(temp_F,temp_p) berechnen, weil die arrays oben bereits initialisiert worden sind. das ist falsch! du initialisierst tatsächlich nur ein feld von temp_p! die felder temp_p[2][] werden erst in den nächsten schleifendurchläufen berechnet. und temp_p[0][] berechnest du überhaupt nicht!

tut mir leid, aber der ganze code ist murks kannst du weitgehend neu schreiben.

[*char*jesus]

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