/*  xquantum, a simple quantum simulator for X windows.
 - this file is in the public domain (acknowledgements of course are welcome)
 - author:   Oliver Knill, knill@math.harvard.edu
 - compile:  gcc -O -o xquantum -L/usr/X11R6/lib xquantum.c -lX11
 - run:      xquantum  file.ppm 
 - version:  22 April 2000,  20 March 2001 add movie export, December 5 2025: removed again
             as it does no more work 2025
 - The x-interface was adpated from "xfires" by Michael Creutz creutz@wind.phy.bnl.gov
   Check out   http://thy.phy.bnl.gov/www/xtoys/xtoys.html for other nifty X programs
   which are (besides incredibly nice "toys") excellent templates for xwindows programs. 
 - The discrete quantum evolution used here was introduced in the paper: 
    O. Knill, "A remark on quantum dynamics", Helvetica Physica Acta, 71, 233-241, 1998
        the idea is to replace exp(i t L) by exp(i t H), where  cos(L)=H. 
        The exponential becomes then a polynomial, the evolution becomes 
        a "coupled map lattice" (= cellular automaton with continuous state space). 
 - the red and blue color channel in file.ppm represents the wave amplitude, the green 
   color channel implements the Dirichlet boundary condition. 
*/


# 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>

# define PLAYTOP 120                      /* for interface    */
# define PLAYLEFT 14                      /* dito             */

#define FILE_INDEX(x) ((x) < 10) ? "00" : (((x) < 100) ? "0" : ""), (x)

float       *image,*image1,*ptr,*ptr1;
char        stringbuffer[256];
Display     *display;
int         screen;
static char *progname;
Window      window,
            quitbutton,reset_button,
            playground,makebutton();
XColor      xcolor,colorcell;
Colormap    cmap;
GC          gc,gcpen;
int         windowwidth,windowheight;
XFontStruct *font=NULL;
int         font_height,font_width;
XSizeHints  size_hints;
int         darkcolor,lightcolor;
int         black,white;
XImage      *spinimage=NULL;
int         reset=0;

#define amp     20.0                      /* [-amp,amp] is the range 0-255    */
float *image,*image1, *ptr, *ptr1;        /* global, wave functions           */
float *current,*next;                     /* current->next->current->next     */
int x_size, y_size, c_size;               /* picture size colormap size       */
int now=1; 
int it=0;                                 /* for counting the iterations      */

void drawbutton(),
     openwindow(),
     makebuttons(),
     repaint(),
     clean_up(),
     showpic(),
     quantum(),
     get_colors();
void read_ppm_image(char *filename);
long color_trans(float r, float g, float b);


int main(argc,argv) int argc; char **argv;{
 int i,j,jj;
 int time=0;                   
 int frames=1000; 
 int interupt = 10; 
 float *ptr,*ptr1; 
 char *filename;
 XEvent report;
 progname=argv[0];

 if( (argc != 2) && (argc !=3) ){ fprintf(stderr,"Usage: ppmquantum <frames> [ppmfile] \n"); exit(1); }
 if(argc == 2) { filename=argv[1]; } else { frames=atoi(argv[1]); filename=argv[2]; } 

 read_ppm_image(filename);
 current=image;  next=image1; 
 openwindow(argc,argv);
 makebuttons();

 while(time/interupt<frames) {
   if (1==reset) { read_ppm_image(filename); now=1; reset=1-reset; }
   if (now==1) { current=image;  next=image1;}  else  {current=image1; next=image;}
   now=1-now;
   time++; 
   quantum();

   if (XPending(display)) {
     XNextEvent(display,&report);                         //  something happend 
     switch (report.type) {
     case Expose:                                         //  picture was covered
        if ((report.xexpose.window)!=window) break; 
        if (report.xexpose.count!=0) break; 
        repaint();  
        break;
     case ConfigureNotify:                                //  user resized picture
        break; 
     case ButtonPress:
        if (report.xbutton.window==quitbutton) clean_up();// user clicked Quit 
         else { if (report.xbutton.window==reset_button) { // user clicked Reset
                   reset=1-reset;
                   drawbutton(reset_button,0,0,60,18,
                   "reset",-(1-2*reset));
                 } else {                                  // clicked in picture
                   i=report.xbutton.x;                     
                   j=report.xbutton.y;
                   jj= j*x_size;
                   ptr  = current + 3 * (jj + i );
                   ptr1 = next    + 3 * (jj + i );
                   *ptr=amp/10;
                   *ptr1=amp/10;  
                 }
         }
    default:
        break;
    } 
 } 

 it++;                              
 if (it %25 ==0) {                                         // show some information
   ptr  = current + 3*(30*x_size + 50);                        
   ptr1 = next    + 3*(30*x_size + 50);
   sprintf(stringbuffer,"%d iterations",it);
   XDrawImageString(display,window,gc,8,62,
           stringbuffer,strlen(stringbuffer));
   sprintf(stringbuffer,"phi=%f,psi=%f",*ptr,*ptr1);
   XDrawImageString(display,window,gc,8,82,
           stringbuffer,strlen(stringbuffer));
 }

}
 return 0; 
}

void showpic() {
 int i,j,jj;
 float *phi,*psi,*eta;
 for (i=0;i<x_size; i++) {
   for (j=0;j<y_size; j++) {
      jj=j*x_size;
      phi   = current + 3 * (jj + i );
      psi   = phi+1;
      eta   = phi+2;
      XPutPixel(spinimage,i,j,color_trans(*phi,*psi,*eta));
   }
 }
 XPutImage(display,playground,gc,
           spinimage,0,0,0,0,x_size,y_size);
 return;
}

void repaint() {
 drawbutton(quitbutton,0,0,68,18,"quit",1);
 drawbutton(reset_button,0,0,60,18,"reset",1-2*reset);
 drawbutton(window,PLAYLEFT-4,PLAYTOP-4,x_size+8,y_size+8,NULL,4);
 sprintf(stringbuffer,"%d by %d lattice  ",x_size-2,y_size-2);
 XDrawImageString(display,window,gc,14,100,stringbuffer,strlen(stringbuffer));
 showpic();  
 return;
}
 
void openwindow(argc,argv) int argc; char **argv; {
 char *window_name="Test";
 char *icon_name="2s";
 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};
 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=WhitePixel(display,screen);
 lightcolor=white=BlackPixel(display,screen);

 if (XAllocNamedColor(display,cmap,"#0000ff",&colorcell,&xcolor))
              darkcolor=colorcell.pixel;
 if (XAllocNamedColor(display,cmap,"#000044",&colorcell,&xcolor))
              lightcolor=colorcell.pixel;

 windowwidth=(x_size+PLAYLEFT+16);
 windowheight=(PLAYTOP+y_size+10);
 window=XCreateSimpleWindow(display,RootWindow(display,screen),
   0,0,windowwidth,windowheight,4,black,lightcolor);
 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=windowheight-160;
#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

 event_mask=ExposureMask|ButtonPressMask|StructureNotifyMask;
 XSelectInput(display,window,event_mask);

 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;

 gc=XCreateGC(display,window,0,NULL);
 XSetFont(display,gc,font->fid);
 XSetForeground(display,gc,black);
 XSetBackground(display,gc,lightcolor); 

 gcpen=XCreateGC(display,window,0,NULL); 
 XSetFont(display,gcpen,font->fid);
 XSetForeground(display,gcpen,darkcolor);
 XSetBackground(display,gcpen,lightcolor); 

 XMapWindow(display,window);
 return;
}

void makebuttons() {
    int i;
    long event_mask;
    XEvent report;
    Cursor cursor;
    XWindowAttributes attr;

    XDestroySubwindows(display,window);
    quitbutton   = makebutton(4,4,68,18);
    reset_button = makebutton(70+8,4,70,18);
    playground   = XCreateSimpleWindow(display,window,
                       PLAYLEFT,PLAYTOP,x_size,y_size,0,black,white);

    event_mask = ExposureMask|ButtonReleaseMask|ButtonPressMask|
                 PointerMotionHintMask|ButtonMotionMask; 
    XSelectInput(display,playground,event_mask);
    XMapWindow(display,playground);

    /* Wait until the playground window is actually exposed */
    i = 1; 
    while (i) {
        XNextEvent(display,&report); 
        switch (report.type) {
            case Expose:
                /* STOP when the expose is for *playground*, not the other way around */
                if (report.xexpose.window == playground)
                    i = 0;
                break;
            default:
                break;
        }
    }

    if (spinimage != NULL) {
        XDestroyImage(spinimage);
        spinimage = NULL;
    }

    /* Get the real size the server ended up using */
    XGetWindowAttributes(display, playground, &attr);
    unsigned int w = attr.width;
    unsigned int h = attr.height;

    spinimage = XGetImage(display, playground,
                          0, 0, w, h, AllPlanes, ZPixmap);
    if (spinimage == NULL) {
        fprintf(stderr,"trouble creating image\n");
        exit(-1);
    }

    cursor = XCreateFontCursor(display,XC_sb_left_arrow); 
    XDefineCursor(display,playground,cursor);
    cursor = XCreateFontCursor(display,XC_hand2);
    XDefineCursor(display,quitbutton,cursor);
    XDefineCursor(display,reset_button,cursor);

    repaint();
}


Window makebutton(xoffset,yoffset,xsize,ysize)
int xoffset,yoffset,xsize,ysize;
{Window bwindow;
 long event_mask;
 bwindow=XCreateSimpleWindow(display,window,xoffset,yoffset,
        xsize,ysize,0,black,lightcolor);
 event_mask=ButtonPressMask|ExposureMask; 
 XSelectInput(display,bwindow,event_mask);
 XMapWindow(display,bwindow);
 return bwindow;
}

void drawbutton(bwindow,xoffset,yoffset,xsize,ysize,text,state)
Window bwindow;
int xoffset,yoffset,xsize,ysize,state;
char * text; {
 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);
 XSetForeground(display,gcpen,cbutton);  
 XFillRectangle(display,bwindow,gcpen,xoffset+j,yoffset+j,xsize-2*j,ysize-2*j);
 XSetForeground(display,gcpen,cleft);
 XFillRectangle(display,bwindow,gcpen,xoffset,yoffset,xsize,j);
 XFillRectangle(display,bwindow,gcpen,xoffset,yoffset,j,ysize);
 XSetForeground(display,gcpen,cright);
 for (i=0;i<j;i++) 
  {XDrawLine(display,bwindow,gcpen,
     xoffset+i,yoffset+ysize-i-1,xoffset+xsize-i-1,yoffset+ysize-i-1);
   XDrawLine(display,bwindow,gcpen,
     xoffset+xsize-i-1,yoffset+i,xoffset+xsize-i-1,yoffset+ysize-i-1);
  }
 if (NULL!=text)
   {textlength=strlen(text);
    XSetForeground(display,gcpen,ctext);
    XDrawString(display,bwindow,gcpen,
          xoffset+(xsize-textlength*font_width)/2,
          yoffset+(ysize+font_height)/2-2,
          text,textlength);
   }
 return;
}

/*---------------------------------------------------------------------------*/
/* functions clean_up                                                        */
/*---------------------------------------------------------------------------*/
void clean_up() {
 XUnloadFont(display,font->fid);
 XFreeGC(display,gc); 
 XFreeGC(display,gcpen); 
 XCloseDisplay(display);
 XDestroyImage(spinimage);
 if (NULL!=image) free((char *)  image);
 if (NULL!=image1) free((char *) image1);
 exit(1);
}

/*---------------------------------------------------------------------------*/
/* functions color_to_real and real_to_color                                 */
/*---------------------------------------------------------------------------*/
float color_to_real(int n){ return( (amp/127.0)*(n-127.0)  ); }

char real_to_color(float r){
  int l;
  l = floor(2*r*127/amp+127.00000001);
  if (l>255) {l=255;}
  if (l<0)   {l=0;  }
  return( (char) l );
}

float dirichlet(float x) { float y = (x>5.0) ? 0.0 : 1.0; return(y); }

/*---------------------------------------------------------------------------*/
/* Procedure `quantum' applies a quantum step                                */
/*---------------------------------------------------------------------------*/

void quantum(){
  int     i,j;
  int     i_right,j_right;
  int     i_left ,j_left;
  int     i_up,j_up;
  int     i_low,j_low;

  float *phi_cen, *phi_up,*phi_low,*phi_left,*phi_right;
  float *psi_cen, *psi_up,*psi_low,*psi_left,*psi_right;
  float *eta_cen, *eta_up,*eta_low,*eta_left,*eta_right;
  float *phi1, *psi1, *eta1;

  for(j=0; j<y_size; j++) {
    for (i=0; i<x_size; i++) {

      i_right   = (i+1) % x_size; j_right=j;
      i_left    = (i-1+x_size) % x_size; j_left =j;
      j_up      = (j+1) % y_size; i_up   =i;
      j_low     = (j-1+x_size) % y_size; i_low  =i;

      phi_cen   = current + 3 * (j      * x_size + i      );
      phi_up    = current + 3 * (j_up   * x_size + i_up   );
      phi_low   = current + 3 * (j_low  * x_size + i_low  );
      phi_left  = current + 3 * (j_left * x_size + i_left );
      phi_right = current + 3 * (j_right* x_size + i_right);

      psi_cen   = phi_cen+1;
      psi_up    = phi_up +1;
      psi_low   = phi_low+1;
      psi_left  = phi_left+1;
      psi_right = phi_right+1;

      eta_cen   = phi_cen+2;
      eta_up    = phi_up+2;
      eta_low   = phi_low+2;
      eta_left  = phi_left+2;
      eta_right = phi_right+2;

      phi1  = next +  3 * (j * x_size + i );
      psi1  = phi1+1;
      eta1  = phi1+2;

      *phi1 = ((*phi_right + *phi_left + *phi_up + *phi_low)/2
               - *eta_cen)*dirichlet(*psi_cen);
      *psi1 = *psi_cen;
      *eta1 = *phi_cen*dirichlet(*psi_cen);

    }
  }

  showpic(); 
}

/*---------------------------------------------------------------------------*/
/* Function `read_ppm_image'                                                 */
/*---------------------------------------------------------------------------*/
void read_ppm_image(char *filename) {
  FILE *in;
  int i;
  int l,m,n;
  char buffer[1025];

  if((in = fopen(filename,"rb")) == NULL){
    fprintf(stderr,"Can't open file %s\n", filename);
    exit(1);
  }

  if (!fgets(buffer,1025,in)) {
    fprintf(stderr,"Unsupported file format (need PPM raw)\n");
    exit(1);
  }

  do {
   if (!fgets(buffer,1025,in)) {
     fprintf(stderr, "Unexpected EOF while readingi PPM size\n");
     exit(1);
    }
  }

  while(*buffer == '#');
  x_size = atoi(strtok(buffer," "));
  y_size = atoi(strtok(NULL," "));

  if (x_size<-0 || y_size<=0) {
     fprintf(stderr,"invalid images size \n");
     exit(1);
  }

  if (!fgets(buffer,1025,in)) {
      fprintf(stderr,"Unexpected EOF \n");
      exit(1);
  }

  c_size = atoi(buffer);
  if (c_size <= 0 || c_size > 255) {
    fprintf(stderr, "Invalid max color value in PPM file\n");
    exit(1);
  }

  if((image = (float *) malloc((3*x_size*y_size+1)*sizeof(float))) == NULL){
    fprintf(stderr,"Can't load image (memory allocation error)\n");
    exit(1); 
  }
  i = 0;
  ptr = image;
  while(!feof(in) && i<3*x_size*y_size){*ptr++=color_to_real(fgetc(in));i++;}
  fclose(in);
  if(i<3*x_size*y_size){ fprintf(stderr,"Premature end of image\n"); exit(1);}

  if((image1 = (float *) malloc((3*x_size*y_size+1)*sizeof(float))) == NULL){
    fprintf(stderr," Memory allocation error \n"); exit(1); 
  }
}

/*---------------------------------------------------------------------------*/
/* Function `color_trans'                                                    */
/*---------------------------------------------------------------------------*/

long color_trans(float r, float g, float b) {
long code; 
float prob;
int red,green,blue;
red= 0x0; green = 0xf; blue= 0x0; 

// 0 black 255=0xff  blue , 256*255 0xffff green 256*256*255 red
// 0xrrggbb    

if (g==color_to_real(255)) { code=(long) ( 0xaa ); }         // boundary dark blue
else {
  prob=r*r+b*b;
  red =(int)  fabs(1*b*b);  if (red>255)   { red=255;   }
  blue=(int)  fabs(1*r*r);  if (blue>255)  { blue=255;  } 
  green=(int) fabs(0.5*prob); if (green>255) { green=255; }
  code=(long) red*256*256 + green*256 + blue; 
 }
 return(code);
}