#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 MAX_Y_GRID	65
#define MAX_X_GRID	514
#define TWOPI		(6.28318)


double 	x[MAX_Y_GRID][MAX_X_GRID],y[MAX_Y_GRID][MAX_X_GRID];
double 	u[MAX_Y_GRID][MAX_X_GRID],v[MAX_Y_GRID][MAX_X_GRID];
double psi[MAX_Y_GRID][MAX_X_GRID],vort[MAX_Y_GRID][MAX_X_GRID];
double old_vort[MAX_Y_GRID][MAX_X_GRID],new_vort[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 	z1[MAX_Y_GRID][MAX_X_GRID],r1[MAX_Y_GRID][MAX_X_GRID],
 	z2[MAX_Y_GRID][MAX_X_GRID],r2[MAX_Y_GRID][MAX_X_GRID],
 	z3[MAX_Y_GRID][MAX_X_GRID],r3[MAX_Y_GRID][MAX_X_GRID],
 	z4[MAX_Y_GRID][MAX_X_GRID],r4[MAX_Y_GRID][MAX_X_GRID];
double 	time_residual,
	residual,
	tolerance={.0005};
int	ix1,iy1,
	ix2,iy2,
	ix3,iy3,
	ix4,iy4;


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

static Frame		base_frame;
static Canvas		canvas;

static Panel_item	x_points,y_points,cycles_slider;
static Panel_item	canvas_auto_clear,
			clear_data,
			iter_item,
			method_item,
			name_item,
			top_depth_item,
			bot_depth_item,
			phase_item,
			omega_item,
			error_item,
			length_item,
			time_step_item,
			reynolds_item,
			strouhal_item,
    			number_item ,
			print_steps_item,
			file_item,
			free_item,
  			asym_item,
    			mean_item,
			con1_item,
			con2_item ,
			con3_item,
			con4_item,
			con5_item,
			con6_item,
			con7_item,
			con8_item	;


int n_h,kcon;

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 eps,omega,k_i,k_o,gam,g1,g1og,period,wmax;
double con[8]={0.,.5,-.5,1.,-1.,1.15,-1.15,1.3};
char 	top_depth_string[20]		={"1.0"},
 	bot_depth_string[20]		={"1.0"},
	phase_string[20]		={"0.0"},
	omega_string[20]		={"1.0"},
	reynolds_string[20]		={"20."},
	strouhal_string[20]		={".02"},
	time_step_string[20]		={"1.0"},
	error_string[20]		={"0.0005"},
	length_string[20]		={"8.0"},
	free_string[20]			={"0.001"},
	asym_string[20]			={"0.0"},
	mean_string[20]			={"0.0"},
	file_string[40]			={"../"},
	iter_string[20]			={"0"};

char 	number_str[20]		={"8"},
	con1_str[20]		={"0."},
	con2_str[20]		={"0.5"},
	con3_str[20]		={"-0.5"},
	con4_str[20]		={"1."},
	con5_str[20]		={"-1."},
	con6_str[20]		={"1.15"},
	con7_str[20]		={"-1.15"},
	con8_str[20]		={"1.3"};
double length,top_depth,bot_depth,phase,
	flat_length,hump_length,free_value={.001};
double reynolds={20.},
	strouhal={.02},
	asym={0.0},
	mean={0.0},
	pe,
	courant,
	dt={.001},
	current_time;
int number={8};
int width,height;
int x_off={40};
int y_off={40};
int xx[2],yy[2];
int print_steps,steps;
	
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();
static void     modify_contours();
static void     enter_contours();
static void 	contours_done();
/* canvs procs */
/* utilities */
char 		*get_name();
static void	proc_contour_data();
static void 	evaluate_coefficients();
static void 	read_coordinates();
static void 	read_pw();
static void 	write_pw();
static void	init_panel();
static void	draw_line();
static void 	plot_routine();
static void 	frame_plot_area();
static int 	x_pix();
static int 	y_pix(),display_sol_active={0},
dimensions_active={0};
static void	display_sol();
static double 	absol();

FILE *dumpf,*dumps;
FILE *dumpu1,*dumpu2,*dumpu3;
FILE *logfile;
main(argc, argv)
int argc;
char **argv;
{int i,j;
	logfile=fopen("steady.log","w");	
    xv_init(XV_INIT_ARGC_PTR_ARGV,&argc,argv,NULL);
	
    base_frame = xv_create(0, FRAME,
	    FRAME_LABEL,    "Navier-Stokes Solver - Central Difference Version",
	    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);
             


    window_main_loop(base_frame);
	fclose(logfile);
    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, (s=(char *)get_name()) ? s : "",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);

    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);
printf("6");
    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, "Time: ",
        PANEL_VALUE_DISPLAY_LENGTH, 12,
        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,		 200,
	    PANEL_MIN_VALUE,		 5,
	    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, 	200,
	    PANEL_MIN_VALUE, 		5,
	    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);
    print_steps_item = xv_create(panel, PANEL_SLIDER,
	    XV_X,xv_col(panel,60),
	    XV_Y,xv_row(panel,1),
            PANEL_LABEL_STRING, 	"Print steps:",
            PANEL_SLIDER_WIDTH,		 49,
	    PANEL_MIN_VALUE,		 1,
	    PANEL_MAX_VALUE, 		 50, 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);
   reynolds_item 	= xv_create(panel, PANEL_TEXT,
 	    XV_X,xv_col(panel,30),
	    XV_Y,xv_row(panel,5)+5 ,
       PANEL_LABEL_STRING, "Re: ",
        PANEL_VALUE_DISPLAY_LENGTH, 8,
        PANEL_VALUE_STORED_LENGTH, 20,
        PANEL_VALUE, reynolds_string,0);
   strouhal_item 	= xv_create(panel, PANEL_TEXT,
 	    XV_X,xv_col(panel,45),
	    XV_Y,xv_row(panel,5)+5 ,
       PANEL_LABEL_STRING, "St: ",
        PANEL_VALUE_DISPLAY_LENGTH, 8,
        PANEL_VALUE_STORED_LENGTH, 20,
        PANEL_VALUE, strouhal_string,0);

  free_item 	= xv_create(panel, PANEL_TEXT,
 	    XV_X,xv_col(panel,25),
	    XV_Y,xv_row(panel,6)+7 ,
       PANEL_LABEL_STRING, "time step: ",
        PANEL_VALUE_DISPLAY_LENGTH, 8,
        PANEL_VALUE_STORED_LENGTH, 20,
        PANEL_VALUE, free_string,0);
   asym_item 	= xv_create(panel, PANEL_TEXT,
 	    XV_X,xv_col(panel,43),
	    XV_Y,xv_row(panel,6)+7 ,
       PANEL_LABEL_STRING, "Symmetry: ",
        PANEL_VALUE_DISPLAY_LENGTH, 12,
        PANEL_VALUE_STORED_LENGTH, 20,
        PANEL_VALUE, asym_string,0);
   file_item 	= xv_create(panel, PANEL_TEXT,
 	    XV_X,xv_col(panel,0),
	    XV_Y,xv_row(panel,7) ,
       PANEL_LABEL_STRING, "File: ",
        PANEL_VALUE_DISPLAY_LENGTH, 25,
        PANEL_VALUE_STORED_LENGTH, 40,
        PANEL_VALUE, file_string,0);
   mean_item 	= xv_create(panel, PANEL_TEXT,
 	    XV_X,xv_col(panel,45),
	    XV_Y,xv_row(panel,7) ,
       PANEL_LABEL_STRING, "Mean: ",
        PANEL_VALUE_DISPLAY_LENGTH, 8,
        PANEL_VALUE_STORED_LENGTH, 20,
        PANEL_VALUE, mean_string,0);

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

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

    xv_create(panel, PANEL_BUTTON, 
	XV_X,xv_col(panel,64),
	XV_Y,xv_row(panel,3),
	PANEL_LABEL_IMAGE, panel_button_image(panel, "QUIT", 4, 0),
	PANEL_NOTIFY_PROC, quit,0);
    xv_create(panel, PANEL_BUTTON, 
	XV_X,xv_col(panel,71),
	XV_Y,xv_row(panel,3),
	PANEL_LABEL_IMAGE, panel_button_image(panel, "READ PW", 7, 0),
	PANEL_NOTIFY_PROC, read_pw,0);
    xv_create(panel, PANEL_BUTTON, 
	XV_X,xv_col(panel,64),
	XV_Y,xv_row(panel,4),
	PANEL_LABEL_IMAGE, panel_button_image(panel, "HELP", 4, 0),
	PANEL_NOTIFY_PROC, help_msg,0);
    xv_create(panel, PANEL_BUTTON, 
	XV_X,xv_col(panel,64),
	XV_Y,xv_row(panel,5),
	PANEL_LABEL_IMAGE, panel_button_image(panel, "RUN", 3, 0),
	PANEL_NOTIFY_PROC, plot_routine,0);

    xv_create(panel, PANEL_BUTTON, 
	XV_X,xv_col(panel,71),
	XV_Y,xv_row(panel,4),
	PANEL_LABEL_IMAGE, panel_button_image(panel, "READ XY", 7, 0),
	PANEL_NOTIFY_PROC, read_coordinates,0);
    xv_create(panel, PANEL_BUTTON, 
	XV_X,xv_col(panel,71),
	XV_Y,xv_row(panel,6),
	PANEL_LABEL_IMAGE, panel_button_image(panel, "CONTOURS", 8, 0),
	PANEL_NOTIFY_PROC, modify_contours,0);
    xv_create(panel, PANEL_BUTTON, 
	XV_X,xv_col(panel,71),
	XV_Y,xv_row(panel,5),
	PANEL_LABEL_IMAGE, panel_button_image(panel, "SAVE PW", 7, 0),
	PANEL_NOTIFY_PROC, write_pw,0);

    window_fit_height(panel);
}

/* panel notify procs */

static void modify_canvas(item, event)
Panel_item	item;
Event		*event;
{    window_set(canvas, 
	CANVAS_AUTO_CLEAR,	panel_get_value(canvas_auto_clear),0);
/*if(((int) panel_get_value(canvas_auto_clear))==TRUE)
	clear_canvas(item,event); */
}

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);
	frame_plot_area(pw,x_off,y_off);
}

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);
	frame_plot_area(pw,x_off,y_off);
}

static void quit(item, event)
Panel_item	item;
Event		*event;
{    /* quit without user confirmation */
    window_set(base_frame, FRAME_NO_CONFIRM, TRUE, 0);
    window_destroy(base_frame);
}

/* canvas procs */
	        
/* utilities */
static void plot_routine(item,event)
Panel_item	item;
Event		*event;
{
Pixwin	*pw	= canvas_pixwin(canvas);
	int i,j,k,n,m,ns,method;
	int u16,u32,u48;
	double omega1,dx,dy,f1,f2,differ,b,xymax,transfer;
	double v1,v2,v3;

	width= (int)	window_get(canvas,CANVAS_WIDTH);
	height=(int)	window_get(canvas,CANVAS_HEIGHT);
	ns= (int) 	panel_get_value(cycles_slider);
	print_steps=(int)panel_get_value(print_steps_item);
	strcpy(reynolds_string,(char *)	panel_get_value(reynolds_item));
	strcpy(strouhal_string,(char *)	panel_get_value(strouhal_item));
	strcpy(free_string,(char *)	panel_get_value(free_item));
	strcpy(error_string,(char *)	panel_get_value(error_item));
	strcpy(mean_string,(char *)	panel_get_value(mean_item));
	strcpy(omega_string,(char *)	panel_get_value(omega_item));
	sscanf(reynolds_string,"%lf",&reynolds);
	sscanf(strouhal_string,"%lf",&strouhal);
	sscanf(free_string,"%lf",&free_value);
	sscanf(error_string,"%lf",&tolerance);
	sscanf(mean_string,"%lf",&mean);
	sscanf(omega_string,"%lf",&omega);
	dt=free_value;
	courant=dt/strouhal;
	pe=dt/(reynolds*strouhal);
        printf("Reynolds number:%8.1lf  Strouhal number:%8.4lf\n",
			reynolds,strouhal);
        printf(" c= %12.4lf pe = %12.4lf\n",courant,pe);
	u16=ix1/4;
	u32=ix1/2;
	u48=3*ix1/4;
	if((int)panel_get_value(clear_data)==1)
		{
		current_time=dt;
		initialise_all();
		}
	k=0;
	time_residual=0;
	while((k<ns)&&(time_residual<100.))
		{	
		for(steps=0;steps<print_steps;steps++)
			{
			current_time=current_time+dt;
			step_vorticity();
			}
		fprintf(logfile,"%12.5lf %12.5le\n",current_time,vort[0][ix1/2]);
		pw_batch_on(pw);
		calculate_asymmetry();
		sprintf(asym_string,"%8.3e",asym);
		printf("%7.4lf %6.3lf %8.3lf  %8.3lf %8.3le\n",
			current_time,psi[iy1][0],wmax,asym,time_residual);
		panel_set_value(asym_item,asym_string);
		sprintf(iter_string,"%8.4f",current_time);
		panel_set_value(iter_item,iter_string);
		if((int)panel_get_value(canvas_auto_clear)==1)
			{clear_canvas();
			draw_canvas();
			}
		frame_plot_area(pw,x_off,y_off);
		draw_walls(pw);
		for(i=0;i<number;i++)
			{f1=psi[iy1][0]*con[i];
 			proc_contour_data(pw,f1);
			}
		pw_batch_off(pw);
	        k++;
		}
}
initialise_all()
{
int i,j;double f,g,z;
	for(i=0;i<=ix1+1;i++)
		{
		for(j=0;j<=iy1;j++)
		{
			z=j*(1./iy1);
			psi[j][i]=(6*z*z-4*z*z*z-1)/24.;
			vort[j][i]=-(1.-2*z)/2;
			old_vort[j][i]=-(1.-2*z)/2;
			new_vort[j][i]=-(1.-2*z)/2;
		}
		}
	calculate_velocities();
/*	calculate_boundary_vorticity();*/
	for(i=0;i<=ix1+1;i++)
		{vort[0][i]=new_vort[0][i];
		vort[iy1][i]=new_vort[iy1][i];
		}
	

}
static double absol(x)
double x;
{double y;
if(x>0.){return(x);}
else{y=0.-x;return(y);}
}
static void help_msg(item, event)
Panel_item	item;Event		*event;
{static char *msg="NAVIER-STOKES SIMULATOR";
	popup_msg(base_frame,msg,count_lines(msg));
}
static void popup_msg(frame, msg, lines)
Frame *frame;
char *msg;
int lines;
{
    init_popup_msg(frame, msg, lines);
    textsw_insert(popup_text, msg, strlen(msg));
    window_set(popup_frame,WIN_SHOW, TRUE,0);
}
static void init_popup_msg(baseframe, msg, lines)
Frame baseframe;
char *msg;
{
    popup_frame = xv_create(baseframe, FRAME,0);
    popup_panel = xv_create(popup_frame, PANEL,
         XV_X, xv_col(baseframe,0), 0);
    xv_create(popup_panel, PANEL_BUTTON,
        PANEL_LABEL_IMAGE,panel_button_image(popup_panel,"Done",4,NULL),
        PANEL_NOTIFY_PROC,popup_textsw_done,0);
    window_fit(popup_panel);
    popup_text = xv_create(popup_frame, TEXTSW,
        WIN_ERROR_MSG, "Can't create window.",
        WIN_ROWS, 30,
        WIN_COLUMNS, max_line(msg)+4,
        TEXTSW_IGNORE_LIMIT, TEXTSW_INFINITY,TEXTSW_CONTENTS, msg,
        TEXTSW_BROWSING, TRUE,TEXTSW_DISABLE_LOAD, TRUE,
        TEXTSW_DISABLE_CD, TRUE,0);
    window_fit(popup_frame);
}
static void popup_textsw_done()
{   window_set(popup_frame, FRAME_NO_CONFIRM, TRUE, 0);
    window_destroy(popup_frame);
    if(display_sol_active==1)display_sol_active=0;
}

static int count_lines(s)
char *s;
{int count = 0;while (*s) {if (*s++ == '\n')count++;}
    return count+1;
}
static int max_line(s)
char *s;
{int lmax = 0, count = 0;while (*s) 
	{if (*s++ == '\n') {
            	if (count > lmax)lmax = count;
            	count = 0;continue;}
        count++;}
    if (count > lmax)lmax = count;
    return (lmax);
}
char *get_name()
{char *s;struct passwd *p;
    if (s = (char *)getlogin())return s;
    p = (struct passwd *)getpwuid(getuid());
    return p->pw_name;
}
static void modify_contours(item, event)
Panel_item	item;Event		*event;
{ double pr_i;
  popup_frame = xv_create(base_frame, FRAME,
	    FRAME_LABEL,    "Contour Heights",
	    FRAME_SHOW_LABEL,TRUE,0);

    popup_panel = xv_create(popup_frame, PANEL,NULL);
    xv_create(popup_panel, PANEL_BUTTON, 
	XV_X,xv_col(popup_panel,1),
	XV_Y,xv_row(popup_panel,0) + 4,
	PANEL_LABEL_IMAGE, panel_button_image(popup_panel, "DONE", 5, 0),
	PANEL_NOTIFY_PROC, contours_done,0);
    xv_create(popup_panel, PANEL_BUTTON, 
	XV_X,xv_col(popup_panel,18),
	XV_Y,xv_row(popup_panel,0) + 4,
	PANEL_LABEL_IMAGE, panel_button_image(popup_panel, "ENTER", 5, 0),
	PANEL_NOTIFY_PROC, enter_contours,0);

    number_item 	= xv_create(popup_panel, PANEL_TEXT,
 	XV_X,xv_col(popup_panel,1),
	XV_Y,xv_row(popup_panel,2) ,
        PANEL_LABEL_STRING, "Number (1-8): ",
        PANEL_VALUE_DISPLAY_LENGTH, 8,
        PANEL_VALUE_STORED_LENGTH, 20,
        PANEL_VALUE, number_str,0);
    con1_item 	= xv_create(popup_panel, PANEL_TEXT,
 	XV_X,xv_col(popup_panel,1),
	XV_Y,xv_row(popup_panel,3)+2 ,
       PANEL_LABEL_STRING, "Height 1    : ",
        PANEL_VALUE_DISPLAY_LENGTH, 8,
        PANEL_VALUE_STORED_LENGTH, 20,
        PANEL_VALUE, con1_str,0);
    con2_item 	= xv_create(popup_panel, PANEL_TEXT,
 	XV_X,xv_col(popup_panel,1),
	XV_Y,xv_row(popup_panel,4) + 4,
       PANEL_LABEL_STRING, "Height 2    : ",
        PANEL_VALUE_DISPLAY_LENGTH, 8,
        PANEL_VALUE_STORED_LENGTH, 20,
        PANEL_VALUE, con2_str,0);
    con3_item 	= xv_create(popup_panel, PANEL_TEXT,
 	XV_X,xv_col(popup_panel,1),
	XV_Y,xv_row(popup_panel,5) + 6,
       PANEL_LABEL_STRING, "Height 3    : ",
        PANEL_VALUE_DISPLAY_LENGTH, 8,
        PANEL_VALUE_STORED_LENGTH, 20,
        PANEL_VALUE, con3_str,0);
    con4_item 	= xv_create(popup_panel, PANEL_TEXT,
 	XV_X,xv_col(popup_panel,1),
	XV_Y,xv_row(popup_panel,6) + 8,
       PANEL_LABEL_STRING, "Height 4    : ",
        PANEL_VALUE_DISPLAY_LENGTH, 8,
        PANEL_VALUE_STORED_LENGTH, 20,
        PANEL_VALUE, con4_str,0);
    con5_item 	= xv_create(popup_panel, PANEL_TEXT,
 	XV_X,xv_col(popup_panel,1),
	XV_Y,xv_row(popup_panel,7) + 8,
       PANEL_LABEL_STRING, "Height 5    : ",
        PANEL_VALUE_DISPLAY_LENGTH, 8,
        PANEL_VALUE_STORED_LENGTH, 20,
        PANEL_VALUE, con5_str,0);
    con6_item 	= xv_create(popup_panel, PANEL_TEXT,
 	XV_X,xv_col(popup_panel,1),
	XV_Y,xv_row(popup_panel,8) + 8,
       PANEL_LABEL_STRING, "Height 6    : ",
        PANEL_VALUE_DISPLAY_LENGTH, 8,
        PANEL_VALUE_STORED_LENGTH, 20,
        PANEL_VALUE, con6_str,0);
    con7_item 	= xv_create(popup_panel, PANEL_TEXT,
 	XV_X,xv_col(popup_panel,1),
	XV_Y,xv_row(popup_panel,9) + 8,
       PANEL_LABEL_STRING, "Height 7    : ",
        PANEL_VALUE_DISPLAY_LENGTH, 8,
        PANEL_VALUE_STORED_LENGTH, 20,
        PANEL_VALUE, con7_str,0);
    con8_item 	= xv_create(popup_panel, PANEL_TEXT,
 	XV_X,xv_col(popup_panel,1),
	XV_Y,xv_row(popup_panel,10) + 8,
       PANEL_LABEL_STRING, "Height 8    : ",
        PANEL_VALUE_DISPLAY_LENGTH, 8,
        PANEL_VALUE_STORED_LENGTH, 20,
        PANEL_VALUE, con8_str,0);
   window_fit(popup_panel);
   window_fit(popup_frame);
sprintf(number_str,"%8d",number);	panel_set_value(number_item,number_str);
sprintf(con1_str,"%8.3f",con[0]);	panel_set_value(con1_item,con1_str);
sprintf(con2_str,"%8.3f",con[1]);	panel_set_value(con2_item,con2_str);
sprintf(con3_str,"%8.3f",con[2]);	panel_set_value(con3_item,con3_str);
sprintf(con4_str,"%8.3f",con[3]);	panel_set_value(con4_item,con4_str);
sprintf(con5_str,"%8.3f",con[4]);	panel_set_value(con5_item,con5_str);
sprintf(con6_str,"%8.3f",con[5]);	panel_set_value(con6_item,con6_str);
sprintf(con7_str,"%8.3f",con[6]);	panel_set_value(con7_item,con7_str);
sprintf(con8_str,"%8.3f",con[7]);	panel_set_value(con8_item,con8_str);
window_set(popup_frame,WIN_SHOW, TRUE,0);
}
static void enter_contours()
{double pr_i;
  strcpy(number_str,(char *)panel_get_value(number_item));
  strcpy(con1_str,(char *)  panel_get_value(con1_item));
  strcpy(con2_str,(char *) panel_get_value(con2_item));
  strcpy(con3_str,(char *) panel_get_value(con3_item));
  strcpy(con4_str,(char *) panel_get_value(con4_item));
  strcpy(con5_str,(char *) panel_get_value(con5_item));
  strcpy(con6_str,(char *) panel_get_value(con6_item));
  strcpy(con7_str,(char *) panel_get_value(con7_item));
  strcpy(con8_str,(char *) panel_get_value(con8_item));
  sscanf(number_str,"%d",&number);
  sscanf(con1_str,"%lf",&con[0]);
  sscanf(con2_str,"%lf",&con[1]);
  sscanf(con3_str,"%lf",&con[2]);
  sscanf(con4_str,"%lf",&con[3]);
  sscanf(con5_str,"%lf",&con[4]);
  sscanf(con6_str,"%lf",&con[5]);
  sscanf(con7_str,"%lf",&con[6]);
  sscanf(con8_str,"%lf",&con[7]);
}
static void contours_done()
{   enter_contours();
   window_set(popup_frame, FRAME_NO_CONFIRM, TRUE, 0);
    window_destroy(popup_frame);
    if(dimensions_active==1)dimensions_active=0;
}
static void proc_contour_data(pw,h)
Pixwin *pw;
double h;
{	double xi,eta,x1,y1,v0,v1,fac;
	int i,j;
	xx[0]=0-2000.;
	yy[0]=0-2000.;
	for(j=0;j<iy1;j++)
	{for(i=0;i<ix1;i++)
	{kcon=0;
	v0=psi[j][i];
	v1=psi[j][i+1];
	if(cut(v0,v1,h)==0)
		{
		if(v1==v0){fac=0.;}else{fac=(h-v0)/(v1-v0);}
		xi=i+fac;
		eta=j;
		transform(xi,eta,&x1,&y1);
		xx[1]=x_pix(x1,xmax,width,x_off);
		yy[1]=y_pix(y1,ymin,ymax,height,y_off);
		proc_plot_data(pw);
		}
	v0=v1;
	v1=psi[j+1][i+1];
	if(cut(v0,v1,h)==0)
		{if(v1==v0){fac=0.;}else{fac=(h-v0)/(v1-v0);}
		xi=i+1;eta=j+fac;
		transform(xi,eta,&x1,&y1);
		xx[1]=x_pix(x1,xmax,width,x_off);
		yy[1]=y_pix(y1,ymin,ymax,height,y_off);		
		proc_plot_data(pw);
		}
	v0=v1;
	v1=psi[j+1][i];
	if(cut(v0,v1,h)==0)
		{if(v1==v0){fac=0.;}else{fac=(h-v0)/(v1-v0);}
		xi=i+1-fac;eta=j+1;
		transform(xi,eta,&x1,&y1);
		xx[1]=x_pix(x1,xmax,width,x_off);
		yy[1]=y_pix(y1,ymin,ymax,height,y_off);		
		proc_plot_data(pw);
		}
	v0=v1;
	v1=psi[j][i];
	if(cut(v0,v1,h)==0)
		{
		if(v1==v0){fac=0.;}else{fac=(h-v0)/(v1-v0);}
		xi=i;eta=j+1-fac;
		transform(xi,eta,&x1,&y1);
		xx[1]=x_pix(x1,xmax,width,x_off);
		yy[1]=y_pix(y1,ymin,ymax,height,y_off);		
		proc_plot_data(pw);
		}
	}}
}
draw_walls(pw)
Pixwin *pw;
{int i;double xi,eta,x1,y1;
	for(i=0;i<ix1;i++)
	{
	xx[0]=x_pix(x[0][i],xmax,width,x_off);
	yy[0]=y_pix(y[0][i],ymin,ymax,height,y_off);	
	xx[1]=x_pix(x[0][i+1],xmax,width,x_off);
	yy[1]=y_pix(y[0][i+1],ymin,ymax,height,y_off);	
	draw_line(pw,xx[0],yy[0],xx[1],yy[1]);
	}
	for(i=0;i<ix1;i++)
	{
	xx[0]=x_pix(x[iy1][i],xmax,width,x_off);
	yy[0]=y_pix(y[iy1][i],ymin,ymax,height,y_off);	
	xx[1]=x_pix(x[iy1][i+1],xmax,width,x_off);
	yy[1]=y_pix(y[iy1][i+1],ymin,ymax,height,y_off);	
	draw_line(pw,xx[0],yy[0],xx[1],yy[1]);
	}
}
transform(xi,eta,x1,y1)
double xi,eta,*x1,*y1;
{
	int i,j;
	double dxi,deta;
	i=xi;
	j=eta;
	dxi=xi-i;
	deta=eta-j;
	(*x1)=(1.-dxi)*(1.-deta)*x[j][i]
		+dxi*(1.-deta)*x[j][i+1]
		+dxi*deta*x[j+1][i+1]
		+deta*(1.-dxi)*x[j+1][i];
	(*y1)=(1.-dxi)*(1.-deta)*y[j][i]
		+dxi*(1.-deta)*y[j][i+1]
		+dxi*deta*y[j+1][i+1]
		+deta*(1.-dxi)*y[j+1][i];
}

	
	

proc_plot_data(pw)
Pixwin *pw;
{	if(kcon==1)
		{draw_line(pw,xx[0],yy[0],xx[1],yy[1]);}
	xx[0]=xx[1];
	yy[0]=yy[1];
	kcon=1-kcon;
}


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",
			&ix1,&iy1,&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);

	ix2=ix1/2;	iy2=iy1/2;	ix3=ix2/2;	iy3=iy2/2;
	ix4=ix3/2;	iy4=iy3/2;
	panel_set_value(x_points,ix1);
	panel_set_value(y_points,iy1);
	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(clear_data,1);
	strcpy(free_string,(char *)	panel_get_value(free_item));
	sscanf(free_string,"%lf",&free_value);
	strcpy(length_string,(char *)	panel_get_value(length_item));
	sscanf(length_string,"%lf",&length);
	

		for(j=0;j<=iy1;j++)
			{for(i=0;i<=ix1;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;
				}
			x[j][ix1+1]=x[j][1]+length;
			y[j][ix1+1]=y[j][1];
				p[j][ix1+1]=p[j][1];
				q[j][ix1+1]=q[j][1];

			}
	fclose(fp);
	xmin=0.;	xstep=1.;	xmax=length;
	ymin=y[0][0]-bot_depth;	ystep=0.5;	ymax=y[iy1][0]+top_depth;

	evaluate_coefficients();		
		}

}

static void read_pw()
{
	FILE *fp;
	int i,j,k,INT;
	double f,g,p1,q1,DOUBLE;
	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 
             {
		fread((char *)(&ix1),sizeof(INT),1,fp);
		fread((char *)(&iy1),sizeof(INT),1,fp);
		fread((char *)(&length),sizeof(DOUBLE),1,fp);
		fread((char *)(&top_depth),sizeof(DOUBLE),1,fp);
		fread((char *)(&bot_depth),sizeof(DOUBLE),1,fp);
		fread((char *)(&phase),sizeof(DOUBLE),1,fp);

		fread((char *)(&reynolds),sizeof(DOUBLE),1,fp);
		fread((char *)(&strouhal),sizeof(DOUBLE),1,fp);
		fread((char *)(&mean),sizeof(DOUBLE),1,fp);
		fread((char *)(&current_time),sizeof(DOUBLE),1,fp);
		fread((char *)(&dt),sizeof(DOUBLE),1,fp);

		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(reynolds_string,"%6.2f",reynolds);
		sprintf(strouhal_string,"%6.2f",strouhal);
		sprintf(mean_string,"%6.2f",strouhal);

		sprintf(iter_string,"%8.4f",current_time);
		sprintf(free_string,"%8.5f",dt);
		free_value=dt;
		panel_set_value(iter_item,iter_string);

		ix2=ix1/2;	iy2=iy1/2;
		ix3=ix2/2;	iy3=iy2/2;
		ix4=ix3/2;	iy4=iy3/2;

		panel_set_value(x_points,ix1);
		panel_set_value(y_points,iy1);
		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(free_item,free_string);
		panel_set_value(reynolds_item,reynolds_string);
		panel_set_value(strouhal_item,strouhal_string);
		panel_set_value(mean_item,mean_string);

		panel_set_value(clear_data,0);

		for(j=0;j<=iy1;j++)
			{for(i=0;i<=ix1;i++)
				{
				fread((char *)(&(x[j][i])),sizeof(DOUBLE),1,fp);
				fread((char *)(&(y[j][i])),sizeof(DOUBLE),1,fp);
				fread((char *)(&(p[j][i])),sizeof(DOUBLE),1,fp);
				fread((char *)(&(q[j][i])),sizeof(DOUBLE),1,fp);
				fread((char *)(&(psi[j][i])),sizeof(DOUBLE),1,fp);
				fread((char *)(&(old_vort[j][i])),sizeof(DOUBLE),1,fp);
				fread((char *)(&(vort[j][i])),sizeof(DOUBLE),1,fp);
				}
			x[j][ix1+1]=x[j][1]+length;
			y[j][ix1+1]=y[j][1];
			p[j][ix1+1]=p[j][1];
			q[j][ix1+1]=q[j][1];
			psi[j][ix1+1]=psi[j][1];
			vort[j][ix1+1]=vort[j][1];
			old_vort[j][ix1+1]=old_vort[j][1];
			}
		fclose(fp);
	        xmin=0.;	xstep=1.;	xmax=length;
	        ymin=y[0][0]-bot_depth; 	ystep=0.5;	ymax=y[iy1][0]+top_depth;
		evaluate_coefficients();		
		}

}
static void write_pw()
{
	FILE *fp;
	int i,j,k,INT;
	double f,g,p1,q1,DOUBLE;
	
	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 
             {	
/*
	output geometry
*/
		fwrite((char *)(&ix1),sizeof(INT),1,fp);
		fwrite((char *)(&iy1),sizeof(INT),1,fp);
		fwrite((char *)(&length),sizeof(DOUBLE),1,fp);
		fwrite((char *)(&top_depth),sizeof(DOUBLE),1,fp);
		fwrite((char *)(&bot_depth),sizeof(DOUBLE),1,fp);
		fwrite((char *)(&phase),sizeof(DOUBLE),1,fp);
/*
	output flow paramters
*/
		fwrite((char *)(&reynolds),sizeof(DOUBLE),1,fp);
		fwrite((char *)(&strouhal),sizeof(DOUBLE),1,fp);
		fwrite((char *)(&mean),sizeof(DOUBLE),1,fp);
		
		fwrite((char *)(&current_time),sizeof(DOUBLE),1,fp);
		fwrite((char *)(&dt),sizeof(DOUBLE),1,fp);
		for(j=0;j<=iy1;j++)
			{for(i=0;i<=ix1;i++)
				{
				fwrite((char *)(&(x[j][i])),sizeof(DOUBLE),1,fp);
				fwrite((char *)(&(y[j][i])),sizeof(DOUBLE),1,fp);
				fwrite((char *)(&(p[j][i])),sizeof(DOUBLE),1,fp);
				fwrite((char *)(&(q[j][i])),sizeof(DOUBLE),1,fp);
				fwrite((char *)(&(psi[j][i])),sizeof(DOUBLE),1,fp);
				fwrite((char *)(&(old_vort[j][i])),sizeof(DOUBLE),1,fp);
				fwrite((char *)(&(vort[j][i])),sizeof(DOUBLE),1,fp);
				}
			}
		fclose(fp);	
		}
}
static void evaluate_coefficients()
{
	int i,j;
	double xpi,xmi,ypi,ymi,xpj,xmj,ypj,ymj;
	for(j=1;j<iy1;j++)
	{
		for(i=1;i<=ix1;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;
		}
		alpha[j][0]=alpha[j][ix1];
		beeta[j][0]=beeta[j][ix1];
		gama[j][0]=gama[j][ix1];
		jacobian[j][0]=jacobian[j][ix1];
		alpha[j][ix1+1]=alpha[j][1];
		beeta[j][ix1+1]=beeta[j][1];
		gama[j][ix1+1]=gama[j][1];
		jacobian[j][ix1+1]=jacobian[j][1];
	}
}
step_vorticity()
{
int i,j,k;double del1,au,bu,cu,av,bv,cv,um,up,vm,vp,jac;
double fx,fy;
	set_flux();
	for(j=1;j<iy1;j++)
		{
		for(i=1;i<=ix1;i++)
		{
		jac=jacobian[j][i];
		au=courant*u[j][i+1]/(2*jac);
		cu=courant*u[j][i-1]/(2*jac);
		av=courant*v[j+1][i]/(2*jac);
		cv=courant*v[j-1][i]/(2*jac);
		del1=pe/(jac*jac);
		fx=2*del1*alpha[j][i]+2*del1*gama [j][i];
		new_vort[j][i]=(1.-omega)*vort[j][i]+
			omega*(-au*vort[j][i+1]+cu*vort[j][i-1]
			       -av*vort[j+1][i]+cv*vort[j-1][i]	
			+del1*(
			      alpha[j][i]*(vort[j][i+1]+vort[j][i-1])
			      -beeta[j][i]*(vort[j+1][i+1]-vort[j-1][i+1]
			                   -vort[j+1][i-1]+vort[j-1][i-1])/2
	   		      +gama[j][i]*(vort[j-1][i]+vort[j+1][i])
	      		      ))/fx;
		}
		new_vort[j][0]=new_vort[j][ix1];
		new_vort[j][ix1+1]=new_vort[j][1];
		}
	mgrid(&residual);
	calculate_boundary_vorticity();
	for(j=0;j<=iy1;j++)
		{
		for(i=0;i<=ix1+1;i++)
		{
		old_vort[j][i]=vort[j][i];
		vort[j][i]=new_vort[j][i];
		}
		}
	calculate_velocities();
}
calculate_boundary_vorticity()
{int i;double jac,g;
for(i=1;i<=ix1;i++)
	{
	jac=jacobian[1][i];
	g=gama[1][i];
	new_vort[0][i]=2*g*(psi[0][i]-psi[1][i])/(jac*jac);
	jac=jacobian[iy1-1][i];
	g=gama[iy1-1][i];
	new_vort[iy1][i]=2*g*(psi[iy1][i]-psi[iy1-1][i])/(jac*jac);
	}
	new_vort[0][0]=new_vort[0][ix1];
	new_vort[iy1][0]=new_vort[iy1][ix1];
	new_vort[0][ix1+1]=new_vort[0][1];
	new_vort[iy1][ix1+1]=new_vort[iy1][1];
	wmax=new_vort[0][0];
	for (i=1;i<ix1;i++)
		{if(new_vort[0][i]>wmax)wmax=new_vort[0][i];}
	

}
calculate_velocities()
{
int i,j;
	for(j=1;j<iy1;j++)
		{
		for(i=1;i<=ix1;i++)
		{
		u[j][i]=(psi[j+1][i]-psi[j-1][i])/2;
		v[j][i]=(psi[j][i-1]-psi[j][i+1])/2;
		}
		u[j][0]=u[j][ix1];
		u[j][ix1+1]=u[j][1];
		v[j][0]=v[j][ix1];
		v[j][ix1+1]=v[j][1];
		}
}
calculate_asymmetry()
{
int i,j;double f,area;

	area=(2.+top_depth+bot_depth)*(length-8.)+16.+4*(top_depth+bot_depth);
	asym=0.;
	time_residual=0.;
	for(j=1;j<=iy1/2;j++)
		{
		for(i=1;i<=ix1;i++)
			{
			f=psi[j][i]+psi[iy1-j][i];
			asym=asym+f*f*jacobian[j][i];
			
			}
		}
	asym=sqrt(asym/area);
	for(j=1;j<=iy1;j++)
		{
		for(i=1;i<=ix1;i++)
			{
			f=vort[j][i]-old_vort[j][i];
			time_residual=time_residual+absol(f);
			}
		}
	time_residual=time_residual/(ix1*iy1);
}
set_flux()
{
int i;double f;
	f=1./24.;
	for(i=0;i<=ix1+1;i++)
		{psi[iy1][i]=f;
		psi[0][i]=-f;
		}
}

static int x_pix(x1,xm,w,ox)
double x1,xm;
int w,ox;
{return(ox+(w-2*ox)*x1/xm);}
static int y_pix(y1,ym,yp,h,oy)
double y1,ym,yp;
int h,oy;
{return(oy+(h-2*oy)*(1-(y1-ym)/(yp-ym)));}
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);
}}
mgrid(difference)
double *difference;
{	int i,j,k;(*difference)=1.;k=0;
	for(j=0;j<=iy1;j++){for(i=0;i<=ix1;i++)
		{z1[j][i]=psi[j][i];r1[j][i]=0.-new_vort[j][i];}}
	while( ((*difference)>tolerance) && (k<20))
		{
		gauss1(1);	residual12();	zero2();
		gauss2(1);	residual23();	zero3();
		gauss3(1);	residual34();	zero4();
		gauss4(4);	interp43();
		gauss3(2);	interp32();
		gauss2(2);	interp21();
		gauss1(2);
		diff(difference);k++;
		}
	if(k==20)
	{printf("warning: multigrid not converged\n");}
	for(j=1;j<iy1;j++){for(i=0;i<=ix1;i++)
		{psi[j][i]=z1[j][i];}psi[j][ix1+1]=psi[j][1];}

}
interp21()
{	int i,j;for(j=0;j<=iy2;j++){for(i=1;i<ix2;i++)
	{z1[2*j][2*i]=z1[2*j][2*i]+z2[j][i];}}
	for(j=0;j<=iy1;j=j+2){for(i=1;i<ix1;i=i+2)
	{z1[j][i]=z1[j][i]+0.5*(z2[j/2][(i+1)/2]+z2[j/2][(i-1)/2]);}}
	for(j=1;j<iy1;j=j+2){for(i=2;i<=ix1-2;i=i+2)
	{z1[j][i]=z1[j][i]+0.5*(z2[(j+1)/2][i/2]+z2[(j-1)/2][i/2]);}}
	for(j=1;j<=iy1-1;j=j+2){for(i=1;i<ix1;i=i+2)
	{z1[j][i]=z1[j][i]+0.25*(z2[(j+1)/2][(i+1)/2]
	+z2[(j-1)/2][(i+1)/2]+z2[(j+1)/2][(i-1)/2]+z2[(j-1)/2][(i-1)/2]);}}
}
interp32()
{	int i,j;for(j=0;j<=iy3;j++){for(i=1;i<ix3;i++)
	{z2[2*j][2*i]=z2[2*j][2*i]+z3[j][i];}}
	for(j=0;j<=iy2;j=j+2){for(i=1;i<ix2;i=i+2)
	{z2[j][i]=z2[j][i]+0.5*(z3[j/2][(i+1)/2]+z3[j/2][(i-1)/2]);}}
	for(j=1;j<iy2;j=j+2){for(i=2;i<=ix2-2;i=i+2)
	{z2[j][i]=z2[j][i]+0.5*(z3[(j+1)/2][i/2]+z3[(j-1)/2][i/2]);}}
	for(j=1;j<=iy2-1;j=j+2){for(i=1;i<ix2;i=i+2)
	{z2[j][i]=z2[j][i]+0.25*(z3[(j+1)/2][(i+1)/2]
	+z3[(j-1)/2][(i+1)/2]+z3[(j+1)/2][(i-1)/2]+z3[(j-1)/2][(i-1)/2]);}}
}
interp43()
{	int i,j;
	for(j=0;j<=iy4;j++){for(i=1;i<ix4;i++)
	{z3[2*j][2*i]=z3[2*j][2*i]+z4[j][i];}}
	for(j=0;j<=iy3;j=j+2){for(i=1;i<ix3;i=i+2)
	{z3[j][i]=z3[j][i]+0.5*(z4[j/2][(i+1)/2]+z4[j/2][(i-1)/2]);}}
	for(j=1;j<iy3;j=j+2){for(i=2;i<=ix3-2;i=i+2)
	{z3[j][i]=z3[j][i]+0.5*(z4[(j+1)/2][i/2]+z4[(j-1)/2][i/2]);}}
	for(j=1;j<=iy3-1;j=j+2){for(i=1;i<ix3;i=i+2)
	{z3[j][i]=z3[j][i]+0.25*(z4[(j+1)/2][(i+1)/2]
	+z4[(j-1)/2][(i+1)/2]+z4[(j+1)/2][(i-1)/2]+z4[(j-1)/2][(i-1)/2]);}}
}
zero2()
{	int i,j;for(j=0;j<=iy2;j++)
	{for(i=0;i<=ix2;i++){z2[j][i]=0.;}}
}
zero3()
{	int i,j;for(j=0;j<=iy3;j++)
	{for(i=0;i<=ix3;i++){z3[j][i]=0.;}}
}
zero4()
{	int i,j;for(j=0;j<=iy4;j++)
	{for(i=0;i<=ix4;i++){z4[j][i]=0.;}}
}
residual12()
{	int i,j;for(j=2;j<iy1;j=j+2)
	{i=0;r2[j/2][i/2]=r1[j][i]-
	(alpha[j][i]*(z1[j][i+1]-2*z1[j][i]+z1[j][ix1-1])
	+gama[j][i]* (z1[j+1][i]-2*z1[j][i]+z1[j-1][i])
	-beeta[j][i]*(z1[j+1][i+1]-z1[j+1][ix1-1]
     	-z1[j-1][i+1]+z1[j-1][ix1-1])/2)/(jacobian[j][i]*jacobian[j][i])
	-p[j][i]*(z1[j][i+1]-z1[j][ix1-1])/2
	-q[j][i]*(z1[j+1][i]-z1[j-1][i])/2;
	for(i=2;i<ix1;i=i+2){r2[j/2][i/2]=r1[j][i]-
	(alpha[j][i]*(z1[j][i+1]-2*z1[j][i]+z1[j][i-1])
	+gama[j][i]* (z1[j+1][i]-2*z1[j][i]+z1[j-1][i])
	-beeta[j][i]*(z1[j+1][i+1]-z1[j+1][i-1]
     	-z1[j-1][i+1]+z1[j-1][i-1])/2)/(jacobian[j][i]*jacobian[j][i])
	-p[j][i]*(z1[j][i+1]-z1[j][i-1])/2
	-q[j][i]*(z1[j+1][i]-z1[j-1][i])/2;}
	i=ix1;r2[j/2][i/2]=r2[j/2][0];}
}
residual23()
{	int i,j;for(j=2;j<iy2;j=j+2){i=0;
	r3[j/2][i/2]=r2[j][i]-
	(alpha[2*j][2*i]*(z2[j][i+1]-2*z2[j][i]+z2[j][ix2-1])/4
	+gama[2*j][2*i]* (z2[j+1][i]-2*z2[j][i]+z2[j-1][i])/4
	-beeta[2*j][2*i]*(z2[j+1][i+1]-z2[j+1][ix2-1]
     	-z2[j-1][i+1]+z2[j-1][ix2-1])/8)/(jacobian[2*j][2*i]*jacobian[2*j][2*i])
	-p[2*j][2*i]*(z2[j][i+1]-z2[j][ix2-1])/4
	-q[2*j][2*i]*(z2[j+1][i]-z2[j-1][i])/4;
	for(i=2;i<ix2;i=i+2){r3[j/2][i/2]=r2[j][i]-
	(alpha[2*j][2*i]*(z2[j][i+1]-2*z2[j][i]+z2[j][i-1])/4
	+gama[2*j][2*i]* (z2[j+1][i]-2*z2[j][i]+z2[j-1][i])/4
	-beeta[2*j][2*i]*(z2[j+1][i+1]-z2[j+1][i-1]
     	-z2[j-1][i+1]+z2[j-1][i-1])/8)/(jacobian[2*j][2*i]*jacobian[2*j][2*i])
	-p[2*j][2*i]*(z2[j][i+1]-z2[j][i-1])/4
	-q[2*j][2*i]*(z2[j+1][i]-z2[j-1][i])/4;}
	i=ix2;r3[j/2][i/2]=r3[j/2][0];}
}
residual34()
{	int i,j;for(j=2;j<iy3;j=j+2)
	{i=0;r4[j/2][i/2]=r3[j][i]-
	(alpha[4*j][4*i]*(z3[j][i+1]-2*z3[j][i]+z3[j][ix3-1])/16
	+gama[4*j][4*i]* (z3[j+1][i]-2*z3[j][i]+z3[j-1][i])/16
	-beeta[4*j][4*i]*(z3[j+1][i+1]-z3[j+1][ix3-1]
     	-z3[j-1][i+1]+z3[j-1][ix3-1])/32)/(jacobian[4*j][4*i]*jacobian[4*j][4*i])
	-p[4*j][4*i]*(z3[j][i+1]-z3[j][ix3-1])/8
	-q[4*j][4*i]*(z3[j+1][i]-z3[j-1][i])/8;
	for(i=2;i<ix3;i=i+2){r4[j/2][i/2]=r3[j][i]-
	(alpha[4*j][4*i]*(z3[j][i+1]-2*z3[j][i]+z3[j][i-1])/16
	+gama[4*j][4*i]* (z3[j+1][i]-2*z3[j][i]+z3[j-1][i])/16
	-beeta[4*j][4*i]*(z3[j+1][i+1]-z3[j+1][i-1]
     	-z3[j-1][i+1]+z3[j-1][i-1])/32)/(jacobian[4*j][4*i]*jacobian[4*j][4*i])
	-p[4*j][4*i]*(z3[j][i+1]-z3[j][i-1])/8
	-q[4*j][4*i]*(z3[j+1][i]-z3[j-1][i])/8;}
	i=ix3;r4[j/2][i/2]=r4[j/2][0];}
}
gauss1(repeats)
int repeats;
{	int i,j,k;double j2,jp,jq,a,b,g;
	for(k=0;k<repeats;k++){for(j=1;j<iy1;j++)
	{i=0;j2=jacobian[j][i]*jacobian[j][i];jp=j2*p[j][i]/2;
	jq=j2*q[j][i]/2;a=alpha[j][i];b=beeta[j][i]/2;g=gama[j][i];
	z1[j][i]=((a+jp)*z1[j][i+1]+(a-jp)*z1[j][ix1-1]
	+(g+jq)*z1[j+1][i]+(g-jq)*z1[j-1][i]-b*(z1[j+1][i+1]-z1[j+1][ix1-1]
  	-z1[j-1][i+1]+z1[j-1][ix1-1])-j2*r1[j][i])/(2*(a+g));
	for(i=1;i<ix1;i++){j2=jacobian[j][i]*jacobian[j][i];
	jp=j2*p[j][i]/2;jq=j2*q[j][i]/2;a=alpha[j][i];b=beeta[j][i]/2;
	g=gama[j][i];
	z1[j][i]=((a+jp)*z1[j][i+1]+(a-jp)*z1[j][i-1]
	+(g+jq)*z1[j+1][i]+(g-jq)*z1[j-1][i]-b*(z1[j+1][i+1]-z1[j+1][i-1]
  	-z1[j-1][i+1]+z1[j-1][i-1])-j2*r1[j][i])/(2*(a+g));}
	z1[j][ix1]=z1[j][0];}}
}
gauss2(repeats)
int repeats;
{	int i,j,k;double j2,jp,jq,a,b,g;
	for(k=0;k<repeats;k++){for(j=1;j<iy2;j++){i=0;
	j2=jacobian[2*j][2*i]*jacobian[2*j][2*i];jp=2*j2*p[2*j][2*i]/2;
	jq=2*j2*q[2*j][2*i]/2;a=alpha[2*j][2*i];b=beeta[2*j][2*i]/2;
	g=gama[2*j][2*i];
	z2[j][i]=((a+jp)*z2[j][i+1]+(a-jp)*z2[j][ix2-1]
	+(g+jq)*z2[j+1][i]+(g-jq)*z2[j-1][i]-b*(z2[j+1][i+1]-z2[j+1][ix2-1]
  	-z2[j-1][i+1]+z2[j-1][ix2-1])-4*j2*r2[j][i])/(2*(a+g));
	for(i=1;i<ix2;i++){j2=jacobian[2*j][2*i]*jacobian[2*j][2*i];
	jp=2*j2*p[2*j][2*i]/2;jq=2*j2*q[2*j][2*i]/2;a=alpha[2*j][2*i];
	b=beeta[2*j][2*i]/2;g=gama[2*j][2*i];
	z2[j][i]=((a+jp)*z2[j][i+1]+(a-jp)*z2[j][i-1]
	+(g+jq)*z2[j+1][i]+(g-jq)*z2[j-1][i]-b*(z2[j+1][i+1]-z2[j+1][i-1]
  	-z2[j-1][i+1]+z2[j-1][i-1])-4*j2*r2[j][i])/(2*(a+g));}
	z2[j][ix2]=z2[j][0];}}
}
gauss3(repeats)
int repeats;
{	int i,j,k;double j2,jp,jq,a,b,g;
	for(k=0;k<repeats;k++){for(j=1;j<iy3;j++){i=0;
	j2=jacobian[4*j][4*i]*jacobian[4*j][4*i];jp=4*j2*p[4*j][4*i]/2;
	jq=4*j2*q[4*j][4*i]/2;a=alpha[4*j][4*i];b=beeta[4*j][4*i]/2;
	g=gama[4*j][4*i];
	z3[j][i]=((a+jp)*z3[j][i+1]+(a-jp)*z3[j][ix3-1]
	+(g+jq)*z3[j+1][i]+(g-jq)*z3[j-1][i]-b*(z3[j+1][i+1]-z3[j+1][ix3-1]
  	-z3[j-1][i+1]+z3[j-1][ix3-1])-16*j2*r3[j][i])/(2*(a+g));
	for(i=1;i<ix3;i++){j2=jacobian[4*j][4*i]*jacobian[4*j][4*i];
	jp=j2*p[4*j][4*i]/2;jq=j2*q[4*j][4*i]/2;a=alpha[4*j][4*i];
	b=beeta[4*j][4*i]/2;g=gama[4*j][4*i];
	z3[j][i]=((a+jp)*z3[j][i+1]+(a-jp)*z3[j][i-1]
	+(g+jq)*z3[j+1][i]+(g-jq)*z3[j-1][i]-b*(z3[j+1][i+1]-z3[j+1][i-1]
  	-z3[j-1][i+1]+z3[j-1][i-1])-16*j2*r3[j][i])/(2*(a+g));}
	z3[j][ix3]=z3[j][0];}}
}
gauss4(repeats)
int repeats;
{	int i,j,k;double j2,jp,jq,a,b,g;
	for(k=0;k<repeats;k++){for(j=1;j<iy4;j++){i=0;
	j2=jacobian[8*j][8*i]*jacobian[8*j][8*i];jp=8*j2*p[8*j][8*i]/2;
	jq=8*j2*q[8*j][8*i]/2;a=alpha[8*j][8*i];b=beeta[8*j][8*i]/2;
	g=gama[8*j][8*i];
	z4[j][i]=((a+jp)*z4[j][i+1]+(a-jp)*z4[j][ix4-1]
	+(g+jq)*z4[j+1][i]+(g-jq)*z4[j-1][i]-b*(z4[j+1][i+1]-z4[j+1][ix4-1]
  	-z4[j-1][i+1]+z4[j-1][ix4-1])-64*j2*r4[j][i])/(2*(a+g));
	for(i=1;i<ix4;i++){j2=jacobian[8*j][8*i]*jacobian[8*j][8*i];
	jp=j2*p[8*j][8*i]/2;jq=j2*q[8*j][8*i]/2;a=alpha[8*j][8*i];
	b=beeta[8*j][8*i]/2;g=gama[8*j][8*i];
	z4[j][i]=((a+jp)*z4[j][i+1]+(a-jp)*z4[j][i-1]
	+(g+jq)*z4[j+1][i]+(g-jq)*z4[j-1][i]-b*(z4[j+1][i+1]-z4[j+1][i-1]
  	-z4[j-1][i+1]+z4[j-1][i-1])-64*j2*r4[j][i])/(2*(a+g));}
	z4[j][ix4]=z4[j][0];}}
}
diff(w)
double *w;
{	double d;int i,j;(*w)=0.;d=0.;for(j=1;j<iy1;j++)
	{i=0;d=(alpha[j][i]*(z1[j][i+1]-2*z1[j][i]+z1[j][ix1-1])
	+gama[j][i]* (z1[j+1][i]-2*z1[j][i]+z1[j-1][i])
	-beeta[j][i]*(z1[j+1][i+1]-z1[j+1][ix1-1]
     	-z1[j-1][i+1]+z1[j-1][ix1-1])/2)/(jacobian[j][i]*jacobian[j][i])
	-p[j][i]*(z1[j][i+1]-z1[j][ix1-1])/2
	-q[j][i]*(z1[j+1][i]-z1[j-1][i])/2-r1[j][i];(*w)=(*w)+absol(d);
	for(i=1;i<ix1;i++){d=
	(alpha[j][i]*(z1[j][i+1]-2*z1[j][i]+z1[j][i-1])
	+gama[j][i]* (z1[j+1][i]-2*z1[j][i]+z1[j-1][i])
	-beeta[j][i]*(z1[j+1][i+1]-z1[j+1][i-1]
     	-z1[j-1][i+1]+z1[j-1][i-1])/2)/(jacobian[j][i]*jacobian[j][i])
	-p[j][i]*(z1[j][i+1]-z1[j][i-1])/2
	-q[j][i]*(z1[j+1][i]-z1[j-1][i])/2
	-r1[j][i];(*w)=(*w)+absol(d);}}d=absol(z1[iy1][0]-z1[0][0]);
	if(d>0.){(*w)=(*w)/(d*ix1*iy1);}else{(*w)=(*w)/(ix1*iy1);}
}
dumpvort()
{	int i,j;
	printf("----------- vort ---------------\n");
	for(i=0;i<=4;i++)
		{
		for(j=0;j<=iy1;j++)
			{printf(" %6.3lf",new_vort[j][i]);}
		printf("\n");
		}
	printf("\n");
	for(i=ix1-3;i<ix1+1;i++)
		{
		for(j=0;j<=iy1;j++)
			{printf(" %6.3lf",new_vort[j][i]);}
		printf("\n");
		}
}
dumppsi()
{	int i,j;printf("----------- psi ---------------\n");
	for(i=0;i<=4;i++){for(j=0;j<=iy1;j++){printf(" %6.3lf",psi[j][i]);
	}printf("\n");}
	printf("\n");
	for(i=ix1-3;i<ix1+1;i++)
		{
		for(j=0;j<=iy1;j++)
			{printf(" %6.3lf",psi[j][i]);}
		printf("\n");
		}
}