A simple graphics library in C with BMP and WMF export

/*
	SIMPLE FLOW SIMULATION AROUND A PLATE
*/

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

#include "libtg.h"

#define PI				3.141592653589793239
#define U				1.0
#define X2(x)			((x)*(x))
#define	TUPI			(2.0*PI)
#define FOPI			(2*TUPI)
#define R2PI			(1/TUPI)
#define R4PI			(1/FOPI)
#define	ALPHA			0.0
#define toradian(x)		((x)*PI/180)

/*----------------- Lumpted Vortex Method -------------------------*/
/* Gaussian elimination */

typedef float** Array2d;
typedef float* Array1d;

void solveMatrix(Array2d A, Array1d x, Array1d b, int nx)
{
	int i, j, count, n = nx;
	float ratio=0, temp;

	for(i=0;i<(n-1);i++){
		for(j=(i+1);j<n;j++){
			ratio = A[j][i] / A[i][i];
				for(count=i;count<n;count++) {
					A[j][count] -= (ratio * A[i][count]);
                }
			b[j] -= (ratio * b[i]);
		}
	}

	/* Back substitution */
	x[n-1] = b[n-1] / A[n-1][n-1];

	for(i=(n-2);i>=0;i--){
		temp = b[i];
		for(j=(i+1);j<n;j++){
			temp -= (A[i][j] * x[j]);
		}
		x[i] = temp / A[i][i];
	}
}

Array2d createArray2D(int row, int col)
{
	int error = 0;
	int i,ii;
	Array2d p=(Array2d)malloc(row*sizeof(Array1d));

	if(!p) return p;

	for(i=0; i<row; i++) {
		p[i] = (Array1d)malloc(col*sizeof(float));
		if(!p[i]) { error = 1; break; }
	}

	if(error) {
		for(ii = 0; ii < i; ii++) free(p[ii]);
		free(p);
		p=NULL;
	}
	return p;
}

void freeArray2D(Array2d p, int row)
{
	int i;
	for(i = 0; i < row; i++) free(p[i]);
	free(p); p = NULL;
}

Array1d createArray1d(int nelem)
{
	Array1d p = (Array1d)malloc(nelem*sizeof(float));
	return p;
}

void freeArray1d(Array1d p)
{
	if(p) free(p); p = NULL;
}

void getVelocity_UPlate(Array1d gam, Array1d xv, float uu, float vv, float xs, float ys, int npan, float* velu, float* velv)
{
	int j=0;
	float r2=0;

	*velu = *velv = 0;

	for(j = 0; j < npan; j++)
	{
		r2 = (xs-xv[j])*(xs-xv[j])+ys*ys;
		*velu += -gam[j]*ys/r2;
		*velv += -gam[j]*(xs-xv[j])/r2;
	}
	*velu = (float)(uu + *velu*R2PI);
	*velv = (float)(vv - *velv*R2PI);
}

void Euler_UPlate(Array1d gam, Array1d xv, float uu, float vv, float xs, float ys, int npan, float dell, float* xm1p, float* ym1p)
{
	float um,vm,qm,ump,vmp,qmp;

	um = vm = 0;
	getVelocity_UPlate(gam, xv, uu, vv, xs, ys, npan, &um, &vm);
	qm = (float)sqrt(um*um + vm*vm);

	/* 1st order forward difference*/
	*xm1p = xs + um/qm * dell;
	*ym1p = ys + vm/qm * dell;

	ump = vmp = 0;
	getVelocity_UPlate(gam, xv, uu, vv, *xm1p, *ym1p, npan, &ump, &vmp);
	qmp = (float)sqrt(ump*ump + vmp*vmp);

	/* central difference method : improved method*/
	*xm1p = (float)(xs + 0.5 * ( um / qm + ump / qmp ) * dell);
	*ym1p = (float)(ys + 0.5 * ( vm / qm + vmp / qmp ) * dell);
}

/*
	dy/dx = f(x,y) ... x = x0, y = y0
	
	k1 = h f(xi, yi)
	k2 = h f(xi + h/2, yi + k1/2)
	k3 = h f(xi + h/2, yi + k2/2)
	k4 = h f(xi + h  , yi + k3)
	delyi = 1/6 ( k1 + 2k2 + 2k3 + k4)
	
	xi+1 = xi + h
	yi+1 = yi + delyi
*/

void RungeKutta4_UPlate(Array1d gam, Array1d xv, float uu, float vv, float xs, float ys, int npan, float dell, float* xm1p, float* ym1p)
{
	float um=0,vm=0,k1=0, k2=0, k3=0, k4=0, delyi;
	float h=dell, xi=xs,yi=ys;

	getVelocity_UPlate(gam, xv, uu, vv, xi, yi, npan, &um, &vm);
	k1 = h*vm/um;

	um=vm=0;
	getVelocity_UPlate(gam, xv, uu, vv, xi+h/2, yi+k1/2, npan, &um, &vm);
	k2 = h*vm/um;

	um=vm=0;
	getVelocity_UPlate(gam, xv, uu, vv, xi+h/2, yi+k2/2, npan, &um, &vm);
	k3 = h*vm/um;

	um=vm=0;
	getVelocity_UPlate(gam, xv, uu, vv, xi+h, yi+k3, npan, &um, &vm);
	k4 = h*vm/um;

	delyi = (float)(1./6.*(k1 + 2*k2 + 2*k3 + k4));
	*xm1p = xi+h;
	*ym1p = yi+delyi;
}

void U_Plate(TG_PlotZone_Ptr pz, float alpha, float minx, float miny,
	float maxx, float maxy, int nsgx, int nsgy)
{
	int npan = 20;
	float chord = 1, rhs=0;
	float alrad = (float)toradian(alpha);
	int i=0,j=0,k=0;
	float xs=0, ys=0, velu=0, velv=0, r2=0, slop=0, uu=0, vv=0;
	float delx = (maxx - minx)/(float)(nsgx-1);
	float dely = (maxy - miny)/(float)(nsgy-1);
	float dell = delx, xm1p=0, ym1p=0;
	TG_Plot_Ptr p = &pz->plot;
	
	Array1d xv, xc;
	Array2d A;
	Array1d b, gam, gridy;

	xv = createArray1d(npan);
	xc = createArray1d(npan);
	A = createArray2D(npan, npan);
	b = createArray1d(npan);
	gam = createArray1d(npan);
	gridy = createArray1d(nsgy);

	rhs = (float)(-2*PI*U*sin(alrad));
	for(i = 0; i < npan; i++) 
	{
		xv[i] = (float)(1+4*i)*chord/(npan*4);
		xc[i] = (float)(3+4*i)*chord/(npan*4);
		b[i]  = (float)rhs;
	}
	
	for(i = 0; i < npan; i++)
	{
		for(j = 0; j < npan; j++)
		{
			A[i][j] = 1/(xc[i]-xv[j]);
		}
	}
	
	solveMatrix(A, gam, b, npan);
	
	for(i = 0; i < nsgy; i++) gridy[i] = miny+dely*i;
	
	uu = (float)(U*cos(alrad));
	vv = (float)(U*sin(alrad));

	p->gsave();
	p->makepen(pz, TG_GetColor(TG_MAGENTA), 0.01);
	p->moveto(pz, xs,ys);
	p->lineto(pz, (float)chord,ys);
	p->deletepen();
	p->grestore();

	for(j = 0; j < nsgy; j++) 
	{
		xs = minx;
		ys = gridy[j];
		p->moveto(pz,xs,ys);

		while(1) 
		{
			Euler_UPlate(gam, xv, uu, vv, xs, ys, npan, dell, &xm1p, &ym1p);
			/*RungeKutta4_UPlate(gam,xv,uu,vv,xs,ys,npan, dell, &xm1p, &ym1p);*/
			if(xm1p < minx || xm1p > maxx || ym1p < miny || ym1p > maxy) break;
			p->lineto(pz,xm1p, ym1p);
			xs = xm1p, ys = ym1p;
		}
	}

	freeArray1d(xv);
	freeArray1d(xc);
	freeArray2D(A, npan);
	freeArray1d(b);
	freeArray1d(gam);
	freeArray1d(gridy);
}

#define BMP

int main()
{
	int xwid=500, xhgt=500;
	TG_PlotZone_Ptr zone=NULL;
	TG_Workspace_Ptr workspace = TG_InitTGLibrary("Letter");
	float minx=-1.0f, maxx=1.0f, miny=-1.0f, maxy=1.0f;
	
	if(!workspace) return 1;
	
	zone = &workspace->zone;
	TG_SetFrame(&zone->frame, 0, 0, 3, 3, minx, maxx, miny, maxy);
	
#ifdef BMP
	TG_BMP_OpenExport(zone, "u-plate.bmp", xwid, xhgt, 1);
#else
	TG_WMF_OpenExport(zone, "u-plate.wmf");
#endif

	U_Plate(zone, 10, minx,miny,maxx,maxy, 50, 40);

#ifdef BMP
	TG_BMP_CloseExport();
#else
	TG_WMF_CloseExport();	
#endif
	
	return 0;
}