A simple graphics library in C with BMP and WMF export

/*
	SIMULATIOIN BASIC POTENTIAL FLOW IN FLUIC MECHANICS
*/

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

#include "libtg.h"

typedef enum PotentialFlowIndex_ 
{
	UNIFORM_SOURCE = 0,
	UNIFORM_VORTEX = 1,
	UNIFORM_DIPOLE = 2,
	UNIFORM_RNKOVL = 3 
} 
	PotentialFlowIndex_t;
	
#define U				1.0
#define Q				1.0
#define M				1.0
#define RK				1.5
#define RKM				0.2
#define GAM				2.2
#define PI				3.141592653589793239
#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)

/*---------- SET DEGREE ENVIRONMENT FOR UNIFORM FLOW ----------*/

static float degree_of_uniformflow = 0;
static float angle_of_attack = 0;

float get_degree(void) { return degree_of_uniformflow; }
void set_degree(float deg) 
{ 
	degree_of_uniformflow = deg; 
	angle_of_attack = (float)toradian(degree_of_uniformflow); 
}

/*---------------UNIFORM + SOURCE------------------------------*/
float dist(float x0, float y0, float x1, float y1) 
{
  return (float)sqrt( (x1-x0)*(x1-x0) + (y1-y0)*(y1-y0) );
}

float uus(float x0, float y0, float x1, float y1) 
{
  float R = (float)dist(x0,y0,x1,y1);
  return (float)( U*cos(angle_of_attack) + Q*(x1-x0)/(TUPI*R*R));
}

float vus(float x0, float y0, float x1, float y1) 
{
  float R = (float)dist(x0,y0,x1,y1);
  return (float)( U*sin(angle_of_attack) + Q*(y1-y0)/(TUPI*R*R));
}

/*--------------UNIFORM + DIPOLE -----------------------------*/

float uud(float x0, float y0, float x1, float y1) 
{
	return (float)(U*cos(angle_of_attack) - 
		M*R4PI*(X2(y1-y0)-2*X2(x1-x0))/pow(X2(y1-y0)+X2(x1-x0), 5/2.));
}

float vud(float x0, float y0, float x1, float y1) 
{
	return (float)(U*sin(angle_of_attack) - 
		3*M*R4PI*((y1-y0)*(x1-x0))/pow(X2(y1-y0)+X2(x1-x0), 5/2.));
}

/*--------------UNIFORM + VORTEX -----------------------------*/

float uuv(float x0, float y0, float x1, float y1) 
{
  return (float)(U*cos(angle_of_attack) + 
	  GAM*R2PI*(y1-y0)/(X2(x1-x0)+X2(y1-y0)));
}

float vuv(float x0, float y0, float x1, float y1) 
{
  return (float)(U*sin(angle_of_attack) - 
	  GAM*R2PI*(x1-x0)/(X2(x1-x0)+X2(y1-y0)));
}

/*---------------UNIFORM + SOURCE + SINK ----------------------*/
/* Reference: FLUID MECHANICS, 2nd Ed. by A.K. Mohanty         */

float usk(float x0, float y0, float x1, float y1)
{
	return (float)(U*cos(angle_of_attack)+
		RK*R2PI*((x1+RKM)/(X2(y1)+X2(x1+RKM))-(x1-RKM)/(X2(y1)+X2(x1-RKM))));
}

float vsk(float x0, float y0, float x1, float y1)
{
	return (float)(U*sin(angle_of_attack)
		-RK*R2PI*(y1/(X2(y1)+X2(x1-RKM))-y1/(X2(y1)+X2(x1+RKM))));
}

float (*uu[])(float,float,float,float)
     = { uus, uuv, uud, usk };
float (*vv[])(float,float,float,float)
     = { vus, vuv, vud, vsk };

void U_Potential(TG_PlotZone_Ptr p, float deg, float minx, float miny,
	float maxx, float maxy, int nsgx, int nsgy, PotentialFlowIndex_t ipot )
{
	int i=0, j=0;
	float delx = (maxx - minx)/(float)(nsgx-1);
	float dely = (maxy - miny)/(float)(nsgy-1);
	float dell = delx, old_degree=0;
	float xm1p, ym1p, um, vm, ump , vmp , qm , qmp;
	float x0, y0, x1, y1, startx, *gridy=NULL;

	gridy = (float*)malloc(sizeof(float)*nsgy);
	if(!gridy) return;

	gridy[0] = miny;
	for(i = 0; i < nsgy; i++) gridy[i] = gridy[0]+(float)i*dely;

	//p->fillrect1(p, minx, miny, maxx, maxy, BLACK, YES, CYAN);
	x0 = 0.0000;
	y0 = 0.0000;
	startx = minx;
	old_degree = get_degree();
	set_degree(deg);
  
	for(j = 0; j < nsgy; j++) 
	{
		x1 = startx;
		y1 = gridy[j];
		p->plot.moveto(p,x1,y1);

		while(1) 
		{
			um = (uu[ipot])(x0, y0, x1, y1);
			vm = (vv[ipot])(x0, y0, x1, y1);
			qm = (float)sqrt(um*um + vm*vm);

			/* 1st order forward difference*/
			xm1p = x1 + um/qm * dell;
			ym1p = y1 + vm/qm * dell;
			ump = (uu[ipot])(x0, y0, xm1p, ym1p);
			vmp = (vv[ipot])(x0, y0, xm1p, ym1p);
			qmp = (float)sqrt(ump*ump + vmp*vmp);

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

			if(xm1p < minx || xm1p > maxx || ym1p < miny || ym1p > maxy) break;
			p->plot.lineto(p,x1, y1);
			x1 = xm1p, y1 = ym1p;
		}
	}
	free(gridy);
	set_degree(old_degree);
}

#define BMP

int main()
{
	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;
	int xwid=500, xhgt=500;

	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, "uniform-source.bmp", xwid, xhgt, 1);
#else
	TG_WMF_OpenExport(zone, "uniform-source.wmf");
#endif

	U_Potential(zone, 0., minx,miny,maxx,maxy, 50, 40, UNIFORM_SOURCE);

#ifdef BMP
    TG_BMP_CloseExport();	
#else
    TG_WMF_CloseExport();
#endif

#ifdef BMP
	TG_BMP_OpenExport(zone, "uniform-vortex.bmp", xwid, xhgt, 1);
#else
	TG_WMF_OpenExport(zone, "uniform-vortex.wmf");
#endif

	U_Potential(zone, 0., minx,miny,maxx,maxy, 50, 40, UNIFORM_VORTEX);

#ifdef BMP
    TG_BMP_CloseExport();	
#else
    TG_WMF_CloseExport();
#endif

#ifdef BMP
	TG_BMP_OpenExport(zone, "uniform-dipole.bmp", xwid, xhgt, 1);
#else
	TG_WMF_OpenExport(zone, "uniform-dipole.wmf");
#endif

	U_Potential(zone, 0., minx,miny,maxx,maxy, 50, 40, UNIFORM_DIPOLE);

#ifdef BMP
    TG_BMP_CloseExport();	
#else
    TG_WMF_CloseExport();
#endif

#ifdef BMP
	TG_BMP_OpenExport(zone, "uniform-rnkovl.bmp", xwid, xhgt, 1);
#else
	TG_WMF_OpenExport(zone, "uniform-rnkovl.wmf");
#endif

	U_Potential(zone, 0., minx,miny,maxx,maxy, 50, 40, UNIFORM_RNKOVL);
	
#ifdef BMP
    TG_BMP_CloseExport();	
#else
    TG_WMF_CloseExport();
#endif

	return 0;
}