[C] Jak dokładnie działa ten kod


(Bulii) #1


(Enterbios) #2

Ciężko będzie bo w programie masz stałe(zakładam że stałe?) zdefiniowane gdzieś indziej np. NX, tablica FLAG, C_DONT, C_BND. Z kolei c nie różni się aż tak bardzo od matlaba.

t+=(PI/(float)180) tego w matlabie nie ma (operatora +=) wiec w matlabie będzie to tak t = t + (PI/180), castowanie na float też możesz usunąć,

warunek jest tak samo,

pętla odpowiada matlabowej for i=0:1:NX, a dokładniej for i=1:NX gdyż matlab zdaje się indeksuje wektory od 1, a w C indeksacja tablic zaczyna się od 0.

Reszta to obliczenia nic mi nie mówiące gdyż nie pofatygowałeś się nawet o napisanie na co jest to algorytm, a chyba nie powiesz mi że znalazłeś jakiś kod w C i chcesz dla zabawy go przepisać do matlaba?

A co do kończenia się pętli:

zapis

for(;:wink:{

wykonuje instrukcje w bloku pomiędzy klamrami

}

zapis bez klamer

for(;:wink:

instrukcja1

instrukcja2

wykonuje w pętli tylko pierwszą instrukcje.


(Bulii) #3

Oto caly kod programu:

#include 

#include 

#include 

#include 


#include 


#define CI_RED COLOR_RED

#define CI_ANTI_ALIAS_GREEN 16

#define CI_ANTI_ALIAS_YELLOW 32

#define CI_ANTI_ALIAS_RED 48



GLenum rgb, doubleBuffer=1, windType;

GLint windW, windH;


#include "tkmap.c"


GLenum mode=0;

GLint size;


#ifndef PI

#define PI 3.14159265358979323846

#endif


#define PX 0

#define PY 1

#define PZ 2


float point[3] = { 0.0, 0.0, 0.0 };


float point1[3] = { 0.0, 0.0, 0.0 };

float point2[3] = { 0.0, 0.0, 0.0 };



float t,dt;


#define WWIDTH 500

#define WHEIGHT	250


#define NX 1000

#define SWIDTH 1000

float DX=(SWIDTH/(float)NX);

float DT=0.05;

float VEL=107.8;


float U[NX+1]={0}; // U - wych.

float V[NX+1]={0}; // velocity

float T[NX+1]={0}; // teta




int	FLAG[NX+1]={0}; // komórki znaczone



#define C_BND 0x00000001 /* boundary cells */

#define C_DONT 0x00000002 /* do not calculate U,T,E for that cells */					




static void Key(unsigned char key, int x, int y)

{


    switch (key) {

	case 'q':

		exit(0);


	break;


	case ' ':

		t=0;


	break;

      default:

	return;

        }


    glutPostRedisplay();

}


static void Init(void)

{

    GLint i;


    glClearColor(0.0, 0.0, 0.0, 0.0);


    glBlendFunc(GL_SRC_ALPHA, GL_ZERO);


    mode = GL_FALSE;

    size = 1;



// warunki brzegowe


	for(i=0;i<19;i++)

	{

		FLAG[i] |= C_DONT;

	 //	FLAG[NX-i] |= C_DONT;

	}


//	for(i=0;i<20;i++)

// FLAG[i] |= C_DONT;


	for(i=0;i<20;i++)

	{

	   	//FLAG[i] |= C_BND;

		FLAG[NX-i] |= C_BND;

	}



}


static void Reshape(int width, int height)

{


    windW = (GLint)width;

    windH = (GLint)height;


    glViewport(0, 0, width, height);


    glMatrixMode(GL_PROJECTION);

    glLoadIdentity();

    //gluOrtho2D(-windW/2, windW/2, -windH/2, windH/2);

    gluOrtho2D(-1, 1, -1, 1);

    glMatrixMode(GL_MODELVIEW);

}



int init=0,line=0,step=0,k=0;



static void Draw(void)

{

   int i;


//if(init==0)



	glClearColor(1.0, 1.0, 1.0, 0.0);      

  	glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);    



	glDisable(GL_BLEND);     

	glDisable(GL_POINT_SMOOTH);     


	glColor3f(0.0,0.0,0.0);


	glPointSize(4);

	glLineWidth(2);     



	step++;


	point1[PX]=-1+(2/(float)NX);


 	glBegin(GL_POINTS);     

		for(i=1;i
		{

			point1[PY]=U[i];

			point1[PX]+=(2/(float)NX);

			glVertex3fv(point1);

	 	}

	glEnd();      




    if (doubleBuffer) {

	glutSwapBuffers();

    }

}


static GLenum Args(int argc, char **argv)

{

    rgb = GL_TRUE;

    doubleBuffer = GL_TRUE;


    return GL_TRUE;

}



void bndconditions(void)

{

	int i;


	for(i=0;i
	if(FLAG[i] == C_BND)

	{


	/* V[i]=V[i-1];

	   	T[i]=T[i-1];

	   	U[i]=U[i-1];

	*/

	   	V[i]=0;

	   	T[i]=0;

	   	U[i]=0;  


	}

}


[code]

void idle(void)

{

	int i;



t+=(PI/(float)180); //t=t+PI/(float)180;

 if(t
  	for(i=0;i
		begin

		if(FLAG[i]==C_DONT)

			begin

			V[i]=sin(t)/8;

			end

		end

end


		for(i=0;i
		if( (FLAG[i] != (C_DONT|C_BND)) && (FLAG[i+1] != (C_DONT|C_BND)) )

		{		  

			V[i] = (V[i+1]+V[i-1])/2 + VEL*(DT/DX)*(T[i+1]-T[i]);

			T[i] = (T[i+1]+T[i-1])/2 + VEL*(DT/DX)*(V[i+1]-V[i]);


		}



 	bndconditions();



		for(i=0;i
		if(FLAG[i] != (C_DONT|C_BND))

					U[i] = U[i]+V[i]*DT;			



	for(i=0;i
		if( (FLAG[i] != (C_DONT|C_BND)) && (FLAG[i+1] != (C_DONT|C_BND)) )

				T[i] = (U[i+1]-U[i-1])/(2*DX);	



	glutPostRedisplay();

}


int main(int argc, char **argv)

{

    glutInit(&argc, argv);


    if (Args(argc, argv) == GL_FALSE) {

	exit(1);

    }


    windW = WWIDTH;

    windH = WHEIGHT;

    glutInitWindowPosition(150, 150); glutInitWindowSize( windW, windH);


    windType = (rgb) ? GLUT_RGB : GLUT_INDEX;

    windType |= (doubleBuffer) ? GLUT_DOUBLE : GLUT_SINGLE;

    glutInitDisplayMode(windType);


    if (glutCreateWindow("Wave 1D") == GL_FALSE) {

	exit(1);

    }


    InitMap();


    Init();


    glutKeyboardFunc(Key);

    glutReshapeFunc(Reshape);

    glutDisplayFunc(Draw);

    glutIdleFunc(idle);

    glutMainLoop();


	return 0;

}

Program jest z książki i działa dobrze. Rozwiązuje on równanie falowe drugiego rzędu(jako układ równań sprzężonych) i wyświetla na ekranie w czasie rzeczywistym biegnącą po ekranie falę w jednym wymiarze (przesuwające się wychylenie na strunie )

Wzorami kolejno:

V = (V[i+1]+V[i-1])/2 + VEL*(DT/DX)*(T[i+1]-T__);

T = (T[i+1]+T[i-1])/2 + VEL*(DT/DX)*(V[i+1]-V__);

U = U +V *DT;

przybliżane są wartości (dla komórki bieżącej) odchylenia kątowego, zmiany wychylenia w czasie i wychylenie w następnym kroku czasowym.

Nie wiem czy dobrze rozumiem tą pętlę czasową.

Czyli ten krok czasowy t idzie od 0 w nieskończoność, a od 0 do pi co krok pi/180 wykonuje się operacja V__=sin(t)/8 ?

Gdy t>pi wykonują się kolejne pętle? i co wtedy daje zmiana parametru t większego od pi w kolejnych krokach? Nie mogę zrozumieć zasady działania tego:/ Byłby ktoś w stanie mi to wytłumaczyć?


(Enterbios) #4

Klasyczny kod starych weteranów C, czyli najlepszy sposób na zmarnowanie sobie całego dnia na analizę :smiley: współczuje że musisz coś takiego analizować.

Z tego co udało mi się zauważyć to tak jak napisałeś t rośnie w nieskończoność, i jest resetowane tylko w przypadku kiedy ktoś kliknie spacje, krok z jakim rośnie t w każdej iteracji jest jak już zauważyłeś PI/180. Z tego co rozumiem wyświetlanych jest NX =1000 punktów,

Tablica jest NX-elementowa tak więc numer elementu tablicy odpowiada współrzędnej x na osi, natomiast wartość n-tego elementu jest współrzędną y.

Flagi C_BND i C_DONT są ustawiane dla pierwszych i ostatnich 20 punktów, a następnie parametry dla tych punktów końcowych (z C_BND)są zerowane w funkcji bndconditions();

nie widzę animacji ale obstawiam że na końcu fali jest odcinek prosty z tego powodu (moglem gdzieś coś pominąć).

Zanim t będzie większe od PI, jest obliczane V dla pierwszych 20 punktów. Generalnie z tymi flagami strasznie ktoś namieszał przez co kod jest zupełnie nie czytelny.

Dla Ciebie najważniejsze jest to że w każdym wywołaniu funkcji idle, dla wszystkich punktów (prócz pierwszej i ostatniej 20) obliczane są parametry, V, z nich T oraz U, potem z nowych U znowu jest wyliczane T (nie rozumiem po co T było wyliczane wcześniej skoro nie było nigdzie użyte, po czym zostało obliczone ponownie tym razem korzystając z U). Po wyliczeniu parametrów następuje przerysowanie okna. Nie wiem czy to dobry pomysł aby przerabiać ten kod, coś mi się wydaje że łatwiej będzie napisać własny w oparciu o wzory, które przecież masz.

edit:

ok czaje po co jeszcze raz jest wyliczane T, do następnej iteracji :smiley:


(Bulii) #5

Ok dzięki za pomoc. Sam pisze tylko ze mi nie wychodzi dlatego wzoruje się na tym kodzie :?