/* Schrodinger equation laboratory  */

/* From an ideaa of Michael Creutz  */
/* creutz@wind.phy.bnl.gov          */

/* to compile:             */
/* cc -O -L/usr/X11R6/lib schrodlab.c -lX11 -lm  */

/* Version  May 3, 1998.  
/* by Claudio Destri, destri@mi.infn.it, with many changes and     */
/* additions with respect to the original by Creutz.               */
/* !!!!!! FFT CHANGED TOO !!!!!!!!!!!!!!!!!!!!!!                   */

/* this solves the Schrodinger equation in a potential             */
/* displays amplitude squared and the potential superposed         */
/* the left mouse button perturbs the wave function                */
/* the right mouse button perturbs the potential                   */
/* the damp button evolves system towards the ground state         */
/* the t -> -t button implements time reversal                     */
/* the show E button lets the avarage energy be shown as a line    */
/* the rotate V button rotates through few preset potentials       */
/* the kick button varies the average momentum                     */
/* the hbar button varies Plank's constant                         */
/* the speed button varies the speed of the evolution              */
/* the fourier button shows the amplitude squared in fourier space */
/* the project button projects out the psi's captured in the       */
/* register buttons 0 to 5 which are pressed                       */


# include <X11/Xlib.h>
# include <X11/Xutil.h>
# include <X11/Xos.h>
# include <X11/Xatom.h> 
# include <X11/cursorfont.h>
# include <stdio.h>
# include <stdlib.h>
# include <string.h>
# include <malloc.h>
# include <math.h>
# include <unistd.h>

#ifndef _COMPLEX_DEFINED
typedef struct
		{
		  double    re;
		  double    im;
		}
		  complex;
#define _COMPLEX_DEFINED
#endif

/* NSITES=number of lattice sites 
   fft (as implemented here) REQUIRES a power of two */
# define NSITES 512
# define NPOINTS 512

/* dimensions for buttons */
# define BUTTONWIDTH 78
# define BUTTONHEIGHT 18

/* some important constants */
# define SIZE 128
# define DT .05 
# define DX (SIZE/(1.*NSITES))
# define MASS 0.1
# define DAMPFACTOR 1.0

/* how many preset potentials? */
# define NPOT 7


complex xpsi[NSITES];                 /* wave function in x space             */
complex kpsi[NSITES];                 /* wave function in k space             */
complex inpsi[NSITES];                /* to monitor norm change when damped   */
complex gpsi[NSITES][6];              /* to store intersting psi's            */ 
complex pfactor[NSITES];              /* for updating potential term          */ 
complex kfactor[NSITES];              /* for updating kinetic term            */
double potential[NSITES];             /* potential in Schrodinger equation    */
double kenergy[NSITES+1];             /* array with kinetic energy            */ 
double ene;                           /* average energy                       */
double sn[NSITES];
int ik[NSITES/4];                     /* relabels k space                     */
int iamp[NPOINTS+1],ipot[NPOINTS+1];  /* X-plottable discrete counterparts    */
XPoint xpoints[2*NPOINTS];            /* for drawing waves                    */
struct timeval timeout;               /* timer for speed control              */

char getstr[1+1];                     /* to repaint capture buttons           */
char commstr[24+1];                   /* to show norm change when damped      */
short int ifgot[6]={0,0,0,0,0,0};     /* needed by the capture routines       */

short int paused=1;    /* is the system paused?                */
short int damped=0;    /* is the system damped?                */
short int trevd=1;     /* time grows (1) or decreases (-1)?    */
short int kickversor=0;/* direction of kicks                   */  
short int kicked=0;    /* signed number of kicks               */
short int hbar=5;      /* between 1 and 10                     */
short int speed=1;     /* updating delay proportional to speed; between 0 and 10 */
short int delay=0;     /* counter for implementing speed control */
short int potident=1;  /* initialize the rotation of the potential */
short int show=0;      /* do not show energy at startup */
short int fourier=0;   /* fourier=1 to show psi in momentum space  */
short int project=0;   /* when 1 project out low lying states stored in gpsi */
short int fixedbc=0;   /* fixedbc=1 to change to Dirichlet bc's    */
int ncalls=0;          /* init counter of fft calls          */
double enscale=0.3;    /* energy scale on graph (max=1)      */ 

void drawbutton(),openwindow(),makebuttons(),update(),repaint(),
      setmom(),rotpot(),cleanup(),drawit(),renorm(),energy(),
      resetbutt(),convert(),makefactor(),doproject(),fastf();

# define line(x1,y1,x2,y2) XDrawLine(display,window,gc,x1,y1,x2,y2)
# define color(i) XSetForeground(display,gc,translate[i])

/* various window stuff */
Display *display;
int screen;
static char *progname;
Window window,quitbutton,pausebutton,potbutton,bcbutton,dampbutton,
       projbutton,trevbutton,enebutton,foubutton,speedbutton,
       getbutton,hbarbutton,kickbutton,commbutton,makebutton();
XColor xcolor,colorcell;
Colormap cmap;
GC gc;
int windowwidth,windowheight,playwidth,playheight;
XFontStruct *font=NULL;
int font_height,font_width;
XSizeHints size_hints;
int darkcolor,lightcolor,black,white;
long translate[256];    /* for converting colors   */

int main(argc,argv)
int argc;
char **argv;
      {int i,j,k,p,cget;
 double psipeg,w;
 unsigned int width, height;
 XEvent report;
 progname=argv[0];
 openwindow(argc,argv);
 makebuttons();

 /* init ik */
 for (i=0;i<NSITES/8;i++)
  {ik[i]=7*NSITES/8+i;
   ik[NSITES/8+i]=i;
  }

 /* initial configuration and init of sn */
 for (i=0;i<NSITES;i++)
  {xpsi[i].re=sin(10*M_PI*i/NSITES);
   xpsi[i].im=0;
   sn[i]=sin(M_PI*i/NSITES);
   potential[i]=0.5; 
  }    
 renorm(xpsi);
 makefactor();     
 fastf(xpsi,kpsi,NSITES,fixedbc,1);

  /* loop forever, looking for events */
 while(1)
  {if ((0==XPending(display))&&(0==paused))
    {if (delay) /* don't update yet */
       {delay--;
        /* this use of select() seems a kludge to me; 
           why can't usleep() be more standard? */
        timeout.tv_sec=0;
        timeout.tv_usec=50000; /* .05 sec per delay unit */
        select(0,NULL,NULL,NULL,&timeout); 
       }
     else
      {delay=speed;
       update();
      }
    }
   else
    {XNextEvent(display,&report); 
     switch (report.type)
      {case Expose: /* the window is exposed and needs to be redrawn */
        if ((report.xexpose.window)!=window) break; /* cuts down flashing,
           but you might remove this line if things aren't being redrawn */
        if (report.xexpose.count!=0) break; /* more in queue, wait for them */
        repaint();  
        break;
       case ConfigureNotify: /* the window has been resized */
        width=report.xconfigure.width;
        height=report.xconfigure.height;
        if ((width<size_hints.min_width)||(height<size_hints.min_height))
            {fprintf(stderr,"%s: window too small to proceed.\n",progname);
             cleanup();
            } 
        else if ((width!=windowwidth)||(height!=windowheight))
            {windowwidth=width;
             windowheight=height;
             playwidth=windowwidth;
             playheight=(windowheight-3*BUTTONHEIGHT);
             makebuttons();
            } 
        break; 
       case ButtonPress:  /* which button was pressed? */
        if (report.xbutton.window==quitbutton)
	  {printf("Total number of fft calls = %d\n",ncalls);
           cleanup();
	  }
        else if (report.xbutton.window==pausebutton)
           {paused=1-paused;
            drawbutton(pausebutton,0,0,BUTTONWIDTH,BUTTONHEIGHT,
             "pause",2-4*paused);
           }
        else if (report.xbutton.window==enebutton)
	   {show=1-show;
	    energy();
	    drawit();
	    drawbutton(enebutton,0,0,BUTTONWIDTH,BUTTONHEIGHT,
			  "show E",2-4*show);
           }
        else if (report.xbutton.window==foubutton)
	   {fourier=1-fourier; 
	    drawbutton(foubutton,0,0,BUTTONWIDTH,BUTTONHEIGHT,
			  "fourier",2-4*fourier);
            drawit();
           }
        else if (report.xbutton.window==trevbutton)
           {trevd=-trevd;
	    for (i=0;i<NSITES;i++) xpsi[i].im=-xpsi[i].im;
            drawbutton(trevbutton,0,0,BUTTONWIDTH,BUTTONHEIGHT,
            "t -> -t",2*trevd);
           }
        else if (report.xbutton.window==dampbutton)
           {damped=1-damped;
            if (damped) commbutton=makebutton(0,0,2*BUTTONWIDTH,BUTTONHEIGHT);
            drawbutton(dampbutton,0,0,BUTTONWIDTH,BUTTONHEIGHT,
              "damp",2-4*damped);
            makefactor();
	    if (1-damped) XDestroyWindow(display, commbutton);
           } 
        else if (report.xbutton.window==getbutton)
           {cget=report.xbutton.x/26;
            ifgot[cget]=1-ifgot[cget];
            if (ifgot[cget]) for (i=0;i<NSITES;i++) gpsi[i][cget]=xpsi[i];
            if (1-ifgot[cget]) 
              {for (i=0;i<NSITES;i++) xpsi[i]=gpsi[i][cget];
               energy();
	       drawit();
	      }
	    sprintf(getstr,"%d",cget);
            drawbutton(getbutton,cget*BUTTONWIDTH/3,0,BUTTONWIDTH/3,
                             BUTTONHEIGHT, getstr,2-4*ifgot[cget]);
	   }
        else if (report.xbutton.window==projbutton)
           {project=1-project;
            drawbutton(projbutton,0,0,BUTTONWIDTH,BUTTONHEIGHT,
              "project",2-4*project);
            if (project) 
              {doproject(); 
               renorm(xpsi);
              }
            energy();
            drawit();
           }
        else if (report.xbutton.window==potbutton)
           {potident+=1;
	    if (potident>NPOT) potident=1;
            rotpot();
            resetbutt();
           }
        else if (report.xbutton.window==bcbutton)
           {fixedbc=1-fixedbc;
            makefactor();
            energy();
            resetbutt();
            drawit();
           } 
        else if (report.xbutton.window==speedbutton)
           {  /* reset speed */
            speed=10-(11*report.xbutton.x)/BUTTONWIDTH;
            if (speed<0) speed=0;
            if (speed>10) speed=10;
            delay=speed;
            drawbutton(speedbutton,0,0,BUTTONWIDTH,BUTTONHEIGHT,"speed",-2);
            drawbutton(speedbutton,1+((10-speed)*(BUTTONWIDTH-2))/11,1,
                      (BUTTONWIDTH-2)/11,BUTTONHEIGHT-2,"",2);
          }
        else if (report.xbutton.window==hbarbutton)
          {  /* reset hbar */
           hbar=1+(10*report.xbutton.x)/BUTTONWIDTH;
           if (hbar<1) hbar=1;
           if (hbar>10) hbar=10;
           makefactor();
	   if (show) energy();
           drawbutton(hbarbutton,0,0,BUTTONWIDTH,BUTTONHEIGHT,"hbar",-2);
           drawbutton(hbarbutton,1+((hbar-1)*(BUTTONWIDTH-2))/10,1,
                      (BUTTONWIDTH-2)/11,BUTTONHEIGHT-2,"",2);
          }
        else if (report.xbutton.window==kickbutton)
          {/* kick the wavefunction */
           i=report.xbutton.x;
           if (i<BUTTONWIDTH/4) 
             {kickversor=-1; kicked+=kickversor;}
           else if (i>3*BUTTONWIDTH/4) 
             {kickversor=1; kicked+=kickversor;}
           else kickversor=0;  
           setmom();
          }
        else
         /* button one changes wave function, two updates, and three 
            modifies potential */
         {i=report.xbutton.x;
          j=report.xbutton.y;
          if (j<playheight) /* add in Gaussian where clicked */
           {p=(i*NSITES)/playwidth; /* site where clicked */
	    psipeg=5*(playheight-j)/(1.0*playheight); 
            if (report.xbutton.button==Button1)
              {int peggy=5+1.5*psipeg; int u;
               for (k=0;k<NSITES;k++)
                {
		 w=exp(-0.00002*(p-k)*(p-k));   /* peaks at p */
		 for (u=0;u<peggy;u++) {w*=w;}  /* sharpen the peak */
                 xpsi[k].re=psipeg*w+(1-w)*xpsi[k].re;
                }
               if (project) doproject(); 
               renorm(xpsi);
	      }
            else if (report.xbutton.button==Button3)
              {int peggy=4+1.5*psipeg; int u;
               for (k=0;k<NSITES;k++)
                {   
		 w=exp(-0.00002*(p-k)*(p-k));   /* peaks at p */
		 for (u=0;u<peggy;u++) {w*=w;}  /* sharpen the peak */
                 potential[k]=(1-enscale)*psipeg*w+potential[k]*(1-w);
                }
               makefactor();
              }            
            else if (report.xbutton.button==Button2) update();
            if (show) energy();
            drawit(); 
	   }
          else update(); /* do a sweep when mouse clicked outside play area*/
         }
        break;
       default:
        break;
      } /* end of switch */
    } /* end of if XPending */
  } /* end of while(1) */
} /* end of main */


void update()
{int i,iteration;
 double temp;
 if (damped) for (i=0;i<NSITES;i++) inpsi[i]=xpsi[i];
 for (iteration=0;iteration<5;iteration++)
  {  /* half-update the potential term */
   for (i=0;i<NSITES;i++)
    {
     temp     =xpsi[i].im*pfactor[i].re+xpsi[i].re*pfactor[i].im;
     xpsi[i].re=xpsi[i].re*pfactor[i].re-xpsi[i].im*pfactor[i].im;
     xpsi[i].im=temp;
    }
     /* go to momentum space and update the kinetic term */
   fastf(xpsi,kpsi,NSITES,fixedbc,1);
   for (i=0;i<NSITES;i++)
    {
     temp     =kpsi[i].im*kfactor[i].re+kpsi[i].re*kfactor[i].im;
     kpsi[i].re=kpsi[i].re*kfactor[i].re-kpsi[i].im*kfactor[i].im;
     kpsi[i].im=temp;
    }
    /* go back to position space and do the other half */
    /* notice the fixedbc-dependent direction in the fft call */
   fastf(kpsi,xpsi,NSITES,fixedbc,2*fixedbc-1);
   for (i=0;i<NSITES;i++)
    {
     temp     =xpsi[i].im*pfactor[i].re+xpsi[i].re*pfactor[i].im;
     xpsi[i].re=xpsi[i].re*pfactor[i].re-xpsi[i].im*pfactor[i].im;
     xpsi[i].im=temp;
    }
  }
 if (project) doproject();
 if (damped+project) renorm(xpsi);
 if (damped) 
   {temp=0;
    for (i=0;i<NSITES;i++)
     {inpsi[i].re=xpsi[i].re-inpsi[i].re;
      inpsi[i].im=xpsi[i].im-inpsi[i].im;
      temp+=(inpsi[i].re*inpsi[i].re+inpsi[i].im*inpsi[i].im);    
     }
    sprintf(commstr,"%1.10g",temp);
    drawbutton(commbutton,0,0,2*BUTTONWIDTH,BUTTONHEIGHT,commstr,2); 
   }
 if (show) energy(); 
 drawit();
 return;
}

void makefactor() /* make the diagonal matrices used in updating */
{int i;
 double temp;
 for (i=0;i<NSITES;i++)
   {/* first the potential factor */
    temp=0.5*DT*(5.0/hbar)*potential[i];
    pfactor[i].re= cos(temp);
    pfactor[i].im=-sin(temp);
    if (damped)
     {temp=exp(-DAMPFACTOR*temp);
      pfactor[i].re=temp;
      pfactor[i].im=0.0;
     }
    /* then kenergy and the kinetic factor (with lattice defs!)  */
    temp=sin(M_PI*(i+fixedbc)/(1.0*NSITES*(1+fixedbc))); 
    temp=(hbar/5.0)*(2.0/DX)*temp; 
    kenergy[i]=(0.5/MASS)*temp*temp;
    temp=DT*(5.0/hbar)*kenergy[i];
    kfactor[i].re= cos(temp);
    kfactor[i].im=-sin(temp);
    if (damped)
     {temp=exp(-DAMPFACTOR*temp);
      kfactor[i].re=temp;
      kfactor[i].im=0.0;
     }
   }
 return;
}

void doproject() /* project out gpsi */
{int k,i;
 complex scalar;
 for (k=0;k<6;k++)
   {if (ifgot[k])
     {scalar.re=0.0;
      scalar.im=0.0;
      for (i=0;i<NSITES;i++)
        {scalar.re+=(gpsi[i][k].re*xpsi[i].re+gpsi[i][k].im*xpsi[i].im);
         scalar.im+=(gpsi[i][k].re*xpsi[i].im-gpsi[i][k].im*xpsi[i].re);
        }
      scalar.re*=DX;
      scalar.im*=DX;
      for (i=0;i<NSITES;i++)
        {xpsi[i].re-=(scalar.re*gpsi[i][k].re-scalar.im*gpsi[i][k].im);
         xpsi[i].im-=(scalar.re*gpsi[i][k].im+scalar.im*gpsi[i][k].re);
        }
     }
   }
 return;
}

void renorm(complex *psi) /* normalizes the wave function */
{int i;
 double norm;
  {norm=0.0;
   for(i=0;i<NSITES;i++)
     norm+=(psi[i].re*psi[i].re+psi[i].im*psi[i].im); 
   if (0.00001<norm) 
    {norm=1.0/sqrt(norm*DX);
     for (i=0;i<NSITES;i++)
      {psi[i].re*=norm;
       psi[i].im*=norm;
      }
    }  
  }
 return;
}

void setmom() /* gives psi some average momentum */
{int i ; double kickunit,kickfactor,cf,sf,temp;
  kickunit=4.0*DX;
  if (kickversor) 
   {kickfactor=kickversor/(1.0*NSITES*(1+fixedbc));
    kickfactor=kickunit*(2.0/DX)*sin(M_PI*kickfactor);
    drawbutton(kickbutton,BUTTONWIDTH/4,0,
                  BUTTONWIDTH/2,BUTTONHEIGHT,"kick",2-4*(kicked != 0));
   }
  else 
   {kickfactor=kicked/(1.0*NSITES*(1+fixedbc));
    kickfactor=-kickunit*(2.0/DX)*sin(M_PI*kickfactor);
    kicked=0;
    drawbutton(kickbutton,BUTTONWIDTH/4,0,
                  BUTTONWIDTH/2,BUTTONHEIGHT,"kick",2);
    drawbutton(kickbutton,0,0,BUTTONWIDTH/4-1,BUTTONHEIGHT,"<",2);
   }
  for (i=0;i<NSITES;i++)
   {
    cf=cos(kickfactor*i);
    sf=sin(kickfactor*i);
    temp=xpsi[i].re;
    xpsi[i].re= cf*temp - sf*xpsi[i].im;
    xpsi[i].im= sf*temp + cf*xpsi[i].im;
   }
 energy();
 return;
}

void energy() /* calculates the average energy */
{int i,j; 
 double mag;
  fastf(xpsi,kpsi,NSITES,fixedbc,1);
  ene=0.0;
  for (i=0;i<NSITES;i++)
    {mag=xpsi[i].re*xpsi[i].re+xpsi[i].im*xpsi[i].im;
     ene+=mag*potential[i];
     mag=kpsi[i].re*kpsi[i].re+kpsi[i].im*kpsi[i].im;
     ene+=mag*kenergy[i];
    }
  ene*=DX;
  if (paused) drawit();
 return;
}

void rotpot() /* defines the preset potentials */
{int i; double g=0.0005/DX;
  if (potident==1)
   {for (i=0;i<NSITES;i++) potential[i]=0.5;
   }
  if (potident==2)
   {for (i=0;i<NSITES;i++)
     {potential[i]=g*DX*DX*(i-NSITES/2)*(i-NSITES/2);
     }
   }
  if (potident==3) 
   {for (i=0;i<NSITES;i++)
     {potential[i]=g*DX*DX*(i-NSITES/2)*(i-NSITES/2);
      potential[i]*=potential[i];
     }
   }
  if (potident==4)
   {for (i=0;i<NSITES/4;i++)
     {potential[i]=potential[NSITES-i-1]=1.75;
      potential[NSITES/4+i]=potential[NSITES/2+i]=0;
     }
   }
  if (potident==5)
   {double a=0.04/DX;
    for (i=0;i<NSITES;i++)
     {potential[i]=g*DX*DX*(i-(0.5+a)*NSITES)*(i-(0.5-a)*NSITES);
      potential[i]*=potential[i];
     }
   }
  if (potident==6)
   {g=0.2/DX;
    for (i=0;i<NSITES;i++)
     {potential[i]=g*(1+cos(40*M_PI*i/NSITES));
     }
   }
  if (potident==7)
   {complex temp[NSITES];
    for (i=0;i<NSITES;i++)
     {temp[i].re=2.0*rand()/RAND_MAX;
      temp[i].im=0.0;
     }
    fft(temp,NSITES,1);
    for (i=0;i<3*NSITES/4;i++)
     {temp[NSITES/8+i].re=0.0;
      temp[NSITES/8+i].im=0.0;
     }
    fft(temp,NSITES,-1);
    for (i=0;i<NSITES;i++)
     {potential[i]=temp[i].re/NSITES;
     }
   }
  if (show) energy();
  drawit();
  makefactor();
 return;
}
 

#define min(a,b) ((a>b)? b : a)
#define max(a,b) ((a>b)? a : b)

void drawit() /* draws the wave function */
{int i,k1,k2,iene;
  if (fourier) convert(kpsi);
  else convert(xpsi);
  /* now draw it        */
  /* first the red part */
 for (i=0;i<NPOINTS;i++)
  {
   k1=iamp[i];
   k2=ipot[i];
   xpoints[i].x=xpoints[2*NPOINTS-1-i].x=(i*playwidth)/(NPOINTS-1);
   xpoints[i].y=max(k1,k2); 
   xpoints[2*NPOINTS-1-i].y=max(k2,min(k1,k2)); /* for later use */
  }
 xpoints[NPOINTS  ].x=playwidth;
 xpoints[NPOINTS  ].y=playheight; 
 xpoints[NPOINTS+1].x=0;
 xpoints[NPOINTS+1].y=playheight;
 color(0);
 XFillPolygon(display,window,gc,xpoints,NPOINTS+2,Nonconvex,CoordModeOrigin);

  /* fix the two points for the magenta part */
 for (i=NPOINTS-2;i<NPOINTS;i++)
  {
   k1=iamp[i];
   k2=ipot[i];
   xpoints[2*NPOINTS-1-i].x=(i*playwidth)/(NPOINTS-1);
   xpoints[2*NPOINTS-1-i].y=max(k2,min(k1,k2));
  }
 color(12);
 XFillPolygon(display,window,gc,xpoints,NPOINTS*2,Nonconvex,CoordModeOrigin);

  /* now the green */
 for (i=0;i<NPOINTS;i++)
  {
   k1=iamp[i];
   k2=ipot[i];
   xpoints[i].y=min(k1,k2);
  }
 color(2);
 XFillPolygon(display,window,gc,xpoints,NPOINTS*2,Nonconvex,CoordModeOrigin);

  /* and the cyan on top */
 xpoints[NPOINTS  ].x=playwidth;
 xpoints[NPOINTS  ].y=0; 
 xpoints[NPOINTS+1].x=0;
 xpoints[NPOINTS+1].y=0;
 color(3);
 XFillPolygon(display,window,gc,xpoints,NPOINTS+2,Nonconvex,CoordModeOrigin);

 /* now draw the energy line */
 if (show)
   {
    iene=playheight*(1-enscale*ene);
    color(1);
    XDrawLine(display,window,gc,0,iene,playwidth,iene);
    XDrawLine(display,window,gc,0,iene+1,playwidth,iene+1);
   }
 return;
}
  
void convert(complex *psi)
/* calculate potential and amplitude to display as integers */
{int i,ii,k1,k2,j,neff,jmax;
 double mag;
  neff=NSITES/(1+3*fourier);
  jmax=NPOINTS/neff;
  j=0;
  for (i=0;i<neff;i++)
   {ii=i;
    if (fourier && (fixedbc==0)) ii=ik[i];
    mag=sqrt(psi[ii].re*psi[ii].re+psi[ii].im*psi[ii].im)/(1+fourier);
    k1=playheight*(1-sqrt(4.0*DX)*mag);
    k2=playheight*(1-enscale*(1-fourier)*potential[i]);
    iamp[j]=max(0,min(playheight,k2));
    ipot[j]=max(0,min(playheight,k1));
    j+=jmax;
   }
  iamp[NPOINTS]=iamp[0];
  ipot[NPOINTS]=ipot[0];

  /* interpolate between NSITES to NPOINTS */
  for (j=1;j<jmax;j++)
   for (i=0;i<neff;i++)
    {
     iamp[i*jmax+j]=((jmax-j)*iamp[i*jmax]+j*iamp[(i+1)*jmax])/jmax;
     ipot[i*jmax+j]=((jmax-j)*ipot[i*jmax]+j*ipot[(i+1)*jmax])/jmax;
    }
 return;
}


void repaint()  /* this fixes the window up whenever it is uncovered */
{int j;color(9);
 drawbutton(enebutton,    0,0,BUTTONWIDTH,BUTTONHEIGHT,"show E", 2-4*show);
 drawbutton(foubutton,    0,0,BUTTONWIDTH,BUTTONHEIGHT,"fourier", 2-4*fourier);
 drawbutton(trevbutton,   0,0,BUTTONWIDTH,BUTTONHEIGHT,"t -> -t", 2*trevd);
 drawbutton(dampbutton,   0,0,BUTTONWIDTH,BUTTONHEIGHT,"damp",   2-4*damped);  
 for (j=0;j<6;j++)
   {sprintf(getstr,"%d",j);
    drawbutton(getbutton,j*BUTTONWIDTH/3,0,BUTTONWIDTH/3,
                             BUTTONHEIGHT, getstr,2-4*ifgot[j]);
   }
 drawbutton(projbutton,   0,0,BUTTONWIDTH,BUTTONHEIGHT,"project", 2-4*project);  

 drawbutton(pausebutton,  0,0,BUTTONWIDTH,BUTTONHEIGHT,"pause",  2-4*paused);
 drawbutton(quitbutton,   0,0,BUTTONWIDTH,BUTTONHEIGHT,"quit",   2);
 drawbutton(potbutton,  0,0,5*BUTTONWIDTH/4,BUTTONHEIGHT,"rotate V", 2);
 drawbutton(bcbutton,  0,0,5*BUTTONWIDTH/4,BUTTONHEIGHT,"change bc", 2);
 drawbutton(speedbutton,0,0,BUTTONWIDTH,BUTTONHEIGHT,"speed",-2);
 drawbutton(speedbutton,1+((10-speed)*(BUTTONWIDTH-2))/11,1,
                      (BUTTONWIDTH-2)/11,BUTTONHEIGHT-2,"",2);
 drawbutton(hbarbutton,0,0,BUTTONWIDTH,BUTTONHEIGHT,"hbar",-2);
 drawbutton(hbarbutton,1+((hbar-1)*(BUTTONWIDTH-2))/10,1,
                      (BUTTONWIDTH-2)/11,BUTTONHEIGHT-2,"",2);
 drawbutton(kickbutton,3*BUTTONWIDTH/4,0,BUTTONWIDTH/4,BUTTONHEIGHT,">",2);
 drawbutton(kickbutton,BUTTONWIDTH/4,0,
                                BUTTONWIDTH/2,BUTTONHEIGHT,"kick",2);
 drawbutton(kickbutton,0,0,BUTTONWIDTH/4-1,BUTTONHEIGHT,"<",2);
 drawbutton(window,4*BUTTONWIDTH+1,playheight+BUTTONHEIGHT,
   playwidth-4*BUTTONWIDTH-1,BUTTONHEIGHT,
   "left mouse button perturbs psi^2 (magenta)",0);
 drawbutton(window,4*BUTTONWIDTH+1,playheight+2*BUTTONHEIGHT,
   playwidth-4*BUTTONWIDTH-1,BUTTONHEIGHT,
   "right mouse button perturbs potential V (green)",0);
 if (damped) drawbutton(commbutton,0,0,2*BUTTONWIDTH,BUTTONHEIGHT,commstr,2);
 resetbutt();
 color(4);
 XFillRectangle(display,window,gc,0,playheight,windowwidth,BUTTONHEIGHT);
  drawit();
 return;
}

void resetbutt() /* fixes potential and boundary conditions buttons */
{
 if (potident==1)
  {drawbutton(window,0,playheight+BUTTONHEIGHT,7*BUTTONWIDTH/4,
                    BUTTONHEIGHT,"V=constant",0);
  }
 if (potident==2)
  {drawbutton(window,0,playheight+BUTTONHEIGHT,7*BUTTONWIDTH/4,
                    BUTTONHEIGHT,"V=harmonic",0);
  }
 if (potident==3) 
  {drawbutton(window,0,playheight+BUTTONHEIGHT,7*BUTTONWIDTH/4,
                     BUTTONHEIGHT,"V=quartic",0);
  }
 if (potident==4)
  {drawbutton(window,0,playheight+BUTTONHEIGHT,7*BUTTONWIDTH/4,
                      BUTTONHEIGHT,"V=square well",0);
  }
 if (potident==5)
  {drawbutton(window,0,playheight+BUTTONHEIGHT,7*BUTTONWIDTH/4,
                       BUTTONHEIGHT,"V=double well",0);
  }
 if (potident==6)
  {drawbutton(window,0,playheight+BUTTONHEIGHT,7*BUTTONWIDTH/4,
                       BUTTONHEIGHT,"V=periodic",0);
  }
 if (potident==7)
  {drawbutton(window,0,playheight+BUTTONHEIGHT,7*BUTTONWIDTH/4,
                       BUTTONHEIGHT,"V=random",0);
  }
 if (fixedbc) drawbutton(window,0,playheight+2*BUTTONHEIGHT,
     7*BUTTONWIDTH/4,BUTTONHEIGHT,"bc=dirichlet",0);
 else drawbutton(window,0,playheight+2*BUTTONHEIGHT,
     7*BUTTONWIDTH/4,BUTTONHEIGHT,"bc=periodic",0);
 return;
}

void openwindow(argc,argv)
/* a lot of this is taken from basicwin in the Xlib Programming Manual */
int argc;
char **argv;
{char *window_name="Schrodinger";
 char *icon_name="waves";
 long event_mask;
 Pixmap icon_pixmap;
 char *display_name=NULL;
 int i;
# define icon_bitmap_width 16
# define icon_bitmap_height 16
 static unsigned char icon_bitmap_bits[] = {
   0x1f, 0xf8, 0x1f, 0x88, 0x1f, 0x88, 0x1f, 0x88, 0x1f, 0x88, 0x1f, 0xf8,
   0x1f, 0xf8, 0x1f, 0xf8, 0x1f, 0xf8, 0x1f, 0xf8, 0x1f, 0xf8, 0xff, 0xff,
   0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff};

/* open up the display */
 if ((display=XOpenDisplay(display_name))==NULL)
  {fprintf(stderr,"%s: cannot connect to X server %s\n",
    progname,XDisplayName(display_name));
   exit(-1);
  }
 screen=DefaultScreen(display);
 cmap=DefaultColormap(display,screen);

 darkcolor=black=BlackPixel(display,screen);
 lightcolor=white=WhitePixel(display,screen);
 translate[0]=white;
 translate[1]=black;
 translate[2]=lightcolor;
 translate[3]=darkcolor;
 for(i=4;i<256;i++)
   translate[i]=translate[i%4]; 
  /* a bunch of colors to play with */
 if (XAllocNamedColor(display,cmap,"salmon",&colorcell,&xcolor))
              darkcolor=colorcell.pixel;
 if (XAllocNamedColor(display,cmap,"wheat",&colorcell,&xcolor))
              lightcolor=colorcell.pixel;
 if (XAllocNamedColor(display,cmap,"red",&colorcell,&xcolor))
              translate[0]=colorcell.pixel;
 if (XAllocNamedColor(display,cmap,"blue",&colorcell,&xcolor))
              translate[1]=colorcell.pixel;
 if (XAllocNamedColor(display,cmap,"green",&colorcell,&xcolor))
              translate[2]=colorcell.pixel;
 if (XAllocNamedColor(display,cmap,"cyan",&colorcell,&xcolor))
              translate[3]=colorcell.pixel;
 if (XAllocNamedColor(display,cmap,"orange",&colorcell,&xcolor))
              translate[4]=colorcell.pixel;
 if (XAllocNamedColor(display,cmap,"purple",&colorcell,&xcolor))
              translate[5]=colorcell.pixel;
 if (XAllocNamedColor(display,cmap,"yellow",&colorcell,&xcolor))
              translate[6]=colorcell.pixel;
 if (XAllocNamedColor(display,cmap,"pink",&colorcell,&xcolor))
              translate[7]=colorcell.pixel;
 if (XAllocNamedColor(display,cmap,"brown",&colorcell,&xcolor))
              translate[8]=colorcell.pixel;
 if (XAllocNamedColor(display,cmap,"grey",&colorcell,&xcolor))
              translate[9]=colorcell.pixel;
 if (XAllocNamedColor(display,cmap,"turquoise",&colorcell,&xcolor))
              translate[10]=colorcell.pixel;
 if (XAllocNamedColor(display,cmap,"gold",&colorcell,&xcolor))
              translate[11]=colorcell.pixel;
 if (XAllocNamedColor(display,cmap,"magenta",&colorcell,&xcolor))
              translate[12]=colorcell.pixel;
 if (XAllocNamedColor(display,cmap,"navy",&colorcell,&xcolor))
              translate[13]=colorcell.pixel;
 if (XAllocNamedColor(display,cmap,"tan",&colorcell,&xcolor))
              translate[14]=colorcell.pixel;
 if (XAllocNamedColor(display,cmap,"violet",&colorcell,&xcolor))
              translate[15]=colorcell.pixel;
 if (XAllocNamedColor(display,cmap,"white",&colorcell,&xcolor))
              translate[16]=colorcell.pixel;
 if (XAllocNamedColor(display,cmap,"black",&colorcell,&xcolor))
              translate[17]=colorcell.pixel;
  /* feel free to type in more colors, I got bored */
 for(i=18;i<256;i++)
   translate[i]=translate[i%18];

    /* make the main window */
 windowwidth=10*BUTTONWIDTH;
 windowheight=500;
 playwidth=windowwidth;
 playheight=(windowheight-3*BUTTONHEIGHT);
 window=XCreateSimpleWindow(display,RootWindow(display,screen),
   0,0,windowwidth,windowheight,4,translate[14],lightcolor);
    /* make the icon */
 icon_pixmap=XCreateBitmapFromData(display,window,
   icon_bitmap_bits,icon_bitmap_width,icon_bitmap_height);

 size_hints.flags=PPosition | PSize | PMinSize;
 size_hints.min_width=windowwidth;
 size_hints.min_height=300;
 /* this is probably silly, and not really tested, but here goes */
#ifdef X11R3
 size_hints.x=x;
 size_hints.y=y;
 size_hints.width=windowwidth;
 size_hints.height=windowheight;
 XSetStandardProperties(display,window,window_name,icon_name,
    icon_pixmap,argv,argc,&size_hints);
#else
 {XWMHints wm_hints;
  XClassHint class_hints;
  XTextProperty windowName, iconName;
  if (XStringListToTextProperty(&window_name,1,&windowName)==0)
   {fprintf(stderr,"%s: structure allocation for windowName failed.\n"
      ,progname);
    exit(-1);
   }
  if (XStringListToTextProperty(&icon_name,1,&iconName)==0)
   {fprintf(stderr,"%s: structure allocation for iconName failed.\n"
       ,progname);
    exit(-1);
   }
  wm_hints.initial_state=NormalState;
  wm_hints.input=True;
  wm_hints.icon_pixmap=icon_pixmap;
  wm_hints.flags=StateHint|IconPixmapHint|InputHint;
  class_hints.res_name=progname;
  class_hints.res_class="Basicwin";
  XSetWMProperties(display,window,&windowName,&iconName,
       argv,argc,&size_hints,&wm_hints,&class_hints);
 }
#endif

   /* pick the events to look for */
 event_mask=ExposureMask|ButtonPressMask|StructureNotifyMask;
 XSelectInput(display,window,event_mask);
   /* pick font: 9x15 is supposed to almost always be there */
 if ((font=XLoadQueryFont(display,"9x15"))==NULL) 
   if ((font=XLoadQueryFont(display,"fixed"))==NULL)
     {fprintf(stderr,"%s: Sorry, having font problems.\n",progname);
      exit(-1);
     }
 font_height=font->ascent+font->descent;
 font_width=font->max_bounds.width;

  /* make graphics context: */
 gc=XCreateGC(display,window,0,NULL);
 XSetFont(display,gc,font->fid);
 XSetForeground(display,gc,black);
 XSetBackground(display,gc,lightcolor); 

   /* show the window */
 XMapWindow(display,window);
 return;
}


void makebuttons()
{    /* first destroy any old buttons */
 XDestroySubwindows(display,window);
     /* now make the new buttons */
 quitbutton=makebutton(3*BUTTONWIDTH,
           windowheight-BUTTONHEIGHT,BUTTONWIDTH,BUTTONHEIGHT);
 pausebutton=makebutton(3*BUTTONWIDTH,
           windowheight-2*BUTTONHEIGHT,BUTTONWIDTH,BUTTONHEIGHT);

 enebutton=makebutton(0,windowheight-3*BUTTONHEIGHT,
                        BUTTONWIDTH,BUTTONHEIGHT);
 foubutton=makebutton(BUTTONWIDTH,windowheight-3*BUTTONHEIGHT,
                        BUTTONWIDTH,BUTTONHEIGHT);
 trevbutton=makebutton(2*BUTTONWIDTH,windowheight-3*BUTTONHEIGHT,
                        BUTTONWIDTH,BUTTONHEIGHT);
 dampbutton=makebutton(3*BUTTONWIDTH,windowheight-3*BUTTONHEIGHT,
                        BUTTONWIDTH,BUTTONHEIGHT);
 getbutton=makebutton(4*BUTTONWIDTH,windowheight-3*BUTTONHEIGHT,
                        2*BUTTONWIDTH,BUTTONHEIGHT);
 projbutton=makebutton(6*BUTTONWIDTH,windowheight-3*BUTTONHEIGHT,
                        BUTTONWIDTH,BUTTONHEIGHT);

 potbutton=makebutton(7*BUTTONWIDTH/4,windowheight-2*BUTTONHEIGHT,
                        5*BUTTONWIDTH/4,BUTTONHEIGHT);
 bcbutton=makebutton(7*BUTTONWIDTH/4,windowheight-BUTTONHEIGHT,
                        5*BUTTONWIDTH/4,BUTTONHEIGHT);

 speedbutton=makebutton(windowwidth-BUTTONWIDTH,
                  windowheight-3*BUTTONHEIGHT,BUTTONWIDTH,BUTTONHEIGHT);
 hbarbutton=makebutton(windowwidth-2*BUTTONWIDTH,
                  windowheight-3*BUTTONHEIGHT,BUTTONWIDTH,BUTTONHEIGHT);
 kickbutton=makebutton(windowwidth-3*BUTTONWIDTH,
                  windowheight-3*BUTTONHEIGHT,BUTTONWIDTH,BUTTONHEIGHT);

 return;
}

void cleanup()
{XUnloadFont(display,font->fid);
 XFreeGC(display,gc); 
 XCloseDisplay(display);
 exit(0);
} 

Window makebutton(xoffset,yoffset,xsize,ysize)
int xoffset,yoffset,xsize,ysize;
/* Puts button of specified dimensions on window.  Nothing is drawn. */
{Window buttonwindow;
 long event_mask;
 buttonwindow=XCreateSimpleWindow(display,window,xoffset,yoffset,
        xsize,ysize,0,black,lightcolor);
 event_mask=ButtonPressMask|ExposureMask; /* look for mouse-button presses */
 XSelectInput(display,buttonwindow,event_mask);
 XMapWindow(display,buttonwindow);
 return buttonwindow;
}

void drawbutton(buttonwindow,xoffset,yoffset,xsize,ysize,text,state)
Window buttonwindow;
int xoffset,yoffset,xsize,ysize,state;
char * text;
/* Draw a button in buttonwindow of specified dimensions with text
   centered.  Color of button determined by sign of "state", 
   size of border determined by magnitude. 
   I don't need no stinking tookit. */
{int textlength,i,j;
 int cdark,clight,cup,cdown;
 int cleft,cright,cbutton,ctext;
 cup=lightcolor;
 cdown=darkcolor;
 cdark=black;
 clight=white;
 if (state<0)
  {cbutton=cdown;
   ctext=clight;
   cleft=cdark;
   cright=clight;
  } 
 else
  {cbutton=cup;
   ctext=cdark;
   cleft=clight;
   cright=cdark;
  }
 j=abs(state); /* width for button borders */
 XSetForeground(display,gc,cbutton);  
 XFillRectangle(display,buttonwindow,gc,xoffset+j,yoffset+j,
        xsize-2*j,ysize-2*j);
 XSetForeground(display,gc,cleft);
 XFillRectangle(display,buttonwindow,gc,xoffset,yoffset,xsize,j);
 XFillRectangle(display,buttonwindow,gc,xoffset,yoffset,j,ysize);
 XSetForeground(display,gc,cright);
 for (i=0;i<j;i++) 
  {XDrawLine(display,buttonwindow,gc,
     xoffset+i,yoffset+ysize-i-1,xoffset+xsize-i-1,yoffset+ysize-i-1);
   XDrawLine(display,buttonwindow,gc,
     xoffset+xsize-i-1,yoffset+i,xoffset+xsize-i-1,yoffset+ysize-i-1);
  }
 if (NULL!=text)
   {textlength=strlen(text);
    XSetForeground(display,gc,ctext);
    XDrawString(display,buttonwindow,gc,
          xoffset+(xsize-textlength*font_width)/2,
          yoffset+(ysize+font_height)/2-2,
          text,textlength);
   }
 XSetForeground(display,gc,black);
 return;
}

complex tmp[NSITES];
void fastf(in,out,n,bc,sign)

complex *in,*out;
int n,sign,bc;

{int i;   
 double A=1.0/sqrt((double) n/(1.0*(1+bc)));  
 complex r1,r2;
 
 ncalls+=1;
 
 if (bc)  
  {tmp[0].re=0.0;  
   tmp[0].im=0.0;
   for (i=1;i<=n/2;i++)   
    {r1.re=sn[i]*(in[i].re + in[n-i].re);
     r1.im=sn[i]*(in[i].im + in[n-i].im);
     r2.re=  0.5*(in[i].re - in[n-i].re);
     r2.im=  0.5*(in[i].im - in[n-i].im);
     tmp[i].re=   0.5*A*(r1.re + r2.re);
     tmp[i].im=   0.5*A*(r1.im + r2.im);
     tmp[n-i].re= 0.5*A*(r1.re - r2.re);
     tmp[n-i].im= 0.5*A*(r1.im - r2.im);
    }     

   fft(tmp,n,1);
    
   out[0].re=0.0;
   out[0].im=0.0;
   out[1]=tmp[0];
   for (i=1;i<n/2;i++)
    {out[2*i].re=   tmp[i].im - tmp[n-i].im;
     out[2*i].im=  -tmp[i].re + tmp[n-i].re;
     out[2*i+1].re= tmp[i].re + tmp[n-i].re + out[2*i-1].re;
     out[2*i+1].im= tmp[i].im + tmp[n-i].im + out[2*i-1].im;
    }
  }  

 else   
  {for (i=0;i<n;i++)   
    {out[i].re=A*in[i].re;  
     out[i].im=A*in[i].im;  
    }  
   fft(out,n,sign);
  }  

 return;  
} 


/*
------------------------------------------------------------------------------
				FFT.C
Fast Fourier Transform
------------------------------------------------------------------------------
*/

/****
 *
 * name 	fft
 *
 * synopsis	fft(data, n, isign);
 *		complex 	*data;
 *		int		n;
 *		int		isign;
 *
 * description: fft evaluates the Fourier transform of the complex array data.
 *		The resulting array is returned in data.
 *		n is the size of data array and must be a power of 2
 *              smaller than the power of 2 stored in btable[16].
 *		The Rader-Brenner algorithm is used.
 *		Because this routines use some memory buffers, they can
 *		terminate and exit if an memory overflow occur.
 *		(Memory requirement: 2 * n * sizeof(complex))
 *
 * REFERENCES:	H.J.NUSSBAUMER - Fast Fourier Transform and Convolution
 *		Algorithms (Springer-Verlag 1981).
 *
 * (G.Tecchiolli October 89).
 *
 ****/

#ifndef M_PI
#define M_PI 3.141592653589793
#endif

static	double	**cosi = NULL;
static	int	nco = 0;

static	int	ndigit = 0, ndata=0;

static	complex *fft_mpool = NULL;
static	int	fft_mcnt   = 0;

static double	sip[4] = {  0.0000000000000000000,
			    0.3826834323650897000,
			    0.7071067811865475500,
			    0.9238795325112867800 };

static double	sim[4] = {  -0.0000000000000000000,
			    -0.3826834323650897000,
			    -0.7071067811865475500,
			    -0.9238795325112867800 };

static double	co[4] = {  1.0000000000000000000,
			   0.9238795325112867800,
			   0.7071067811865475500,
			   0.3826834323650897000 };

static int	btable[16] = {	0x0001, 0x0002, 0x0004, 0x0008,
				0x0010, 0x0020, 0x0040, 0x0080,
				0x0100, 0x0200, 0x0400, 0x0800,
				0x1000, 0x2000, 0x4000, 0x8000  };

static fft16(x, sign)

complex *x;
int	sign;

{

complex t[27];
complex m[18];
complex s[21];
double	*si;

	if (sign < 0)
		si = &(sip[0]);
	else
		si = &(sim[0]);

	t[ 1].re = x[ 0].re + x[ 8].re;
	t[ 2].re = x[ 4].re + x[12].re;
	t[ 3].re = x[ 2].re + x[10].re;
	t[ 4].re = x[ 2].re - x[10].re;
	t[ 5].re = x[ 6].re + x[14].re;
	t[ 6].re = x[ 6].re - x[14].re;
	t[ 7].re = x[ 1].re + x[ 9].re;
	t[ 8].re = x[ 1].re - x[ 9].re;
	t[ 9].re = x[ 3].re + x[11].re;
	t[10].re = x[ 3].re - x[11].re;
	t[11].re = x[ 5].re + x[13].re;
	t[12].re = x[ 5].re - x[13].re;
	t[13].re = x[ 7].re + x[15].re;
	t[14].re = x[ 7].re - x[15].re;
	t[15].re = t[ 1].re + t[ 2].re;
	t[16].re = t[ 3].re + t[ 5].re;
	t[17].re = t[15].re + t[16].re;
	t[18].re = t[ 7].re + t[11].re;
	t[19].re = t[ 7].re - t[11].re;
	t[20].re = t[ 9].re + t[13].re;
	t[21].re = t[ 9].re - t[13].re;
	t[22].re = t[18].re + t[20].re;
	t[23].re = t[ 8].re + t[14].re;
	t[24].re = t[ 8].re - t[14].re;
	t[25].re = t[10].re + t[12].re;
	t[26].re = t[12].re - t[10].re;

	t[ 1].im = x[ 0].im + x[ 8].im;
	t[ 2].im = x[ 4].im + x[12].im;
	t[ 3].im = x[ 2].im + x[10].im;
	t[ 4].im = x[ 2].im - x[10].im;
	t[ 5].im = x[ 6].im + x[14].im;
	t[ 6].im = x[ 6].im - x[14].im;
	t[ 7].im = x[ 1].im + x[ 9].im;
	t[ 8].im = x[ 1].im - x[ 9].im;
	t[ 9].im = x[ 3].im + x[11].im;
	t[10].im = x[ 3].im - x[11].im;
	t[11].im = x[ 5].im + x[13].im;
	t[12].im = x[ 5].im - x[13].im;
	t[13].im = x[ 7].im + x[15].im;
	t[14].im = x[ 7].im - x[15].im;
	t[15].im = t[ 1].im + t[ 2].im;
	t[16].im = t[ 3].im + t[ 5].im;
	t[17].im = t[15].im + t[16].im;
	t[18].im = t[ 7].im + t[11].im;
	t[19].im = t[ 7].im - t[11].im;
	t[20].im = t[ 9].im + t[13].im;
	t[21].im = t[ 9].im - t[13].im;
	t[22].im = t[18].im + t[20].im;
	t[23].im = t[ 8].im + t[14].im;
	t[24].im = t[ 8].im - t[14].im;
	t[25].im = t[10].im + t[12].im;
	t[26].im = t[12].im - t[10].im;

	m[ 0].re =  t[17].re + t[22].re;
	m[ 1].re =  t[17].re - t[22].re;
	m[ 2].re =  t[15].re - t[16].re;
	m[ 3].re =  t[ 1].re - t[ 2].re;
	m[ 4].re =  x[ 0].re - x[ 8].re;
	m[ 5].re = (t[19].re - t[21].re) * co[2];
	m[ 6].re = (t[ 4].re - t[ 6].re) * co[2];

	m[ 0].im =  t[17].im + t[22].im;
	m[ 1].im =  t[17].im - t[22].im;
	m[ 2].im =  t[15].im - t[16].im;
	m[ 3].im =  t[ 1].im - t[ 2].im;
	m[ 4].im =  x[ 0].im - x[ 8].im;
	m[ 5].im = (t[19].im - t[21].im) * co[2];
	m[ 6].im = (t[ 4].im - t[ 6].im) * co[2];

	m[ 7].re = (t[24].re + t[26].re) * co[3];
	m[ 8].re = (co[1] + co[3]) * t[24].re;
	m[ 9].re = (co[3] - co[1]) * t[26].re;
	s[ 7].re =  m[8].re - m[7].re;
	s[ 8].re =  m[9].re - m[7].re;

	m[ 7].im = (t[24].im + t[26].im) * co[3];
	m[ 8].im = (co[1] + co[3]) * t[24].im;
	m[ 9].im = (co[3] - co[1]) * t[26].im;
	s[ 7].im =  m[8].im - m[7].im;
	s[ 8].im =  m[9].im - m[7].im;

	if (sign < 0)
		{
		  m[10].re = t[18].im - t[20].im;
		  m[10].im = t[20].re - t[18].re;
		  m[11].re = t[ 3].im - t[ 5].im;
		  m[11].im = t[ 5].re - t[ 3].re;
		  m[12].re = x[ 4].im - x[12].im;
		  m[12].im = x[12].re - x[ 4].re;
		}
	else
		{
		  m[10].re = -t[18].im + t[20].im;
		  m[10].im = -t[20].re + t[18].re;
		  m[11].re = -t[ 3].im + t[ 5].im;
		  m[11].im = -t[ 5].re + t[ 3].re;
		  m[12].re = -x[ 4].im + x[12].im;
		  m[12].im = -x[12].re + x[ 4].re;
		}

	m[13].re = (t[19].im + t[21].im) * si[2];
	m[13].im = (-t[19].re - t[21].re) * si[2];
	m[14].re = (t[4].im + t[6].im) * si[2];
	m[14].im = (-t[4].re - t[6].re) * si[2];

	m[15].re = (t[23].im + t[25].im) * si[3];
	m[15].im = (-t[23].re - t[25].re) * si[3];
	m[16].re = -t[23].im * (si[3] - si[1]);
	m[16].im = t[23].re * (si[3] - si[1]);
	m[17].re = t[25].im * (si[3] + si[1]);
	m[17].im = -t[25].re * (si[3] + si[1]);
	s[15].re = m[15].re + m[16].re;
	s[15].im = m[15].im + m[16].im;
	s[16].re = m[15].re - m[17].re;
	s[16].im = m[15].im - m[17].im;

	s[ 1].re = m[ 3].re + m[ 5].re;
	s[ 2].re = m[ 3].re - m[ 5].re;
	s[ 3].re = m[11].re + m[13].re;
	s[ 4].re = m[13].re - m[11].re;
	s[ 5].re = m[ 4].re + m[ 6].re;
	s[ 6].re = m[ 4].re - m[ 6].re;
	s[ 9].re = s[ 5].re + s[ 7].re;
	s[10].re = s[ 5].re - s[ 7].re;
	s[11].re = s[ 6].re + s[ 8].re;
	s[12].re = s[ 6].re - s[ 8].re;
	s[13].re = m[12].re + m[14].re;
	s[14].re = m[12].re - m[14].re;
	s[17].re = s[13].re + s[15].re;
	s[18].re = s[13].re - s[15].re;
	s[19].re = s[14].re + s[16].re;
	s[20].re = s[14].re - s[16].re;

	s[ 1].im = m[ 3].im + m[ 5].im;
	s[ 2].im = m[ 3].im - m[ 5].im;
	s[ 3].im = m[11].im + m[13].im;
	s[ 4].im = m[13].im - m[11].im;
	s[ 5].im = m[ 4].im + m[ 6].im;
	s[ 6].im = m[ 4].im - m[ 6].im;
	s[ 9].im = s[ 5].im + s[ 7].im;
	s[10].im = s[ 5].im - s[ 7].im;
	s[11].im = s[ 6].im + s[ 8].im;
	s[12].im = s[ 6].im - s[ 8].im;
	s[13].im = m[12].im + m[14].im;
	s[14].im = m[12].im - m[14].im;
	s[17].im = s[13].im + s[15].im;
	s[18].im = s[13].im - s[15].im;
	s[19].im = s[14].im + s[16].im;
	s[20].im = s[14].im - s[16].im;


	x[ 0].re = m[ 0].re;
	x[ 1].re = s[ 9].re + s[17].re;
	x[ 2].re = s[ 1].re + s[ 3].re;
	x[ 3].re = s[12].re - s[20].re;
	x[ 4].re = m[ 2].re + m[10].re;
	x[ 5].re = s[11].re + s[19].re;
	x[ 6].re = s[ 2].re + s[ 4].re;
	x[ 7].re = s[10].re - s[18].re;
	x[ 8].re = m[ 1].re;
	x[ 9].re = s[10].re + s[18].re;
	x[10].re = s[ 2].re - s[ 4].re;
	x[11].re = s[11].re - s[19].re;
	x[12].re = m[ 2].re - m[10].re;
	x[13].re = s[12].re + s[20].re;
	x[14].re = s[ 1].re - s[ 3].re;
	x[15].re = s[ 9].re - s[17].re;

	x[ 0].im = m[ 0].im;
	x[ 1].im = s[ 9].im + s[17].im;
	x[ 2].im = s[ 1].im + s[ 3].im;
	x[ 3].im = s[12].im - s[20].im;
	x[ 4].im = m[ 2].im + m[10].im;
	x[ 5].im = s[11].im + s[19].im;
	x[ 6].im = s[ 2].im + s[ 4].im;
	x[ 7].im = s[10].im - s[18].im;
	x[ 8].im = m[ 1].im;
	x[ 9].im = s[10].im + s[18].im;
	x[10].im = s[ 2].im - s[ 4].im;
	x[11].im = s[11].im - s[19].im;
	x[12].im = m[ 2].im - m[10].im;
	x[13].im = s[12].im + s[20].im;
	x[14].im = s[ 1].im - s[ 3].im;
	x[15].im = s[ 9].im - s[17].im;

}

static fft_init(nd)

int	nd;

{

int	i,j;

	if (cosi != NULL)
		{
		  for (i = 0; i < nco; ++i)
			  free((char *) cosi[i]);
		  free((char *) cosi);
		}

	nco  = nd - 3;
	cosi = (double **) malloc(nco * sizeof(double *));
	if (cosi == NULL)
		{
		  fprintf(stderr,"No memory available in fft_init\n");
		  exit(1);
		}

	for (i = 0; i < nco; ++i)
		{
		  cosi[i] = (double *) malloc(btable[i + 3] * sizeof(double));
		  if (cosi[i] == NULL)
			{
			  fprintf(stderr,"No memory available in fft_init\n");
			  exit(1);
			}
		  for (j = 0; j < btable[i + 3]; ++j)
			if (j == (btable[i + 2]))
				cosi[i][j] = 0;
			else
				cosi[i][j] = 0.5 / cos(2.0 * M_PI * j / btable[i

 + 4]);
		}


	if (fft_mpool != NULL)
		free((char *) fft_mpool);

	fft_mpool = (complex *) malloc(btable[nd + 1] * sizeof(complex));
	if (fft_mpool == NULL)
		{
		  fprintf(stderr, "No memory available in fft_init\n");
		  exit(1);
		}

	fft_mcnt  = 0;

}

static complex *fft_malloc(n)

int	n;

{

complex *w = fft_mpool + fft_mcnt;

	fft_mcnt += n;

#ifdef	DEBUG_MEM
	printf("Malloc: %d\n", fft_mcnt);
#endif

	return w;

}

fft_free(w,n)

complex *w;
int	n;

{

	fft_mcnt -= n;

#ifdef	DEBUG_MEM
	printf("Free  : %d\n",fft_mcnt);
#endif

}

static RBfft(data, Nd, sign)

complex data[];
int	Nd, sign;

{

register complex	*wp1,*wp2,*wp3;
register double 	*w;
register int		m;

int	N;
int	N2;
int	N4;
complex *x;
complex *A;
complex v0,v1;

	if (Nd == 4)
		{
		  fft16(data, sign);
		  return;
		}

	N  = (1 << Nd);
	N2 = N	>> 1;
	N4 = N2 >> 1;

	x = fft_malloc(N2);
	A = fft_malloc(N2);

	wp1 = x;
	wp2 = data;
	wp3 = data + N2;
	for (m = 0; m < N2; ++m, ++wp1, ++wp2, ++wp3)
		{
		  wp1->re = wp2->re + wp3->re;
		  wp1->im = wp2->im + wp3->im;
		}

	w   = cosi[Nd - 4] + 1;
	wp1 = A + 1;
	wp2 = data + 1;
	wp3 = data + 1 + N2;
	for (m = 1; m < N4; ++m, ++wp1, ++wp2, ++wp3, ++w)
		{
		  wp1->re = (wp2->re - wp3->re) * (*w);
		  wp1->im = (wp2->im - wp3->im) * (*w);
		}
	++m, ++wp1, ++wp2, ++wp3, ++w;
	for (; m < N2; ++m, ++wp1, ++wp2, ++wp3, ++w)
		{
		  wp1->re = (wp2->re - wp3->re) * (*w);
		  wp1->im = (wp2->im - wp3->im) * (*w);
		}

	A[0].re	= 0.0;
	A[0].im	= 0.0;
	A[N4].re = 0.0;
	A[N4].im = 0.0;

	RBfft(A, Nd-1, sign);
	RBfft(x, Nd-1, sign);

	if(sign < 0)
		{
		  v0.re = data[0].re - data[N2].re + data[N4].im - data[N - N4].im;
		  v0.im = data[0].im - data[N2].im - data[N4].re + data[N - N4].re;
		  v1.re = data[0].re - data[N2].re - data[N4].im + data[N - N4].im;
		  v1.im = data[0].im - data[N2].im + data[N4].re - data[N - N4].re;
		}
	else
		{
		  v0.re = data[0].re - data[N2].re - data[N4].im + data[N - N4].im;
		  v0.im = data[0].im - data[N2].im + data[N4].re - data[N - N4].re;
		  v1.re = data[0].re - data[N2].re + data[N4].im - data[N - N4].im;
		  v1.im = data[0].im - data[N2].im - data[N4].re + data[N - N4].re;
		}

	wp1 = data;
	wp2 = x;
	for (m = 0; m < N2; ++m, wp1 += 2, ++wp2)
		*wp1 = *wp2;

	wp1 = data + 1;
	wp2 = A;
	for (m = 1; m < N4; ++m, wp1 += 2, ++wp2)
		{
		  wp1->re = wp2->re + (wp2 + 1)->re + v0.re;
		  wp1->im = wp2->im + (wp2 + 1)->im + v0.im;

		  wp1 += 2;
		  ++wp2;

		  wp1->re = wp2->re + (wp2 + 1)->re + v1.re;
		  wp1->im = wp2->im + (wp2 + 1)->im + v1.im;
		}

	wp1->re = wp2->re + (wp2 + 1)->re + v0.re;
	wp1->im = wp2->im + (wp2 + 1)->im + v0.im;

	data[N - 1].re = A[N2 - 1].re + A[0].re + v1.re;
	data[N - 1].im = A[N2 - 1].im + A[0].im + v1.im;

	fft_free(A, N2);
	fft_free(x, N2);

}

fft(data, n, sign)

complex *data;
int	n, sign;

{

int	i;
int	j = 0;

if (n != ndata)
	{
	  ndata = n;

	  i = 1;
	  while (n > i)
		{
		  i <<= 1;
		  ++j;
		}

	  fft_init(j);
	  ndigit = j;
	}

RBfft(data, ndigit, sign);

}