/*

	X program to generate conformal mesh for
        a channel flow. It compiles on my linux
        setup (installed from Suse distribution)
        but note the need for xview to be installed

	Once compiled (use makefile) it should be
        fairly clear what different buttons do.

        The furrows have a top and bottom displacement
        and have flat sections joined by sinusoidal sections

        Some controls are not implimented (phase for instance)
        but it sohuld be simple to activate them if you want to.

*/



#include <xview/xview.h>
#include <xview/panel.h>
#include <xview/canvas.h>
#include <xview/openmenu.h>
#include <xview/scrollbar.h>
#include <sys/ioctl.h>
#include <stdio.h>
#include <xview/textsw.h>
#include <pwd.h>
#include <math.h>


#define cut(a,b,c)	(((a==c)||(((a-c)*(b-c))<0.))?0:1)
#define HUMP_LENGTH	(1.5)
#define MAX_Y_GRID	129
#define MAX_X_GRID	514
#define PI		(3.14159265358979)

Pixwin           *mypw;
Display * display;
GC        gc;
XGCValues gc_val;
XID       xid;
XFontStruct *font_info;
char *fontname="9x15";
Rect 			base_rect;

Frame		base_frame;
Canvas		canvas;

Panel_item	x_points,y_points,cycles_slider;
Panel_item	canvas_auto_clear,
			clear_data,
			iter_item,
			method_item,
			name_item,
			top_depth_item,
			bot_depth_item,
			hump_length_item,
			flat_length_item,
			phase_item,
			omega_item,
			error_item,
			length_item,
    			number_item ,
			file_item,
			free_item;

static short icon_data[] = {
#include "use.icon"
};
static mpr_static(icon_pr, 64, 64, 1, icon_data);

/* panel notify procs */
static void	modify_canvas();
static void	clear_canvas();
static void	draw_canvas();
static void	quit();
static void     help_msg();
/* canvs procs */
static void	repaint_canvas();
static void	handle_event();
/* utilities */
static void	proc_draw_coordinates();
static void	init_panel();
static void	draw_line();
static void 	file_coordinates();
static void 	read_coordinates();
static void 	read_conditions();
static void 	plot_routine();
static void 	frame_plot_area();
static int 	x_pix();
static int 	y_pix();
static int display_sol_active={0};
static double 	absol();
static void 	evaluate_coefficients();
static void 	relax();
static void 	initialise_walls();
static void 	initialise_x_y();

static Frame popup_frame;
static Textsw popup_text;
static Panel popup_panel;
static Panel_item popup_msg_item;
static void popup_msg(), init_popup_msg(), popup_textsw_done();
char *xlabel={"X"};
char *title={"Y"};
static double xmax,ymax,xmin,ymin,xstep,ystep;
double omega;
char 	top_depth_string[20]		={"1.0"},
 	bot_depth_string[20]		={"1.0"},
 	hump_length_string[20]		={"2.0"},
 	flat_length_string[20]		={"2.0"},
	phase_string[20]		={"0.0"},
	omega_string[20]		={"1.0"},
	error_string[20]		={"0.0"},
	length_string[20]		={"8.0"},
	free_string[20]			={"0.0"},
	file_string[60]			={"../"},
	iter_string[20]			={"0"};
int width,height;
int x_off={40};
int y_off={40};
char *get_name();

double 	x_upper_wall[MAX_X_GRID],y_upper_wall[MAX_X_GRID];
double 	x_lower_wall[MAX_X_GRID],y_lower_wall[MAX_X_GRID];
double 	x[MAX_Y_GRID][MAX_X_GRID],y[MAX_Y_GRID][MAX_X_GRID];
double 	alpha[MAX_Y_GRID][MAX_X_GRID],
	beeta[MAX_Y_GRID][MAX_X_GRID],
	gama[MAX_Y_GRID][MAX_X_GRID];
double 	jacobian[MAX_Y_GRID][MAX_X_GRID];
double 	p[MAX_Y_GRID][MAX_X_GRID],
	q[MAX_Y_GRID][MAX_X_GRID];
double length,top_depth,bot_depth,flat_length,hump_length,phase,free_value;
int m,n,n_iterates;


main(argc, argv)
int argc;
char **argv;
{  
	xv_init(XV_INIT_ARGC_PTR_ARGV,&argc,argv,NULL);
  base_frame = xv_create(0, FRAME,
	    XV_LABEL,    "Boundary Fitted Coordinates",
	    FRAME_ICON,	    icon_create(ICON_IMAGE, &icon_pr, 0),
	    FRAME_ARGS,     argc, argv,0);
	base_rect.r_top=base_rect.r_left=50;
	base_rect.r_width=1000;
	base_rect.r_height=600;
	frame_set_rect(base_frame,&base_rect);

    init_panel(base_frame);    

	canvas=xv_create(base_frame,CANVAS,
            WIN_CONSUME_EVENT,  LOC_DRAG,
            WIN_IGNORE_EVENTS,
                WIN_ASCII_EVENTS,WIN_UP_ASCII_EVENTS,WIN_UP_EVENTS,
            WIN_EVENT_PROC,             modify_canvas,
            CANVAS_WIDTH,               1000,
            CANVAS_HEIGHT,              600,
            CANVAS_AUTO_SHRINK,         TRUE,
            CANVAS_FIXED_IMAGE,         FALSE,
            CANVAS_X_PAINT_WINDOW,      TRUE,
		0);
	xv_set(canvas,CANVAS_RESIZE_PROC,modify_canvas,0);
	display=(Display *)xv_get(base_frame,XV_DISPLAY);
        xid=(XID)xv_get(canvas_paint_window(canvas),XV_XID);
        gc=XCreateGC(display,xid,0,&gc_val);
        XSetForeground(display,gc,0);
	window_fit(base_frame);
	mypw=(Pixwin *)canvas_pixwin(canvas);
	font_info=XLoadQueryFont(display,fontname);

    read_conditions();
    window_main_loop(base_frame);
    exit(0);
}

static void init_panel(base_frame)
Frame	base_frame;
{   Panel	panel;
	char *s;

    panel = xv_create(base_frame, PANEL,  NULL);

    name_item 	 = xv_create(panel, PANEL_TEXT,
        PANEL_LABEL_STRING, "Name: ",
        PANEL_VALUE_DISPLAY_LENGTH, 8,
        PANEL_VALUE_STORED_LENGTH, 20,
        PANEL_VALUE, "Fred",0);
    top_depth_item 	= xv_create(panel, PANEL_TEXT,
        PANEL_LABEL_STRING, "Furrows: Top",
        PANEL_VALUE_DISPLAY_LENGTH, 5,
        PANEL_VALUE_STORED_LENGTH, 20,
        PANEL_VALUE, top_depth_string,0);
    bot_depth_item 	= xv_create(panel, PANEL_TEXT,
        PANEL_LABEL_STRING, ":Bottom",
        PANEL_VALUE_DISPLAY_LENGTH, 5,
        PANEL_VALUE_STORED_LENGTH, 20,
        PANEL_VALUE, bot_depth_string,0);
    hump_length_item 	= xv_create(panel, PANEL_TEXT,
        PANEL_LABEL_STRING, ":hump",
        PANEL_VALUE_DISPLAY_LENGTH, 5,
        PANEL_VALUE_STORED_LENGTH, 20,
        PANEL_VALUE, hump_length_string,0);
    flat_length_item 	= xv_create(panel, PANEL_TEXT,
        PANEL_LABEL_STRING, ":flat",
        PANEL_VALUE_DISPLAY_LENGTH, 5,
        PANEL_VALUE_STORED_LENGTH, 20,
        PANEL_VALUE, hump_length_string,0);
    phase_item 	= xv_create(panel, PANEL_TEXT,
        PANEL_LABEL_STRING, "Phase(deg): ",
        PANEL_VALUE_DISPLAY_LENGTH, 5,
        PANEL_VALUE_STORED_LENGTH, 20,
        PANEL_VALUE, phase_string,0);
    omega_item 	= xv_create(panel, PANEL_TEXT,
        PANEL_LABEL_STRING, "Omega:",
        PANEL_VALUE_DISPLAY_LENGTH, 5,
        PANEL_VALUE_STORED_LENGTH, 20,
        PANEL_VALUE, omega_string,0);
    error_item 	= xv_create(panel, PANEL_TEXT,
        PANEL_LABEL_STRING, "Error:",
        PANEL_VALUE_DISPLAY_LENGTH, 10,
        PANEL_VALUE_STORED_LENGTH, 20,
        PANEL_VALUE, error_string,0);
     length_item    = xv_create(panel, PANEL_TEXT,
        PANEL_LABEL_STRING, "Length: ",
        PANEL_VALUE_DISPLAY_LENGTH, 5,
        PANEL_VALUE_STORED_LENGTH, 20,
        PANEL_VALUE, length_string,0);
   iter_item 	= xv_create(panel, PANEL_TEXT,
        PANEL_LABEL_STRING, "Iterate: ",
        PANEL_VALUE_DISPLAY_LENGTH, 4,
        PANEL_VALUE_STORED_LENGTH, 20,
        PANEL_VALUE, iter_string,0);

    x_points = xv_create(panel, PANEL_SLIDER,
	    XV_X,xv_col(panel,0),
	    XV_Y,xv_row(panel,2),
            PANEL_LABEL_STRING, 	"X grid points    : ",
            PANEL_SLIDER_WIDTH,		 31,
	    PANEL_MIN_VALUE,		 16,
	    PANEL_MAX_VALUE, 		(MAX_X_GRID-2),0);

    y_points =  xv_create(panel, PANEL_SLIDER,
	    XV_X,xv_col(panel,0),
	    XV_Y,xv_row(panel,3),
	    PANEL_LABEL_STRING,		"Y grid points    :  ",
	    PANEL_SLIDER_WIDTH, 	241,
	    PANEL_MIN_VALUE, 		8,
	    PANEL_MAX_VALUE, 		(MAX_Y_GRID-1), 0);

    cycles_slider = xv_create(panel, PANEL_SLIDER,
	    XV_X,xv_col(panel,0),
	    XV_Y,xv_row(panel,4),
            PANEL_LABEL_STRING, 	"Number iterations:",
            PANEL_SLIDER_WIDTH,		 200,
	    PANEL_MIN_VALUE,		 1,
	    PANEL_MAX_VALUE, 		30000, 0);
    canvas_auto_clear = xv_create(panel, PANEL_CYCLE,
	    XV_X,xv_col(panel,0),
	    XV_Y,xv_row(panel,5)+5,
	    PANEL_LABEL_STRING, "Auto Clear: ",
	    PANEL_CHOICE_STRINGS, "NO", "YES", 0,
	    PANEL_VALUE, 1,  0);
   clear_data= xv_create(panel, PANEL_CYCLE,
	    XV_X,xv_col(panel,0),
	    XV_Y,xv_row(panel,6)+7,
	    PANEL_LABEL_STRING, "Auto zero array: ",
	    PANEL_CHOICE_STRINGS, "NO", "YES", 0,
	    PANEL_VALUE, 1,  0);

/*    method_item = xv_create(panel, PANEL_CYCLE,
	    XV_X,xv_col(panel,30),
	    XV_Y,xv_row(panel,5)+5,
	    PANEL_LABEL_STRING, "Method: ",
	    PANEL_CHOICE_STRINGS, "Jacobi", "Gauss-Siedel",
		"SOR", 0,
	    PANEL_VALUE, 0,  0);*/
   free_item 	= xv_create(panel, PANEL_TEXT,
 	    XV_X,xv_col(panel,30),
	    XV_Y,xv_row(panel,6)+7 ,
       PANEL_LABEL_STRING, "Param: ",
        PANEL_VALUE_DISPLAY_LENGTH, 8,
        PANEL_VALUE_STORED_LENGTH, 20,
        PANEL_VALUE, free_string,0);
   file_item 	= xv_create(panel, PANEL_TEXT,
 	    XV_X,xv_col(panel,0),
	    XV_Y,xv_row(panel,7)+9 ,
       PANEL_LABEL_STRING, "File: ",
        PANEL_VALUE_DISPLAY_LENGTH, 25,
        PANEL_VALUE_STORED_LENGTH, 40,
        PANEL_VALUE, file_string,0);

    xv_create(panel, PANEL_BUTTON, 
	    XV_X,xv_col(panel,71),
	    XV_Y,xv_row(panel,2)+2 ,
	    PANEL_LABEL_IMAGE, panel_button_image(panel, "CLEAR", 5, 0),
	    PANEL_NOTIFY_PROC, clear_canvas,0);

    xv_create(panel, PANEL_BUTTON, 
	    XV_X,xv_col(panel,64),
	    XV_Y,xv_row(panel,2) + 2,
	    PANEL_LABEL_IMAGE, panel_button_image(panel, "DRAW", 5, 0),
	    PANEL_NOTIFY_PROC, draw_canvas,0);

    xv_create(panel, PANEL_BUTTON, 
	XV_X,xv_col(panel,64),
	XV_Y,xv_row(panel,3) + 4,
	PANEL_LABEL_IMAGE, panel_button_image(panel, "QUIT", 5, 0),
	PANEL_NOTIFY_PROC, quit,0);
    xv_create(panel, PANEL_BUTTON, 
	XV_X,xv_col(panel,64),
	XV_Y,xv_row(panel,4) + 6,
	PANEL_LABEL_IMAGE, panel_button_image(panel, "HELP", 5, 0),
	PANEL_NOTIFY_PROC, help_msg,0);
    xv_create(panel, PANEL_BUTTON, 
	XV_X,xv_col(panel,64),
	XV_Y,xv_row(panel,5) + 8,
	PANEL_LABEL_IMAGE, panel_button_image(panel, "RUN", 5, 0),
	PANEL_NOTIFY_PROC, plot_routine,0);

    xv_create(panel, PANEL_BUTTON, 
	XV_X,xv_col(panel,71),
	XV_Y,xv_row(panel,4) + 6,
	PANEL_LABEL_IMAGE, panel_button_image(panel, "READ", 5, 0),
	PANEL_NOTIFY_PROC, read_coordinates,0);
    xv_create(panel, PANEL_BUTTON, 
	XV_X,xv_col(panel,71),
	XV_Y,xv_row(panel,5) + 8,
	PANEL_LABEL_IMAGE, panel_button_image(panel, "FILE", 5, 0),
	PANEL_NOTIFY_PROC, file_coordinates,0);

    window_fit_height(panel);
}
/* panel notify procs */

static void modify_canvas(item, event)
Panel_item	item;
Event		*event;
{ Pixwin *pw=   canvas_pixwin(canvas);
	window_set(canvas, 
	CANVAS_AUTO_CLEAR,	panel_get_value(canvas_auto_clear),0);
	read_conditions();
	pw_batch_on(pw);
	proc_draw_coordinates(pw);
	pw_batch_off(pw);
}

static void clear_canvas(item, event)
Panel_item	item;
Event		*event;
{
    Pixwin	*pw	= canvas_pixwin(canvas);
   pw_writebackground(pw, 0, 0, window_get(canvas, CANVAS_WIDTH), 
	window_get(canvas, CANVAS_HEIGHT), PIX_CLR);
}

static void draw_canvas(item, event)
Panel_item	item;
Event		*event;
{    Pixwin	*pw	= canvas_pixwin(canvas);
	width= (int)	window_get(canvas,CANVAS_WIDTH);
	height=(int)	window_get(canvas,CANVAS_HEIGHT);
	ystep=0.5;
	xstep=4.;
	frame_plot_area(pw,x_off,y_off);
}

static void quit(item, event)
Panel_item	item;
Event		*event;
{      window_set(base_frame, FRAME_NO_CONFIRM, TRUE, 0);
    window_destroy(base_frame);
}
static void plot_routine(item,event)
Panel_item	item;
Event 		event;
{
int k,ns;
double omega1,difference;
Pixwin	*pw	= canvas_pixwin(canvas);
	read_conditions();
	for(k=1;k<=n_iterates;k++)
	{
		pw_batch_on(pw);
		sprintf(iter_string,"%4d",k);
		panel_set_value(iter_item,iter_string);
		if((int)panel_get_value(canvas_auto_clear)==1)
			{clear_canvas();
			draw_canvas();
			}
		relax(&difference);
		sprintf(error_string,"%8.3le",difference);
		panel_set_value(error_item,error_string);
		proc_draw_coordinates(pw);
		pw_batch_off(pw);
	}
}
static void read_conditions()
{Pixwin	*pw	= canvas_pixwin(canvas);
	width= (int)	window_get(canvas,CANVAS_WIDTH);
	height=(int)	window_get(canvas,CANVAS_HEIGHT);
	m	= (int) 	panel_get_value(x_points);
	n	= (int) 	panel_get_value(y_points);
	n_iterates= (int) 	panel_get_value(cycles_slider);
	strcpy(top_depth_string,(char *)panel_get_value(top_depth_item));
	strcpy(bot_depth_string,(char *)panel_get_value(bot_depth_item));
	strcpy(hump_length_string,(char *)panel_get_value(hump_length_item));
	strcpy(flat_length_string,(char *)panel_get_value(flat_length_item));
	strcpy(phase_string,(char *)	panel_get_value(phase_item));
	strcpy(length_string,(char *)	panel_get_value(length_item));
	strcpy(omega_string,(char *)	panel_get_value(omega_item));
	strcpy(free_string,(char *)	panel_get_value(free_item));
	sscanf(top_depth_string,"%lf",&top_depth);
	sscanf(bot_depth_string,"%lf",&bot_depth);
	sscanf(hump_length_string,"%lf",&hump_length);
	sscanf(flat_length_string,"%lf",&flat_length);
	sscanf(phase_string,"%lf",&phase);
	sscanf(length_string,"%lf",&length);
	sscanf(omega_string,"%lf",&omega);
	sscanf(free_string,"%lf",&free_value);
	xmin=0.;
	xmax=length;
	if(bot_depth>0)
		{ymin=0.-bot_depth;}
	else
		{ymin=0.;}
	ymax=1.+top_depth;
	intialise_pq();
	if((int)panel_get_value(clear_data)==1)
		{
		pw_batch_on(pw);
		clear_canvas();
		draw_canvas();
		initialise_walls();
		initialise_x_y();
		proc_draw_coordinates(pw);
		pw_batch_off(pw);
		}

}
static double absol(x)
double x;
{double y;
if(x>0.){return(x);}
else{y=0.-x;return(y);}
}
static int x_pix(x,xm,w,ox)
double x,xm;
int w,ox;
{return(ox+(w-2*ox)*x/xm);}
static int y_pix(y,y0,y1,h,oy)
double y,y0,y1;
int h,oy;
{
	return(oy+(h-2*oy)*(1-(y-y0)/(y1-y0)));
}
static void draw_line(pw, x1, y1, x2, y2)
Pixwin	*pw;
int	x1, y1, x2, y2;
{
	pw_vector(pw, x1, y1, x2, y2, PIX_SRC, 1);
}
static void frame_plot_area(pw,ox,oy)
Pixwin	*pw;int ox,oy;
{int i0,j0,i,j,k,w,h;
double x,y,dx,dy;char str[20];
	w=width;
	h=height;
	i0=ox;
	j0=h-oy;
	i=w-ox;
	j=j0;
	draw_line(pw,i0,j0,i,j);
	i0=i;
	j0=j;
	j=oy;
	draw_line(pw,i0,j0,i,j);
	j0=j;
	i=ox;
	draw_line(pw,i0,j0,i,j);
	i0=ox;
	j0=oy;
	j=h-oy;
	draw_line(pw,i0,j0,i,j);
	i0=w/4;
	j0=h-5;
	pw_text(pw,i0,j0,PIX_SRC,NULL,xlabel);
	i0=10;	
	j0=20;
	pw_text(pw,i0,j0,PIX_SRC,NULL,title);
	for(x=xmin;x<=xmax;x+=xstep)
		{i0=x_pix(x,xmax,w,ox);
		j0=h-oy+5;
		j=h-oy;
		draw_line(pw,i0,j0,i0,j);
		sprintf(str,"%3.1f",x);
		pw_text(pw,i0-10,j0+10,PIX_SRC,NULL,str);
		}
	for(y=ymin;y<=ymax;y+=ystep)
		{j0=y_pix(y,ymin,ymax,h,oy);
		i0=ox-5;
		i=ox;
		draw_line(pw,i0,j0,i,j0);
		sprintf(str,"%3.1f",y);
		pw_text(pw,i0-30,j0+5,PIX_SRC,NULL,str);
		}
}
static void help_msg(item, event)
Panel_item	item;Event		*event;
{static char *msg="Void at present";
}

static void proc_draw_coordinates(pw)
Pixwin *pw;
{	int i,j,i0,j0,i1,j1;
	for(i=0;i<=m;i++)
	{for(j=0;j<n;j++)
		{i0=x_pix(x[j][i],xmax,width,x_off);
		j0=y_pix(y[j][i],ymin,ymax,height,y_off);
		i1=x_pix(x[j+1][i],xmax,width,x_off);
		j1=y_pix(y[j+1][i],ymin,ymax,height,y_off);
		draw_line(pw,i0,j0,i1,j1);
		}
	}
	for(j=0;j<=n;j++)
		{
		for(i=0;i<m;i++)
		{i0=x_pix(x[j][i],xmax,width,x_off);
		j0=y_pix(y[j][i],ymin,ymax,height,y_off);
		i1=x_pix(x[j][i+1],xmax,width,x_off);
		j1=y_pix(y[j][i+1],ymin,ymax,height,y_off);
		draw_line(pw,i0,j0,i1,j1);
		}
	}
}

static void evaluate_coefficients()
{
int i,j;
double xpi,xmi,ypi,ymi,xpj,xmj,ypj,ymj;
	for(j=1;j<n;j++)
	{
		for(i=1;i<=m;i++)
		{xpi=x[j][i+1];
		xmi=x[j][i-1];
		xpj=x[j+1][i];
		xmj=x[j-1][i];
		ypi=y[j][i+1];
		ymi=y[j][i-1];
		ypj=y[j+1][i];
		ymj=y[j-1][i];
		alpha[j][i]=((xpj-xmj)*(xpj-xmj)+(ypj-ymj)*(ypj-ymj))/4;
		beeta[j][i]=((xpi-xmi)*(xpj-xmj)+(ypi-ymi)*(ypj-ymj))/4;
		gama[j][i]=((xpi-xmi)*(xpi-xmi)+(ypi-ymi)*(ypi-ymi))/4;
		jacobian[j][i]=((xpi-xmi)*(ypj-ymj)-(xpj-xmj)*(ypi-ymi))/4;
		}
	}
}
static void relax(difference)
double *difference;
{	int i,j;
	double omega1,j2,a,b,g,new,dif;
	omega1=1.-omega;
	dif=0.;
	evaluate_coefficients();
	calc_walls();
/*	for(j=1;j<n;j++)
	{
		for(i=1;i<=m;i++)
		{
		j2=jacobian[j][i]*jacobian[j][i];
		a=alpha[j][i];
		b=beeta[j][i];
		g=gama[j][i];
		new=omega1*x[j][i]+omega*(
		   (a+j2*p[j][i]/2)*x[j][i+1]+(a-j2*p[j][i]/2)*x[j][i-1]
		   +(g+j2*q[j][i]/2)*x[j+1][i]+(g-j2*q[j][i]/2)*x[j-1][i]
		   -b*(x[j+1][i+1]-x[j-1][i+1]-x[j+1][i-1]+x[j-1][i-1])/2
		       )/(2*(a+g));
		dif=dif+absol(new-x[j][i]);
		x[j][i]=new;
		new=omega1*y[j][i]+omega*(
		   (a+j2*p[j][i]/2)*y[j][i+1]+(a-j2*p[j][i]/2)*y[j][i-1]
		   +(g+j2*q[j][i]/2)*y[j+1][i]+(g-j2*q[j][i]/2)*y[j-1][i]
		   -b*(y[j+1][i+1]-y[j-1][i+1]-y[j+1][i-1]+y[j-1][i-1])/2
		       )/(2*(a+g));
		dif=dif+absol(new-y[j][i]);
		y[j][i]=new;
		}
		x[j][0]=x[j][m]-length;
		x[j][m+1]=x[j][1]+length;
		y[j][0]=y[j][m];
		y[j][m+1]=y[j][1];
	}
*/
	for(j=1;j<n;j++)
	{
		for(i=1;i<=m;i++)
		{
		j2=jacobian[j][i]*jacobian[j][i];
		a=alpha[j][i];
		b=beeta[j][i];
		g=gama[j][i];
		new=omega1*y[j][i]+omega*(
		   (a+j2*p[j][i]/2)*y[j][i+1]+(a-j2*p[j][i]/2)*y[j][i-1]
		   +(g+j2*q[j][i]/2)*y[j+1][i]+(g-j2*q[j][i]/2)*y[j-1][i]
		   -b*(y[j+1][i+1]-y[j-1][i+1]-y[j+1][i-1]+y[j-1][i-1])/2
		       )/(2*(a+g));
		dif=dif+absol(new-y[j][i]);
		y[j][i]=new;
		}
		y[j][0]=y[j][m];
		y[j][m+1]=y[j][1];
	}
	evaluate_coefficients();
	for(j=1;j<n;j++)
	{
		for(i=1;i<=m;i++)
		{
		j2=jacobian[j][i]*jacobian[j][i];
		a=alpha[j][i];
		b=beeta[j][i];
		g=gama[j][i];
		new=omega1*x[j][i]+omega*(
		   (a+j2*p[j][i]/2)*x[j][i+1]+(a-j2*p[j][i]/2)*x[j][i-1]
		   +(g+j2*q[j][i]/2)*x[j+1][i]+(g-j2*q[j][i]/2)*x[j-1][i]
		   -b*(x[j+1][i+1]-x[j-1][i+1]-x[j+1][i-1]+x[j-1][i-1])/2
		       )/(2*(a+g));
		dif=dif+absol(new-x[j][i]);
		x[j][i]=new;
		}
		x[j][0]=x[j][m]-length;
		x[j][m+1]=x[j][1]+length;
	}

	(*difference)=dif/(m*n);
}
calc_walls()
{
int i,j;
double la,dx,w;
	for(i=1;i<=m-1;i++)
	{
		la=(y[n][i+1]-y[n][i-1])/(x[n][i+1]-x[n][i-1]);
		x[n][i]=x[n-1][i]+la*(y[n-1][i]-y[n][i]);
		top_wall_value(x[n][i],&w);
		y[n][i]=w;

		la=(y[0][i+1]-y[0][i-1])/(x[0][i+1]-x[0][i-1]);
		x[0][i]=x[1][i]+la*(y[1][i]-y[0][i]);
		bot_wall_value(x[0][i],&w);
		y[0][i]=w;
	}
		x[n][0]=0.;
		top_wall_value(x[n][0],&w);
		y[n][0]=w;
		x[0][0]=0.;
		bot_wall_value(x[0][0],&w);
		y[0][0]=w;
		x[n][m]=x[n][0]+length;
		y[n][m]=y[n][0];
		x[0][m]=x[0][0]+length;
		y[0][m]=y[0][0];
		x[n][m+1]=x[n][1]+length;
		y[n][m+1]=y[n][1];
		x[0][m+1]=x[0][1]+length;
		y[0][m+1]=y[0][1];
	
}

static void initialise_x_y()
{int i,j;double dx,dy;
	dx=length/m;
	dy=1./n;
	for(i=0;i<=(m+1);i++)
	{
		for(j=1;j<n;j++)
		{
		x[j][i]=i*dx;
		y[j][i]=(y[n][i]-y[0][i])*j*dy+y[0][i];
		}
	}
}
static void initialise_walls()
{int i,j;double dx,eta,f,t,dz,z,w,arc,c_arc,darc,sqrt();
	dx=length/m;
	for(i=0;i<=(m+1);i++)
	{
	x[n][i]=i*dx;
	top_wall_value(x[n][i],&w);
	y[n][i]=w;
	x[0][i]=x[n][i];
	bot_wall_value(x[0][i],&w);
	y[0][i]=w;
	}
	intialise_pq();
}
intialise_pq()
{int i,j;double eta,t,f;
	for(i=0;i<m+2;i++)
		{for(j=0;j<=n;j++)
			{p[j][i]=0.;
			eta=(2.*j-n)/n;
			t=free_value*PI*PI*sin(PI*eta);
			f=1+free_value*PI*cos(PI*eta);

			q[j][i]=n*t/(2.*f*f*f);
			}
		}
}
top_wall_value(z,w)
double z,*w;
{
	double f1,f2,l1,l2,l3;
	
	f1=(1-cos(PI*(z-flat_length)/hump_length))/2;
	f2=(1-cos(PI*(z-length+flat_length)/hump_length))/2;
	if(z<=flat_length)				{(*w)=1.;}
	else if(z<=flat_length+hump_length)		{(*w)= 1.+top_depth*f1;}
	else if(z<=(length -hump_length-flat_length))	{(*w)=1.+top_depth;}
	else if(z<=(length -flat_length))		{(*w)= 1.+top_depth*f2;}
	else						{(*w)=1.;}

}
bot_wall_value(z,w)
double z,*w;
{
	double f1,f2;
	f1=(1-cos(PI*(z-flat_length)/hump_length))/2;
	f2=(1-cos(PI*(z-length+flat_length)/hump_length))/2;
	if(z<=flat_length)				{(*w)=0.;}
	else if(z<=flat_length+hump_length)		{(*w)= 0.-bot_depth*f1;}
	else if(z<=(length -hump_length-flat_length))	{(*w)=0.-bot_depth;}
	else if(z<=(length -flat_length))		{(*w)= 0.-bot_depth*f2;}
	else						{(*w)=0.;}

}
static void file_coordinates()
{
	FILE *fp;
	int i,j,k;
	double f;
	k=n/2;
	strcpy(file_string,(char *)	panel_get_value(file_item));

	if((fp=fopen(file_string,"w"))==NULL)
             {printf("can't open file : %s\n",file_string);}
	else 
             {/* symmetrise data */
/*		for(i=0;i<=m+1;i++)
		    {if(n-2*k==0)*/	/*n even: centerline there*/ /*
			{y[k][i]=0.;k=k-1;}
		    for(j=1;j<=k;j++)
			{f=0.5*(y[j][i]-y[n-j][i]);
			y[j][i]=f;
			y[n-j][i]=0.-f;
			f=0.5*(x[j][i]+x[n-j][i]);
			x[j][i]=f;
			x[n-j][i]=f;
			}
		    }*/
	fprintf(fp,"%6d %6d %6.2f %6.2f %6.2f\n",
			m,n,length,top_depth,bot_depth);
	fprintf(fp,"%6.2f %6.2f %6.2f\n",
			phase,flat_length,hump_length);
		for(j=0;j<=n;j++)
			{for(i=0;i<=m;i++)
			{fprintf(fp," %12.7f %12.7f %12.7f %12.7f\n",
				x[j][i],y[j][i],p[j][i],q[j][i]);
			}
			}
		fclose(fp);

		}
}
static void read_coordinates()
{
	FILE *fp;
	int i,j,k;
	double f,g,p1,q1;
	strcpy(file_string,(char *)	panel_get_value(file_item));

	if((fp=fopen(file_string,"r"))==NULL)
             {printf("can't open file : %s\n",file_string);}
	else 
             {	
	fscanf(fp,"%d %d %lf %lf %lf",
			&m,&n,&length,&top_depth,&bot_depth);
	fscanf(fp,"%lf %lf %lf",&phase,&flat_length,&hump_length);
	sprintf(length_string,"%6.2f",length);
	sprintf(top_depth_string,"%6.2f",top_depth);
	sprintf(bot_depth_string,"%6.2f",bot_depth);
	sprintf(phase_string,"%6.0f",phase);
	sprintf(flat_length_string,"%6.2f",flat_length);
	sprintf(hump_length_string,"%6.2f",hump_length);

	panel_set_value(x_points,m);
	panel_set_value(y_points,n);
	panel_set_value(length_item,length_string);
	panel_set_value(top_depth_item,top_depth_string);
	panel_set_value(bot_depth_item,bot_depth_string);
	panel_set_value(phase_item,phase_string);
	panel_set_value(flat_length_item,flat_length_string);
	panel_set_value(hump_length_item,hump_length_string);
	panel_set_value(clear_data,0);
		for(j=0;j<=n;j++)
			{for(i=0;i<=m;i++)
			{fscanf(fp," %lf %lf %lf %lf",&f,&g,&p1,&q1);
			x[j][i]=f;y[j][i]=g;
			p[j][i]=p1;q[j][i]=q1;
			}
			}
	fclose(fp);
		
		}
}