/*-----------------------------------------*/
/* xphased by Thomas P. Witelski 1995/1996 */
/* witelski@math.mit.edu                   */
/* This software is provided "as is" without express or implied warranty. */
/*-----------------------------------------*/
#include<stdio.h>
#include<math.h>
#include<stdlib.h>
#include<ctype.h>
#include "fcn.h"
#include "myX.h"
#define sqr(x) ((x)*(x))
#define fabs(x) ((x)>0.0 ? (x) : -(x))
#define max(a,b) ((a)>(b)? (a):(b))
#define min(a,b) ((b)>(a)? (a):(b))
#define sgn(x) ((x)>=0.0 ? 1 : -1)
#define mymod(a,b) (fmod(a,b)<0.0 ? (b)+fmod(a,b):fmod(a,b))
#define pi M_PI

/*----------------------------------*/
/* XPR and XWD external commands used */
/* either cc -DXPR=\"`which xpr`\" */
/* OR assume that xpr and xwd are in the path for system */
/*----------------------------------*/

#define fzero(a) (fabs(a)<1e-12? 1:0)


char xpr1[]="xwd -nobdrs -name \"";
char xpr2[]="\" | xpr -portrait -compact -device ps -output ";
#define BAR1 "----------------------------------------------------------------"
#define MAIN "--------------------Main Menu (KEYBOARD)------------------------"
#define LBAR "--------------------"
#define MYCLEAR "\033[H\033[2J"
char f_text[1024];
char f_map[1024];
void *f;
double ff();
double ffx();
double ffy();
char g_text[1024];
char g_map[1024];
void *g;
double gg();
double ggx();
double ggy();
double dz=1e-6;
double dz2=2e-6;

FILE *fout;

double XMAX,YMAX;
double XMIN,YMIN;
double TMIN,TMAX;
double Mx,Bx;
double My,By;
double Mt,Bt;
double xa();
double ya();
double ta();

/*------------RK STUFF------------*/
double h=0.02;
double h0=0.02;
double h00=0.02; /* 0.002 */
/*------------RK STUFF------------*/
double hin=0.01;
int method=0;
int F();
double x[2];
double work[6];
/*------------RK STUFF------------*/
int error();
myXwin *W;
myXwin *Wx;
myXwin *Wy;
myXwin *Wp;

char xrktext[50];
double xrkx[2];
double xrky[2];

double a_val=0.0;
double b_val=0.0;
double p_time=2.0*pi/1.2;

/*======GLOBAL SWITCHES======*/
int TFORCE=0;
int FILEWRITE=0;
int SIDEDRAW=1;
int PSECTION=0;
int PERIODIC=0;
/*===========================*/



void main()
{
        restart();

        puts("\t+--------------------------------------------+");
        puts("\t| XphaseD-2 : Tom's phase plane plotter      |");
        puts("\t+--------------------------------------------+");
        puts("\t| 2nd order differential equations version   |");
        puts("\t| Email comments, suggestions and bugs to:   |");
        puts("\t| Thomas P. Witelski (witelski@math.duke.edu)|");
        puts("\t| Dept of Mathematics, Duke University       |");
        puts("\t+--------------------------------------------+");
        W=myXopen("X-Y Phase Plane",0.45,0.45,0.0,0.0,1);
        XSync(W->display,True);
        startWxy();

        drawaxes();
        puts(BAR1);
        getfcns();
        menu();

}

startWxy()
{
        Wx=myXopen("X(t) profile",0.45,0.45,0.0,0.5,0);
        Wy=myXopen("Y(t) profile",0.45,0.45,0.5,0.0,0);
}
startWp()
{
        Wp=myXopen("Poincare section",0.45,0.45,0.5,0.1);
}
endWp()
{
        myXclose(Wp);
}
endWxy()
{
        myXclose(Wx);
        myXclose(Wy);
}

menu()
{
        int c;
        double a1;
char text[80];

        while(1)
        {
                myXclean(W);
                if(SIDEDRAW)
                {
                        myXclean(Wx);
                        myXclean(Wy);
                }
        puts(MYCLEAR);
        puts(BAR1);
        puts("Differential equations:");
        printf("\tdx\n\t-- = %s\n\tdt\n\n",f_map);
        printf("\tdy\n\t-- = %s\n\tdt\n\n",g_map);

sprintf(text,"[%g:%g]",XMIN,XMAX);
printf("X-range %30s\t Parameter a=%g\n",text,a_val);
sprintf(text,"[%g:%g]",YMIN,YMAX);
printf("Y-range %30s\t Parameter b=%g\n",text,b_val);
sprintf(text,"[%g:%g]",TMIN,TMAX);
printf("T-range %30s\n",text);

puts(MAIN);
puts("E)quations           \t\tR)esize domain limits");
puts("D)irection field     \t\tN)umerical method");
puts("T)rajectories        \t\tP)rint graph");
puts("S)ingular points     \t\tC)lear graph");
puts("V)alues of parameters\t\tQ)uit");
        puts(BAR1);
INPUT:
        printf("Enter letter\t: ");
        gets(text);
        c=tolower(text[0]);
        switch(c)
        {
                case 'c':
                        drawaxes();
                        break;
                case 'p':
                        printme();
                        break;
                case 'r':
                        getxyt();
                        break;
                case 'v':
                        getparam();
                        break;
                case 'd':
                        dirplot();
                        break;
                case 'e': /* equations */
                        puts(MYCLEAR);
                        getfcns();
                        break;
                case 't': /* trajectories */
                        trajectory();
                        break;
                /* singular points */
                case 's':
                        singular();
                        break;
                case 'n':  /* numerical method */
                        getnum();
                        break;
                case 'q': /* quit */
                        myXclose(W);
                        endWxy();
                        exit(1);
                        break;
                default: goto INPUT;
                        break;
        }
}
}
getyt()
{
        int button;
        double xm,ym;
        puts(MYCLEAR);
        printf("%sGet (t,y(t)) coordinates--------\n",LBAR);
        printf("Move cursor in Y(t) window\n");
        printf("Press any mouse button for main menu\n");
        do
        {
                button=myXmouse(Wy,&xm,&ym);
                if(button==MouseMove)
                        Ycoords(xm,ym);
        }
        while(button<1);
}

getparam()
{
        double a1;
        int c;
        int tf,tg;
        char text[80];
        puts(MYCLEAR);
        printf("%sSet Parameter values (KEYBOARD)---------------\n",LBAR);
        printf("a) a=%f\n",a_val);
        printf("b) b=%f\n",b_val);
        puts(BAR1);
        printf("\nEnter letter\t : ");
        gets(text);
        c=tolower(text[0]);
        if(c!='a' && c!='b')
                return;
        printf("Enter value, %c = ",c);
        gets(text);
        if(sscanf(text,"%lf",&a1)==1)
        {
                if(c=='a')
                        a_val=a1;
                else
                        b_val=a1;
                fcn_kill(f);
                f=fcn_make(&tf,f_text);
                fcn_kill(g);
                g=fcn_make(&tg,g_text);
                TFORCE=(tf||tg);
                if(method==4)
                        TFORCE=2;
                drawaxes();
        }
}
/*----------------------------------------------------------*/
getout()
{
        char text[80];
        double input;
        int c;
puts(MYCLEAR);
printf("%sSet output option(KEYBOARD)---------\n",LBAR);
puts("P)rint graph");
puts("W)rite to file");
puts("E)nd file-write");
puts("C)lear graph");
puts("S)ection: Poincare window");
puts("T)ime profile: XY windows");
if(SIDEDRAW)
        puts("Y)-T cursor");
puts(BAR1);

        printf("\nEnter letter\t: ");
        gets(text);
        c=tolower(text[0]);
        switch(c)
        {
                case 't':
                        if(SIDEDRAW)
                        {
                                SIDEDRAW=0;
                                endWxy();
                        }
                        else
                        {
                                SIDEDRAW=1;
                                startWxy();
                                drawaxes();
                        }
                        break;
                case 's':
                        if(PSECTION)
                        {
                                PSECTION=0;
                                endWp();
                        }
                        else
                        {
BADT:
                                printf("Enter Poincare time Period\t:");
                                gets(text);
                                sscanf(text,"%lf",&p_time);
                                if(p_time<0.0)
                                        goto BADT;
                                PSECTION=1;
                                startWp();
                                drawaxes();
                        }
                        break;
                case 'y':
                        if(SIDEDRAW)
                                getyt();
                        break;
                case 'c':
                        drawaxes();
                        break;
                case 'p':
                        printme();
                        break;
                case 'w':
                        writeme();
                        break;
                case 'e':
                        fclose(fout);
                        FILEWRITE=0;
                        break;
                default:
                        break;
        }
}
/*---------------------------------------------------*/
getxyt()
{
        char text[80];
        double input;
        int c;
puts(MYCLEAR);
        printf("The x range is [%g:%g]\n",XMIN,XMAX);
        printf("The y range is [%g:%g]\n",YMIN,YMAX);
        printf("The t range is [%g:%g]\n",TMIN,TMAX);
printf("%sSet domain ranges(KEYBOARD)---------\n",LBAR);
puts("Z)oom x-y");
puts("D)ouble x-y");
puts("M)anual x-y");
puts("T)imescale t");
puts(BAR1);

        printf("Enter letter\t: ");
        gets(text);
        c=tolower(text[0]);
        switch(c)
        {
                case 'z':
                        zoom();
                        break;
                case 'd':
                        xpand();
                        break;
                case 't':
        /*=================================*/
        input= -1.0;
        while(input<=0.0)
        {
        printf("\n\nEnter Timescale (Tmax)\t:");
                gets(text);
                sscanf(text,"%lf",&input);
                if(input<=0.0)
                        printf("Error: Tmax must be positive!\n");
        }
                TMAX=input;
        drawaxes();
        break;
        /*=================================*/
                case 'm':
        /*=================================*/
        printf("\n\nEnter Xmin\t:");
        gets(text);
        sscanf(text,"%lf",&input);
        XMIN=input;
        printf("Enter Xmax\t:");
        gets(text);
        sscanf(text,"%lf",&input);
        XMAX=input;
        printf("Enter Ymin\t:");
        gets(text);
        sscanf(text,"%lf",&input);
        YMIN=input;
        printf("Enter Ymax\t:");
        gets(text);
        sscanf(text,"%lf",&input);
        YMAX=input;
        drawaxes();
        break;
        /*=================================*/
                default:
                        break;
        }
}
getnum()
{
        char text[80];
        int c;
        int tf,tg;
puts(MYCLEAR);
printf("%sSet Numerical Integration Scheme(KEYBOARD)---------\n",LBAR);
printf("S)mart Runge-Kutta\n");
printf("E)uler\n");
printf("I)mproved Euler\n");
printf("R)unge-Kutta 4th order\n");
puts(BAR1);

        printf("Enter letter\t: ");
        gets(text);
        c=tolower(text[0]);
        switch(c)
        {
                case 's':
                        method=0;
                        break;
                case 'e':
                        method=1;
                        break;
                case 'i':
                        method=2;
                        break;
                case 'r':
                        method=3;
                        break;
                default:
                        method=0;
        }
BADH:
        if(method!=0 && method!=4)
        {
                printf("\n\nEnter stepsize (h)\t:");
                gets(text);
                sscanf(text,"%lf",&hin);
                if(hin<0.0)
                {
                        printf("Error: h must be positive!\n");
                        goto BADH;
                }
        }
}

/*----------------------------------------------------------*/
xpand()
{
        double xmid,xhlf;
        double ymid,yhlf;
        xmid=0.5*(XMAX+XMIN);
        ymid=0.5*(YMAX+YMIN);
        yhlf=0.5*(YMAX-YMIN);
        xhlf=0.5*(XMAX-XMIN);

        XMAX=xmid+2.0*xhlf;
        XMIN=xmid-2.0*xhlf;
        YMAX=ymid+2.0*yhlf;
        YMIN=ymid-2.0*yhlf;
        drawaxes();
}
/*----------------------------------------------------------*/
zoom()
{
        int button;
        double xm,ym;
        double xc,yc;
        double xb,yb;

myXupdate(W);
puts(MYCLEAR);
printf("%sZoom in (MOUSE)--------------------------\n",LBAR);
printf("Click LEFT button in new upper left corner of detail area\n");
                do
                {
                        button=myXmouse(W,&xm,&ym);
                        if(button==MouseMove)
                                Xcoords(xm,ym);
                }
                while(button!=MouseLeft);
        myXupdate(W);
                        xc=(XMAX-XMIN)*xm+XMIN;
                        yc=(YMAX-YMIN)*ym+YMIN;
printf("Click LEFT button in new lower right corner of detail area\n");
                do
                {
                        button=myXmouse(W,&xm,&ym);
                        if(button==MouseMove)
                                Xzoombox(xc,yc,xm,ym);
                }
                while(button!=MouseLeft);
                        xb=(XMAX-XMIN)*xm+XMIN;
                        yb=(YMAX-YMIN)*ym+YMIN;

if(fzero(xc-xb)) return;
if(fzero(yc-yb)) return;

myXprint(W,0.2,0.5,"Click to continue");
                do
                {
                        button=myXmouse(W,&xm,&ym);
                }
                while(button!=MouseLeft);


        XMIN=min(xc,xb);
        XMAX=max(xc,xb);
        YMIN=min(yc,yb);
        YMAX=max(yc,yb);
        drawaxes();
}


/*----------------------------------------------------------*/
/*----------------------------------------------------------*/
/* xwd -name \"X-Y Phase Plane" ... | xpr - device ps -output ... */
/*----------------------------------------------------------*/
/*----------------------------------------------------------*/
writeme()
{
        char file[20];
        if(FILEWRITE)
        {
                printf("Output to file already set!\n");
                return;
        }
        FILEWRITE=1;
        drawaxes();
        strcpy(file,"trjXXXXXX");
        mktemp(file);
        fout=fopen(file,"w");
        printf("------------------Save data to file-------\n");
        printf("Subsequent trajectories (t,x,y) will be saved\n");
        printf("to the file \"%s\"\n",file);
        puts(BAR1);
        printf("Press RETURN to continue\n");
        gets(file);
}

printme()
{
        char file[20];
        char xcmd[200];
        char win[50];
        char text[80];
        int c;
        puts(MYCLEAR);
        strcpy(file,"pltXXXXXX");
        myXsetwrite(W,0);
        myXsetcolor(W,0);
        myXprint(W,0.05,0.05,xrktext);
        myXsetcolor(W,1);
        myXupdate(W);
        printf("%sPrint plot (KEYBOARD)--------------------\n",LBAR);
if(SIDEDRAW||PSECTION)
{
        puts("P)hase plane");
if(SIDEDRAW)
{
        puts("X)(t) profile");
        puts("Y)(t) profile");
}
if(PSECTION)
        puts("S)ection: Poincare");
puts(BAR1);
        printf("\nEnter letter\t: ");
        gets(text);
        c=tolower(text[0]);
        switch(c)
        {
                case 'p':
                        strcpy(win,"X-Y Phase Plane");
                        break;
                case 'x':
                        if(SIDEDRAW)
                                strcpy(win,"X(t) profile");
                        else
                                return;
                        break;
                case 'y':
                        if(SIDEDRAW)
                                strcpy(win,"Y(t) profile");
                        else
                                return;
                        break;
                case 's':
                        if(PSECTION)
                                strcpy(win,"Poincare section");
                        else
                                return;
                        break;
                default:
                        return;
        }
}
else
        strcpy(win,"X-Y Phase Plane");

        printf("Saving....");
        strcpy(xcmd,xpr1);
        strcat(xcmd,win);
        strcat(xcmd,xpr2);
        mktemp(file);
        strcat(xcmd,file);
        system(xcmd);
        printf("Postscript plot of %s saved to \"%s\"\n",win,file);
        printf("Print this file with the lpr command.\n");
        puts(BAR1);
printf("Press RETURN to continue\n");
gets(file);
}
/*----------------------------------------------------------*/
singular()
{
        int button;
        double xm,ym;
        double xx,yy;
puts(MYCLEAR);
printf("%sCritical points (MOUSE)-----------------\n",LBAR);
printf("Buttons:\tLeft -- start/stop plot \tRight -- return to menu\n\n");
while(1)
{
                button=0;
                xm=ym=0.0;
                while(button<1)
                {
                        button=myXmouse(W,&xm,&ym);
                        if(button==MouseMove)
                                Xcoords(xm,ym);
                }
                if(button==MouseRight) return;
                if(button==MouseMiddle)
                        drawaxes();
                if(button==MouseLeft)
                {
                        xx=(XMAX-XMIN)*xm+XMIN;
                        yy=(YMAX-YMIN)*ym+YMIN;
                        getsing(&xx,&yy);
                }
}
}
/*----------------------------------------------------------*/
trajectory()
{
        int button;
        double xm,ym;
        double xx,yy;
puts(MYCLEAR);
        printf("%sTrajectory plot (MOUSE)-------------------\n",LBAR);

printf("Buttons:\tLeft -- start/stop plot\n\t\tMiddle -- clear\n\t\tRight -- back to menu\n");
while(1)
{
                button=0;
                xm=ym=0.0;
                while(button<1)
                {
                        button=myXmouse(W,&xm,&ym);
                        if(button==MouseMove)
                                Xcoords(xm,ym);
                }
                if(button==MouseMiddle)
                        drawaxes();
                if(button==MouseRight) return;
                if(button==MouseLeft)
                {
                        xx=(XMAX-XMIN)*xm+XMIN;
                        yy=(YMAX-YMIN)*ym+YMIN;
                        if(SIDEDRAW)
                        {
                                myXclear(Wx);
                                myXclear(Wy);
                                drawxy();
                        }
                        if(Xrk(xx,yy,F,1)) menu();
                        if(!TFORCE)
                        if(Xrk(xx,yy,F,-1)) menu();
                }
}
}
/*----------------------------------------------------------*/

drawaxes()
{
        Mx=1.0/(XMAX-XMIN);
        Bx= -Mx*XMIN;
        My=1.0/(YMAX-YMIN);
        By= -My*YMIN;
/*----------------------*/
        if(TFORCE>0)
                TMIN=0.0;
        else
                TMIN= -TMAX;
        Mt=1.0/(TMAX-TMIN);
        Bt= -Mt*TMIN;
/*----------------------*/
        strcpy(xrktext,"");
        myXclear(W);
        myXclean(W);
        myXsetcolor(W,1);
        myXsetline(W,0);
        myXsetwrite(W,1);
        myXline(W,xa(XMIN),ya(0.0),xa(XMAX),ya(0.0));
        myXline(W,xa(0.0),ya(YMIN),xa(0.0),ya(YMAX));
        myXsetline(W,1);
        myXupdate(W);
        myXreset(W);


        if(SIDEDRAW)
                drawxy();
        if(PSECTION)
                drawp();

        h0=min(h00,min(XMAX-XMIN,YMAX-YMIN)/100.0);
}
drawp()
{
        myXclear(Wp);
        myXclean(Wp);
        myXsetcolor(Wp,1);
        myXsetline(Wp,0);
        myXsetwrite(Wp,1);
        myXline(Wp,xa(XMIN),ya(0.0),xa(XMAX),ya(0.0));
        myXline(Wp,xa(0.0),ya(YMIN),xa(0.0),ya(YMAX));
        myXsetline(Wp,1);
        myXupdate(Wp);
        myXreset(Wp);
}

drawxy()
{
        myXclear(Wx);
        myXclear(Wy);
        myXclean(Wx);
        myXclean(Wy);
        myXsetline(Wx,0);
        myXsetline(Wy,0);
        myXsetwrite(Wx,1);
        myXsetwrite(Wy,1);
        myXline(Wx,xa(0.0),0.0,xa(0.0),1.0);
        myXline(Wy,0.0,ya(0.0),1.0,ya(0.0));
        myXsetline(Wx,1);
        myXsetline(Wy,1);
        myXupdate(Wx);
        myXupdate(Wy);
        myXreset(Wx);
        myXreset(Wy);

}
/*----------------------------------------------------------*/
restart()
{
        method=0;
        PERIODIC=0;
        PSECTION=0;
        a_val=0.0;
        b_val=0.0;
        h0=h00;
        hin=0.01;
        XMAX=2.5;
        XMIN= -2.5;
        YMAX=2.0;
        YMIN= -2.0;
        TMAX=20.0;  /* */
        TMIN= -TMAX;
}
/*----------------------------------------------------------*/

getfcns()
{
        int err;
        int tf,tg;
        restart();
        printf("%sSet Equations(KEYBOARD)-------------\n",LBAR);
do
{
        err=0;
        printf("Enter dx/dt=F(x,y)= ");
        gets(f_text);

        fcn_kill(f);
        strcpy(f_map,f_text);
        f=fcn_make(&tf,f_text);
        err=fcn_check(f);
        if(err) printf("Please re-enter the function\n");
}while(err);

do
{
        err=0;
        printf("Enter dy/dt=G(x,y)= ");
        gets(g_text);

        fcn_kill(g);
        strcpy(g_map,g_text);
        g=fcn_make(&tg,g_text);
        err=fcn_check(g);
        if(err) printf("Please re-enter the function\n");
}while(err);
TFORCE=(tf||tg);

drawaxes();
}


/*----------------------------------------------------------*/


double xa(x)
double x;
{
        return(Mx*x+Bx);
}
double ya(y)
double y;
{
        return(My*y+By);
}
double ta(t)
double t;
{
        return(Mt*t+Bt);
}

Xzoombox(xc,yc,xm,ym)
double xc,yc;
double xm,ym;
{
        double xb=(XMAX-XMIN)*xm+XMIN;
        double yb=(YMAX-YMIN)*ym+YMIN;

        myXclean(W);
        myXsetcolor(W,1);
        myXsetwrite(W,0);
        sprintf(xrktext,"(%g,%g)",xb,yb);
        myXprint(W,0.05,0.05,xrktext);

        myXreset(W);
        myXdrawto(W,xa(xc),ya(yc));
        myXdrawto(W,xa(xb),ya(yc));
        myXdrawto(W,xa(xb),ya(yb));
        myXdrawto(W,xa(xc),ya(yb));
        myXdrawto(W,xa(xc),ya(yc));
}
Ycoords(xm,ym)
double xm,ym;
{
        double xx=(TMAX-TMIN)*xm+TMIN;
        double yy=(YMAX-YMIN)*ym+YMIN;
        double Dx,Dy,r;
        double dx,dy,kx,ky;
        int sx,sy;
        double sxy,slope;
dx=(TMAX-TMIN)/20.0;
dy=(YMAX-YMIN)/20.0;
        myXclean(Wy);
        myXclean(W);
        myXsetcolor(W,1);
        myXsetwrite(W,0);
        myXsetcolor(Wy,1);
        myXsetwrite(Wy,0);
        sprintf(xrktext,"(%g,%g)",xx,yy);
        myXprint(Wy,0.05,0.05,xrktext);
        myXline(Wy,ta(xx)-0.05,ya(yy),ta(xx)+0.05,ya(yy));
        myXline(Wy,ta(xx),ya(yy)-0.05,ta(xx),ya(yy)+0.05);
        myXline(W,xa(XMIN),ya(yy),xa(XMAX),ya(yy));
        myXupdate(Wy);
        myXupdate(W);
}
Xcoords(xm,ym)
double xm,ym;
{
        double xx=(XMAX-XMIN)*xm+XMIN;
        double yy=(YMAX-YMIN)*ym+YMIN;
        double Dx,Dy,r;
        double dx,dy,kx,ky;
        int sx,sy;
        double sxy,slope;
dx=(XMAX-XMIN)/20.0;
dy=(YMAX-YMIN)/20.0;
sxy=dy/dx;

kx=2.5*dx;
ky=2.5*dy;
ky=kx;
        myXclean(W);
        myXsetcolor(W,1);
        myXsetwrite(W,0);
        sprintf(xrktext,"(%g,%g)",xx,yy);
        myXprint(W,0.05,0.05,xrktext);
        myXline(W,xa(xx)-0.05,ya(yy),xa(xx)+0.05,ya(yy));
        myXline(W,xa(xx),ya(yy)-0.05,xa(xx),ya(yy)+0.05);
        getks(xx,yy,&kx,&ky,&sx,&sy);
        kx*=2.5;
        ky*=2.5;
        myXline(W,xa(xx),ya(yy),xa(xx+sx*kx),ya(yy+sy*ky));
        myXupdate(W);
        if(SIDEDRAW)
        {
                myXclean(Wx);
                myXclean(Wy);
                myXsetwrite(Wx,0);
                myXsetwrite(Wy,0);
                myXline(Wx,xa(xx),ta(TMIN),xa(xx),ta(TMAX));
                myXline(Wy,ta(TMIN),ya(yy),ta(TMAX),ya(yy));
                myXupdate(Wx);
                myXupdate(Wy);
        }
}


Xrk(xx,yy,F,dir) /* phase plane trajectory in a box */
double xx,yy;
int (*F)();
int dir;
{
        double tt;
        double tp,mt;
        int i;

        double XX,YY;
        double ix,iy,it;
        int done=0;
        int button;


/*-------------------forward plot --------------*/
        if(!method)
                h=h0;
        else
                h=hin;
        i=0;
        tp=tt=0.0;


        XX=x[0]=xx;
        YY=x[1]=yy;

        myXsetcolor(W,0);
        myXsetwrite(W,0);
        myXprint(W,0.05,0.05,xrktext);
        myXsetcolor(W,1);
        sprintf(xrktext,"(%g,%g)",xx,yy);
        myXprint(W,0.05,0.05,xrktext);

        myXreset(W);
        myXpoint(W,xa(xx),ya(yy));
        if(FILEWRITE)
                fprintf(fout,"\n0.0 %g %g\n",xx,yy);
        if(SIDEDRAW)
        {
                myXreset(Wx);
                myXreset(Wy);
                myXpoint(Wx,xa(xx),1.0-ta(tt));
                myXpoint(Wy,ta(tt),ya(yy));
        }

        while(!done)
        {
                switch(method)
                {
                        case 0: /* smart RK4 */
                        case 3: /* normal RK4 */
                                rk(x,dir*tt,dir*h,2,F,work);
                        break;
                        case 1: /* euler */
                                euler(x,dir*tt,dir*h,2,F,work);
                        break;
                        case 2: /* improved euler */
                                imp_euler(x,dir*tt,dir*h,2,F,work);
                        break;
                        case 4: /* maps */
                                map(x,dir*tt,dir*h,2,F,work);
                        break;

                }
                if(PERIODIC==2)
                        x[1]=mymod(x[1],YMAX-YMIN);
                if(PERIODIC>0)
                {
                        if(x[0]>XMAX || x[0]<XMIN)
                        {
                                myXreset(W);
                                myXreset(Wx);
                                myXreset(Wy);
                        }
                        x[0]=mymod(x[0],XMAX-XMIN);
                }
                if(x[0]>XMAX || x[0]<XMIN) done=1;
                if(x[1]>YMAX || x[1]<YMIN) done=1;
                tt+=h;
                tp+=h;
                if(tt>TMAX) done=1;

                ix=xa(x[0]);
                iy=ya(x[1]);
                it=ta(dir*tt);

                myXsetwrite(W,1);

                if(method==4)
                        myXpoint(W,ix,iy);
                else
                        myXdrawto(W,ix,iy);
                myXupdate(W);
                if(SIDEDRAW)
                {
                        if(method==4)
                        {
                                myXpoint(Wx,ix,1.0-it);
                                myXpoint(Wy,it,iy);
                        }
                        else
                        {
                                myXdrawto(Wx,ix,1.0-it);
                                myXdrawto(Wy,it,iy);
                        }
                        myXupdate(Wx);
                        myXupdate(Wy);
                }
                if(FILEWRITE)
                        fprintf(fout,"%g %g %g\n",dir*tt,x[0],x[1]);
        if(PSECTION && tp>p_time)
        {
                mt=tp-p_time;
                ix=xa(XX+(x[0]-XX)*mt);
                iy=ya(YY+(x[1]-YY)*mt);
                tp-=p_time;
                myXpoint(Wp,ix,iy);
                myXupdate(Wp);
        }
        XX=x[0];
        YY=x[1];

                        button=myXmouse(W,&ix,&iy);
                        if(button==MouseMiddle)
                                drawaxes();
                        if(button==MouseLeft) return(0);
                        if(button==MouseRight) return(1);
        }
/*===========================*/
                ix=xa(x[0]);
                iy=ya(x[1]);
                it=ta(dir*tt);
        myXdrawto(W,ix,iy);
        myXupdate(W);
        if(FILEWRITE)
                fflush(fout);
        if(SIDEDRAW)
        {
                myXdrawto(Wx,ix,1.0-it);
                myXdrawto(Wy,it,iy);
                myXupdate(Wx);
                myXupdate(Wy);
        }
        return(0);
}

/*--------------------------------------------------------------------------*/
map(x,tt,h,n,F,work) /* discrete map */
double *x;
double *work; /* array of n doubles */
double tt,h;
int (*F)();
int n;
{
        int i;
        double t=tt;
        double *k0;

        k0=work;
        F(x,k0,t);

        for(i=0;i<n;++i)
                x[i]+=k0[i];
}
/*--------------------------------------------------------------------------*/
euler(x,tt,h,n,F,work) /* forward euler for a system of ODEs */
double *x;
double *work; /* array of n doubles */
double tt,h;
int (*F)();
int n;
{
        int i;
        double t=tt;
        double *k0;

        k0=work;
        F(x,k0,t);

        for(i=0;i<n;++i)
                x[i]+=h*k0[i];
}
/*--------------------------------------------------------------------------*/
imp_euler(x,tt,h,n,F,work)/* Improved Euler for a system of ODEs */
double *x;
double *work; /* array of 2*n doubles */
double tt,h;
int (*F)();
int n;
{
        int i;
        double t;

        double *k1,*k0;
        double *x1;
        k0=work;
        k1=k0+n;
        x1=k0+2*n;

        t=tt;
        F(x,k0,t);

        for(i=0;i<n;++i)
                x1[i]=x[i]+h*k0[i];
        t=tt;
        F(x1,k1,t);
        for(i=0;i<n;++i)
                x[i]+=0.5*h*(k0[i]+k1[i]);
}

/*--------------------------------------------------------------------------*/
rk(x,tt,h,n,F,work)/* fourth order general Runge-Kutta for a system of ODEs */
double *x;
double *work; /* array of 3*n doubles */
double tt,h;
int (*F)();
int n;
{
        int i;
        double t;
        double h2=h/2.0;
        double h6=h/6.0;

        double *k1,*k0;
        double *x1;
        k0=work;
        k1=k0+n;
        x1=k0+2*n;

        t=tt;
        F(x,k1,t);

        for(i=0;i<n;++i)
        {
                x1[i]=x[i]+h2*k1[i];
                k0[i]=k1[i];
        }
        t=tt+h2;
        F(x1,k1,t);

        for(i=0;i<n;++i)
        {
                x1[i]=x[i]+h2*k1[i];
                k0[i]+=2.0*k1[i];
        }
        t=tt+h2;
        F(x1,k1,t);

        for(i=0;i<n;++i)
        {
                x1[i]=x[i]+h*k1[i];
                k0[i]+=2.0*k1[i];
        }
        t=tt+h;
        F(x1,k1,t);

        for(i=0;i<n;++i)
                x[i]+=h6*(k0[i]+k1[i]);
  /* dont update time */

}
/*--------------------------------------------------------------------------*/
F(X,KK,t)
double *X,*KK,t;
{
        KK[0]=ff(X[0],X[1],t);
        KK[1]=gg(X[0],X[1],t);
}
/*----------------------------------------------------*/
double ff(x,y,t)
double x,y;
double t;
{
        return(fcn_eval(f,x,y,t));
}
double ffx(x,y)
double x,y;
{
        return((ff(x+dz,y,0.0)-ff(x-dz,y,0.0))/dz2);
}
double ggx(x,y)
double x,y;
{
        return((gg(x+dz,y,0.0)-gg(x-dz,y,0.0))/dz2);
}
double ffy(x,y)
double x,y;
{
        return((ff(x,y+dz,0.0)-ff(x,y-dz,0.0))/dz2);
}
double ggy(x,y)
double x,y;
{
        return((gg(x,y+dz,0.0)-gg(x,y-dz,0.0))/dz2);
}
double gg(x,y,t)
double x,y;
double t;
{
        return(fcn_eval(g,x,y,t));
}

/*------------------------------------------------------*/
/*------------------------------------------------------*/
/*------------------------------------------------------*/
getsing(xx,yy)
double *xx,*yy;
{
        double x= *xx;
        double y= *yy;
        double dx,dy;
        double fx,fy,f,gx,gy,g;
        double J;

        double b,root,desc;
        double norm;
        int i;
        int sign;
        double l1,l2,prod;
        double a;
        double dr=0.1*h0;

        double ex,ey;

        for(i=0;i<24;++i)
        {
                f=ff(x,y,0.0);
                g=gg(x,y,0.0);
                fx=ffx(x,y);
                fy=ffy(x,y);
                gx=ggx(x,y);
                gy=ggy(x,y);
if(fzero(f) && fzero(g)) break;
                J=fx*gy-gx*fy;
if(fzero(J))
{
        printf("Degenerate system, nonunique solution\n");
        return;
}
                dx=(g*fy-f*gy)/J;
                dy=(f*gx-fx*g)/J;
                x+=dx;
                y+=dy;

if(x>XMAX || x< XMIN) return;
if(y>YMAX || y< YMIN) return; /* singular point outside range */
        }
        *xx=x;
        *yy=y;
printf("Critical point at (%f,%f)\t",x,y);

        b=0.5*(fx+gy);
        desc=sqr(fx-gy)+4.0*gx*fy;
        root=0.5*sqrt(fabs(desc));

if(SIDEDRAW)
        drawxy();
if(fzero(desc)) /* double root */
{
        l1=b;

        fx-=l1;
        gy-=l1;

        norm=fabs(fx)+fabs(fy)+fabs(gx)+fabs(gy);
        if(fzero(norm))
        {
                if(l1>0.0)
                        printf("Unstable proper node\n");
                else
                        printf("Stable proper node\n");
                if(Xrk(x+dr,y,F,sgn(l1))) menu();
                if(Xrk(x-dr,y,F,sgn(l1))) menu();
                if(Xrk(x,y+dr,F,sgn(l1))) menu();
                if(Xrk(x,y-dr,F,sgn(l1))) menu();
        }
        else
        {
                if(l1>0.0)
                        printf("Unstable defective proper node\n");
                else
                        printf("Stable defective proper node\n");
                ex=0.0;
                ey=0.0;
                if(fabs(fx)>0.0 || fabs(gx)>0.0) ey=1.0;
                if(fabs(fy)>0.0 || fabs(gy)>0.0) ex=1.0;
                if(Xrk(x+ex*dr,y+ey*dr,F,sgn(l1))) menu();
                if(Xrk(x-ex*dr,y-ey*dr,F,sgn(l1))) menu();
        }
        return;
}
if(desc>0.0)
/*--------------------------------real roots-----------*/
{
        l1=b+root;
        l2=b-root;
        prod=l1*l2;
        if(prod>0.0)
                sign=1;
        else
        {
                if(prod<0.0)
                        sign= -1;
                else
                        sign=0;
        }
        switch(sign)
        {
                case -1:
                        printf("Saddle point\n");
                        break;
                case 0:
                        printf("Nonlinear critical point\n");
                        break;
                case 1:
                        if(l1>0.0)
                                printf("Unstable improper node\n");
                        else
                                printf("Stable improper node\n");
                        break;
        }
        evector(&ex,&ey,fx-l1,fy,gx,gy-l1);
        if(Xrk(x+ex*dr,y+ey*dr,F,sgn(l1))) menu();
        if(Xrk(x-ex*dr,y-ey*dr,F,sgn(l1))) menu();
        evector(&ex,&ey,fx-l2,fy,gx,gy-l2);
        if(Xrk(x+ex*dr,y+ey*dr,F,sgn(l2))) menu();
        if(Xrk(x-ex*dr,y-ey*dr,F,sgn(l2))) menu();
}/*-----------------complex conjugate roots----------------*/
else
{
        if(fzero(b))
        {
                printf("Center point\n");
                if(Xrk(x+5.0*dr,y,F,1)) return;
        }
        else
        {
                if(b>0.0)
                {
                        printf("Unstable spiral\n");
                        if(Xrk(x+3.0*dr,y,F,1)) menu();
                        if(Xrk(x-3.0*dr,y,F,1)) menu();
                }
                else
                {
                        printf("Stable spiral\n");
                        if(Xrk(x+3.0*dr,y,F,-1)) menu();
                        if(Xrk(x-3.0*dr,y,F,-1)) menu();
                }
        }
}
}
evector(x,y,a,b,c,d)
double *x,*y;
double a,b,c,d;
{
        if(!fzero(a*b))
        {
                *x=1.0; *y= -a/b;
                return(0);
        }
        if(!fzero(c*d))
        {
                *x=1.0; *y= -c/d;
                return(0);
        }
        *x=0.0;
        *y=0.0;
        if(!fzero(a)|| !fzero(c)) *y=1.0;
        if(!fzero(b) || !fzero(d)) *x=1.0;
        return(-1);
}



dirplot()
{
        int i,j;
        double dx,dy;
        double xx,yy;
        double kx,ky;
        int sx,sy;  /* signs */

myXsetwrite(W,1);

dx=(XMAX-XMIN)/20.0;
dy=(YMAX-YMIN)/20.0;

        for(i=0;i<21;++i)
        {
                xx=XMIN+i*dx;
                for(j=0;j<21;++j)
                {
                        yy=YMIN+j*dy;
                        getks(xx,yy,&kx,&ky,&sx,&sy);
                        kx*=0.225;
                        ky*=0.225;
                myXline(W,xa(xx-sx*kx),ya(yy-sy*ky),
                                        xa(xx+sx*kx),ya(yy+sy*ky));
                }
                myXupdate(W);
        }
        myXclean(W);
}

/*----------------------*/

getks(xx,yy,kx,ky,sx,sy)
double xx,yy;
double *kx,*ky;
int *sx,*sy;
{

double dx,dy;
double Dx,Dy;
double r,slope,sxy;

dx=(XMAX-XMIN)/20.0;
dy=(YMAX-YMIN)/20.0;
sxy=dy/dx;

*ky= *kx=dx;

        Dx=ff(xx,yy,0.0);
        Dy=gg(xx,yy,0.0);
        if(Dx>0.0) *sx=1; else *sx= -1;
        if(Dy>0.0) *sy=1; else *sy= -1;
        r=sqrt(sqr(Dx)+sqr(Dy));
        if(r>1e-10)
        {
                if(fabs(Dx)>0.0)
                {
                        slope=fabs(Dy/Dx);
                        if(slope>sxy)
                        {
                                *ky=dy;
                                *kx=fabs(dy/Dy*Dx);
                        }
                        else
                        {
                                *kx=dx;
                                *ky=fabs(dx/Dx*Dy);
                        }
                }
                else
                {
                        *kx=0.0;
                        if(fabs(Dy)>0.0)
                                *ky=dy;
                        else
                                *ky=0.0;
                }
        }
        else
        {
                *kx=0.0;
                *ky=0.0;
        }
}


/*==============*/