problem mit function

[ZuK]

Schade dass du uns den code nicht zeigen kannst. Sonst hätten wir vielleicht helfen können

[*char*jesus]

Hallo,

ok ich häng den Code jetzt mal an, ich hoffe ... ihr mögt mich danach noch ... das Beispiel hier funktioniert ... das nicht funktionierende ist auskommentiert .. wenn ich die funktion main umbenenne nach RunSimulation und dann ganz unten eine neue main einfüge (so wie auskommentiert), so geht es nicht mehr ... ich versteh's net ...

jesus

//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 int steps = 1000000.0;  	 // Anzahl Schritte
long int 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 int timediff;
	double restzeit=0., restminuten=0.;
	int 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,((float)(stop-prgstart) /CLOCKS_PER_SEC));
		timetostr(string3,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 (( stepcount % (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]

ach so vielleicht ist's wichtig zu wissen ...
(a) die ausgabedatei ist ein paar MB groß ...
(b) ich benutze Linux ... daher stehen im Code am Anfang und am Ende je eine EscapeSequenz für die Konsole (schaltet den Cursor aus und wieder an) ...

jesus

[ZuK]

Got some bad news.
wenn du main so umbaust funktioniert's

int main()
{
    int ergebnis;
    char bbb[9];   // muss grösser 8 sein sonst stürzts bei mir ab
    ergebnis=RunSimulation();
    printf("\n%d\n",ergebnis);
    return 0;   

}

Das lässt auf einen array-bounds Fehler schliessen. Happy debugging.
Kurt

[*char*jesus]

was ändert das jetzt am code, wenn du einen String definierst, der gar nicht benutzt oder übergeben wird ?

jesus

[ZuK]

Es ändert auf jedenfall das layout der daten am Stack. Und wenn nur die Deklaration einer Dummy Variablen das Ergebnis deiner Berechnung verändert dann bedeutet das für mich dass irgendwo auf einen falschen Pointer zugreifst oder über Arraygrenzen hinaus zugegriffen wird.
Du kannst ja einmal ausprobieren ob es sich bei dir genauso auswirkt.
Das sind natürlich alles nur Vermutungen bei mir läufts auch wenn ich die Dummy-Variable so einfüge

double tempKraft[particles][3]; //temporaere Teichenkraft
    double temp_p[particles][3];    //temporaere Impulse (t+delta/2 *t)
    char bbb[9];   // wenn hier dann geht's
    double temp_F[particles][3];    // temporaere Kräfte (t+delta/2 *t)
//    char bbb[9]; // da nicht mehr
    double temp_omz;

Wenn es sich bei dir genauso verhält wär das ein Indiz dass es mit den Variablen in dieser Gegend zu tun hat.
BTW mein compiler ist gcc 3.3.5 und das Ergebnis ist immer das gleiche auch bei verschiedenen Optimierungsparametern -G1 .. -G3. Ob die Berechnungen stimmen kann ich nicht beurteilen.

[*char*jesus]

mmmh bei mir ändert es gar nichts ob ich diese Variable deklariere oder nicht ... ich hab's genauso wie du vorschlägst ... ich benutze übrigens auch gcc 3.3.5

jetzt habe ich aber eine andere Möglichkeit gefunden, die geht ... ich habe alle Variablen und Arrays global definiert und plötzlich geht es auch als Funktion ... wie ist denn das : wird der Speicher nur dann richtig reserviert, wenn ich alle Varialben und Arrays in main oder global definiere ?

edit: nee jetzt kriege ich immer beim 2. Start vom Prog einen Segmenation Fault ... *heul*

jesus

[*char*jesus]

was passiert denn wenn ich einen String (also char test[20]) nehme und mit fprintf einen Double als Zeichenkette hineinschreibe ? wenn die zeichenkette größer als 20 zeichen ist, wird dann abgeschnitten ? Oder greift sich der String benachbarten Speicher, den er nicht haben darf ?

jesus

[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