/*
NEW: 4mar99: change >0 to >= 0 in following line
       if(geo.iside[i]>=0){ // indicates side is just series of C atoms

Need to fine tune perspective -
   Problem is particulaly noticable in shad
 Fixed by increasing cmeyedist from 60 to 120-
  However, maybe the eyedist should vary with the size of the protein that is being drawn
*/

#include "drawwindow.h"  // the data that needs to be seen by other routines
#include "screen.h"
#include "cube.h"
#include "contour.h" 
//#include "core.h"
#include "bone.h"
#include "peptide.h"
#include "math.h"
#include "sphere.h"
#include "trace.h"
#include "pdbtolev.h"
#include "model.h"
#include "menu.h"
#include "assignseq.h"
#include "fit.h"
#include "extendfit.h"
#include "editbones.h"

extern struct extbones ext; // All the data for edited bones used to extend one end of a given fit
extern struct screenvalue  sc; // Need global variable for X-menu calls
extern struct maiddatafile maidfile;// Need global variable for X-menu calls
extern struct griddata map[NMAPS]; // all the map (denfile),contour data associated with map
extern struct coredata core[NMAPS];
extern struct bonedata bone[NMAPS];
extern struct geometry geo;
extern struct spheredata sp[NSPHERES];
extern struct rotomer rot[30];
extern Widget toplevel;
extern Widget spheremenu;
extern Widget  fitmenu; // Global since needed in fit_cb
extern struct pepfit fit[NFIT]; // defined in peptide.h, see trace.c++ for definitin
extern int fitnum[NFIT];// see trace.c++ for definition
extern int side[AA][NBRN][NPTS];  // defined globally (originally in pdbtolev.c++)  because needed in assignseg.c++ 
extern float rotpos[NROT][NATOMS][3];// position of atoms in rotomer
extern int rottype;  // Type of rotomer to draw
extern int splitscreen; // if = 1 then use standard splitscreen, if = 0 use wholescreen
extern  struct st ss;  // Define before in menu.c++
extern int verbose; // if = 1,print some stuff;  = 2 print a lot of stuff

void update_view(int mstate, int ox, int nx, int oy, int ny); 
void draw_scene(Widget drawwindow,float eyepos);
//void drawpdb(struct pdbfile& pdb1);
void drawresidue(struct pdbfile& pdb1,int is,int i,
       int iprev);
void allocatelev(struct pdbfile& pdb1); 
void pick(int mx,int my);
void processhits(GLint hits,GLuint buffer[]);
void findpdbwind();
void findbonewind();
void draweditbone(struct extbones ext);

static float objmat[16] = {   // global for all files in drawwindow.c++
         1.0, 0.0, 0.0, 0.0,
         0.0, 1.0, 0.0, 0.0,
         0.0, 0.0, 1.0, 0.0,
         0.0, 0.0, 0.0, 1.0,
};


GLuint base;

void makeRasterFont(Display *dpy)  // Make helvetical font 
{
    XFontStruct *fontInfo;
    Font id;
    unsigned int first, last;

    fontInfo = XLoadQueryFont(dpy,
             "-adobe-helvetica-medium-r-normal--17-120-100-100-p-88-iso8859-1");
    if (fontInfo == NULL) {
        printf ("no font found\n");
        exit (0);
    }
    id = fontInfo->fid;
    first = fontInfo->min_char_or_byte2;
    last = fontInfo->max_char_or_byte2;

    base = glGenLists(last+1);
    if (base == 0) {
        printf ("out of display lists\n");
    exit (0);
    }
    glXUseXFont(id, first, last-first+1, base+first);
}

void printString(char *s)
{
    glPushAttrib (GL_LIST_BIT);
    glListBase(base);
    glCallLists(strlen(s), GL_UNSIGNED_BYTE,
                (unsigned char *)s);
    glPopAttrib ();
}

void drawgrid(struct griddata& den1)
/*  Draws and labels x,y,z,borders of map
*/
{
  int i,pt0[3],pt1[3],pt2[3];
  float fpos0[3],fpos1[3];

  glLineWidth(3);
  glColor3fv(greenv);  //This the simpler, vector version
   for(i=0;i<=2;i++){
     pt0[i]=0;  //origin
     pt1[i]=pt0[i];
     pt2[i]=den1.amax[i]-den1.amin[i]; // far corner
   }
   imcv(den1,pt0,fpos0);
  
   glColor3fv(writecolv);
   glRasterPos3f(fpos0[0],fpos0[1],fpos0[2]);
   printString("O");
   pt1[0]= pt2[0];     
   imcv(den1,pt1,fpos1);
   pt1[0]=pt0[0];
   glRasterPos3f(fpos1[0],fpos1[1],fpos1[2]);
   printString("X");
   pt1[1]= pt2[1];     
   imcv(den1,pt1,fpos1);
   pt1[1]=pt0[1];
   glRasterPos3f(fpos1[0],fpos1[1],fpos1[2]);
   printString("Y");
   pt1[2]= pt2[2];     
   imcv(den1,pt1,fpos1);
   pt1[2]=pt0[2];
   glRasterPos3f(fpos1[0],fpos1[1],fpos1[2]);
   printString("Z");
   

   for(i=0;i<=2;i++){
     pt1[i]= pt2[i];     
     imcv(den1,pt1,fpos1);
     glColor3fv(greenv);  //This the simpler, vector version
     glBegin(GL_LINE_STRIP);
        glVertex3f(fpos0[0],fpos0[1],fpos0[2]);
        glVertex3f(fpos1[0],fpos1[1],fpos1[2]);
     glEnd();
     pt1[i]=pt0[i];
   }

   pt0[1]=0; // y
   pt0[0] = pt2[0];
   pt0[2] = pt2[2];
   for(i=0;i<=2;i++){
     pt1[i]=pt0[i];
   }
   imcv(den1,pt0,fpos0);
   pt1[0]=0;
   imcv(den1,pt1,fpos1);
   glBegin(GL_LINE_STRIP);
        glVertex3f(fpos0[0],fpos0[1],fpos0[2]);
        glVertex3f(fpos1[0],fpos1[1],fpos1[2]);
   glEnd();
   pt1[0]=pt0[0];
   pt1[2]=0;
   imcv(den1,pt1,fpos1);
   glBegin(GL_LINE_STRIP);
        glVertex3f(fpos0[0],fpos0[1],fpos0[2]);
        glVertex3f(fpos1[0],fpos1[1],fpos1[2]);
   glEnd();
   pt1[2]=pt0[2];
   pt1[1]=pt2[1];
   imcv(den1,pt1,fpos1);
   glBegin(GL_LINE_STRIP);
        glVertex3f(fpos0[0],fpos0[1],fpos0[2]);
        glVertex3f(fpos1[0],fpos1[1],fpos1[2]);
   glEnd();

   pt0[2]=0; // z
   pt0[1] = pt2[1];
   pt0[0] = pt2[0];
   for(i=0;i<=2;i++){
     pt1[i]=pt0[i];
   }
   imcv(den1,pt0,fpos0);
   pt1[0]=0;
   imcv(den1,pt1,fpos1);
   glBegin(GL_LINE_STRIP);
        glVertex3f(fpos0[0],fpos0[1],fpos0[2]);
        glVertex3f(fpos1[0],fpos1[1],fpos1[2]);
   glEnd();
   pt1[0]=pt0[0];
   pt1[1]=0;
   imcv(den1,pt1,fpos1);
   glBegin(GL_LINE_STRIP);
        glVertex3f(fpos0[0],fpos0[1],fpos0[2]);
        glVertex3f(fpos1[0],fpos1[1],fpos1[2]);
   glEnd();
   pt1[1]=pt0[1];
   pt1[2]=pt2[2];
   imcv(den1,pt1,fpos1);
   glBegin(GL_LINE_STRIP);
        glVertex3f(fpos0[0],fpos0[1],fpos0[2]);
        glVertex3f(fpos1[0],fpos1[1],fpos1[2]);
   glEnd();

#if 1
   pt0[0]=0; // x
   pt0[1] = pt2[1];
   pt0[2] = pt2[2];
   for(i=0;i<=2;i++){
     pt1[i]=pt0[i];
   }
   imcv(den1,pt0,fpos0);
   pt1[1]=0;
   imcv(den1,pt1,fpos1);
   glBegin(GL_LINE_STRIP);
        glVertex3f(fpos0[0],fpos0[1],fpos0[2]);
        glVertex3f(fpos1[0],fpos1[1],fpos1[2]);
   glEnd();
   pt1[1]=pt0[1];
   pt1[2]=0;
   imcv(den1,pt1,fpos1);
   glBegin(GL_LINE_STRIP);
        glVertex3f(fpos0[0],fpos0[1],fpos0[2]);
        glVertex3f(fpos1[0],fpos1[1],fpos1[2]);
   glEnd();
   pt1[2]=pt0[2];
   pt1[0]=pt2[0];
   imcv(den1,pt1,fpos1);
   glBegin(GL_LINE_STRIP);
        glVertex3f(fpos0[0],fpos0[1],fpos0[2]);
        glVertex3f(fpos1[0],fpos1[1],fpos1[2]);
   glEnd();
#endif
}

Boolean draw_screen(XtPointer clientdata)
// NOTE: clientdata is not used, but needed to compile - for some reason
{
/*  Use cut planes equildistance from origin - 
 sc.cut changes position of planes    = 1.0  maximum depth of field,  0 - minimum
*/
   float eyepos,cutfront,cutback;
   glMatrixMode(GL_PROJECTION);
   glLoadIdentity();
 
   cutfront = 0.001+sc.cmeyedist*(1-sc.cut);
   cutback = sc.cmeyedist+sc.cmeyedist*sc.cut; 
   gluPerspective(sc.degfovy, (float)sc.winwidth/sc.winheight,
         cutfront,cutback);// keep front and back planes equal distance from 0
   GLwDrawingAreaMakeCurrent(sc.drawwind1, sc.glxcontext);
   eyepos = sc.cmeyesep/2.0; /* left window -crosseye- viewpoint at +x pos */
   draw_scene(sc.drawwind1,eyepos);
   GLwDrawingAreaSwapBuffers(sc.drawwind1);
   if(splitscreen==1){ // use standard split screen
      GLwDrawingAreaMakeCurrent(sc.drawwind2, sc.glxcontext);
      eyepos = -sc.cmeyesep/2.0;
      draw_scene(sc.drawwind2,eyepos);
      GLwDrawingAreaSwapBuffers(sc.drawwind2);
   }
   return True;  // draws just once  (called by XtAppAddWorkProc)
}

void drawgeo(struct geometry geo)
/* Draws geo - 
   If geo.fixend=3, specialcase, draw the nterm adnd
      cterm ends separate
*/
{
   int origbeg,origend;
   if(geo.fixend!=3)
       drawpeptide(geo);
   else{ // fixend = 3
      origbeg= geo.nend[0];
      origend=geo.nend[1];
      copyfpostopos(geo,geo.fnum,0);// "0" for nterm moves
      geo.nend[1]= geo.fnum;
      drawpeptide(geo); // draw the nterm end 
      geo.nend[1]=origend;
      copyfpostopos(geo,geo.fnum,1);// "1" for cterm moves
      geo.nend[0]= geo.fnum-1;
      drawpeptide(geo); // draw the nterm end 
      geo.nend[0]=origbeg;
   } // fixend = 3
}

void draw_scene(Widget drawwindow,float eyepos)
/*
NOTE: Order of drawing is important since "pick" in this
   order:  Always draw contour last so do not "hit" it during pick
Order:
   1)  glLoadName(2); 
        drawmodel(sc,maidfile,sc.windowselect);
   2)   glLoadName(3); 
        drawpeptide(geo);
   3)    glLoadName(4);   draw bones

To keep proper perspective etc;
For simplicity - assume all number in unit of cm.
  1)Actual screen cnwidth = 38cm, cmheight =29 cm 
       Total screen pixwidth = 1264 pixels, pixheight =  956 pixels
    Typical eye to screen distance = cmeyedist = 60 cm 
    Typical eye seperation = cmeyesep = 8 cm
  2) If alway keep perspective ratio = winwidht/winheight (pixels)
      The object keeps real shape.
  3) Assume object centered at centerpos[x,y,z] 
  4) Place eye at x= centerposx +-cmeyesep/2.0; y=centerposy; z= centerposz +cmeyedist)
  5) Look at centerpos.
  6)Define cmmaxh = maximum height of object
       angle in perspective = fovy =2 * asin(cmmaxh/2.0/cmeyedist)
 */
 {

    int i;
    glMatrixMode(GL_MODELVIEW);
    glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
    glPushMatrix();
    gluLookAt(eyepos,0.0,sc.cmeyedist,0.0,0.0,0.0,0.0,1.0,0.0);/*Look origin (0,0,0) first 3 positions are
         the eye position, next 3 are point along line of sight, last 3 are 
         direction that corresponds to up direction on screen = +x */
    glTranslatef(sc.trans[0], sc.trans[1],sc.trans[2]);
    glMultMatrixf(objmat);
    glPolygonMode(GL_FRONT_AND_BACK, GL_LINE);
    glTranslatef(-sc.centerpos[0],-sc.centerpos[1],-sc.centerpos[2]); //Place centerpos at origin
    if(maidfile.start>0){  // if just starting - just draw acube
       drawcube();
       copyright();
    }
    else{
        glInitNames(); /* reset pick buffer */
        glPushName(-1); // put -1 on the top of the buffer
        glLoadName(2); /* put 2 on top of buffer, indicates residue */
        drawmodel(geo,sc,maidfile,sc.windowselect);  
        glPopName();
        if(sc.drawgrid==1)
              drawgrid(map[1]);            
        if(sc.drawfit == 1){
           glInitNames(); /* reset pick buffer */
           glPushName(-1); // put -1 on the top of the buffer
            glLoadName(6); /* put 6 on top of buffer, indicates all fit
                 If pickfit =1, then turn on drawgeo, and do not draw fit = geo */
           drawfit(geo,fit,fitnum);
           glPopName();      
        }
        if(sc.drawgeo == 1){    
           glInitNames(); /* reset pick buffer */
           glPushName(-1); // put -1 on the top of the buffer
           glLoadName(3); /* put 3 on top of buffer, indicates geo (peptide) */
           drawgeo(geo);
           glPopName();      
        }
        if(sc.drawrot == 1){
          glInitNames(); /* reset pick buffer */
          glPushName(-1); // put -1 on the top of the buffer
          glLoadName(10); /* put 2 on top of buffer, indicates residue */
          for(i=1;i<=rot[rottype].rotnum;i++){
             glPushName(i); /*put pdbfile # on stack*/
             drawrot(rottype,rotpos[i],side);
             glPopName();
          }
          glPopName();      
        }  
        glInitNames(); /* reset pick buffer */
        glPushName(-1); // put -1 on the top of the buffer
        glLoadName(5); /* put 5 on top of buffer, indicates sphere */
        drawspheres();
        glPopName();      
        if(sc.printres == 1){
           labelres();
        }
        glInitNames(); /* reset pick buffer */
        glPushName(-1); // put -1 on the top of the buffer
        glLoadName(4); /* put 4 on top of buffer, indicates bone maxima (X) */
        for(i=1;i<= maidfile.nden;i++){
#if 0
           if(core[i].iexist == 2){
               glPushName(i); /*put pdbfile # on stack*/
               drawcore(map[i],core[i]);
               glPopName();
           }
#endif
           if( (bone[i].iexist == 2)||(bone[i].iexist == 4) ){// =4 new minden picked
               glPushName(i); /*put bonefile # on stack*/
               drawbone(map[i],bone[i],sc.windowselect);
               glPopName();
           }
        }// i
        glPopName();      
        if(sc.drawcontour ==1){
           for(i=1;i<= maidfile.nden;i++){
              if(map[i].contonoff ==1)
                 drawsurface(map[i]);
           }
        }
        if(sc.draweditbones ==1) {
          glInitNames(); /* reset pick buffer */
          glPushName(-1); // put -1 on the top of the buffer
          glLoadName(8); /* indicates drawing edited bones */
          draweditbone(ext);
          glPopName();
        }
    }// else 
    glPopMatrix();
    glFlush();
}

void expose(Widget w, XtPointer client_data, XtPointer call) 
{
    GLwDrawingAreaCallbackStruct *call_data;
    call_data = (GLwDrawingAreaCallbackStruct *) call;
    sc.winwidth=call_data->width;
    sc.winheight=call_data->height;
    glViewport(0, 0, sc.winwidth, sc.winheight);
    draw_screen(client_data);
}
  
void resize(Widget w, XtPointer client_data, XtPointer call) 
{
   GLwDrawingAreaCallbackStruct *call_data;
   call_data = (GLwDrawingAreaCallbackStruct *) call;
   sc.winwidth=call_data->width;
   sc.winheight=call_data->height;
   glViewport(0, 0, sc.winwidth, sc.winheight);
   draw_screen(client_data);
}

void update_view(int mstate, int ox, int nx, int oy, int ny) {
//New Version: just push and release a,s,d,q,and e keys (do not hold down)
XtPointer client_data;
    int dx = ox - nx;
//printf("ox = %d nx = %d dx = %d\n",ox,nx,dx);
    int dy = ny - oy;
    float transmult = sc.cmmaxh/5.0;//scaling factor to use when translating
//    if( (mstate==2)||(mstate==3)||(mstate==4)||(mstate==1)||(mstate==7) ){
    if( (mstate==2)||(mstate==3)||(mstate==4)||(mstate==1) ){
       glPushMatrix();
       glLoadIdentity();
    }
    switch(mstate) {
    case 1:	/* Use right mouse button to rotate along x,y */
	sc.rot[0] = (dy * 180.0)/500.;
	sc.rot[1] = -(dx * 180.0)/500.;
        glRotatef( sc.rot[0],1.0,0.0,0.0); // rotate - angle in degrees
        glRotatef( sc.rot[1],0.0,1.0,0.0); // rotate - angle in degrees
	break;
    case 2:	/*letter s pushed - rotate about x */
	sc.rot[0] = (dx * 180.0)/500.;
        glRotatef( sc.rot[0],1.0,0.0,0.0); // rotate - angle in degrees
	break;
    case 3:	/* rotate about y*/
	sc.rot[1] = (dx * 180.0)/500.;
        glRotatef( sc.rot[1],0.0,1.0,0.0);
	break;
    case 4:	/* rotate about z*/
	sc.rot[2] = (dx * 180.0)/500.;
        glRotatef( sc.rot[2],0.0,0.0,1.0);
	break;
    case 5:	/* Change depth of field*/
	sc.cut = sc.cut+ 0.005*dx;
        if(sc.cut >1.0)
          sc.cut = 1.0;
        else if(sc.cut < 0.01)
           sc.cut = 0.01;
//printf("sc.cut = %f\n",sc.cut);
	break;
    case 6:	/* translate along y */
	sc.trans[1] -= transmult*dx/100.;
	break;
    case 7:	/* Zoom - change perspective angle */
	sc.mag -= dx/500.; // slow down
        sc.degfovy = sc.mag*(180.0/PI)*2.0*asin(sc.cmmaxh/(2.0*sc.cmeyedist));
        if( sc.degfovy > 170.0)
            sc.degfovy=170.0;
        else if( sc.degfovy < 1.0)
            sc.degfovy=1.0;
	break;
    case 8:	/* Use middle mouse button to translate along x,y */
	sc.trans[0] -= transmult*dx/100.;
	sc.trans[1] -= transmult*dy/100.;
	break;
    case 10:	/* zoom */
	sc.trans[2] -= (dx+dy)/100.;
	break;
   case 11:	/* change z of sphere center (x movment) and radius (y movement*/
	sp[sc.isphere].cenpos[2] = sp[sc.isphere].cenpos[2] + transmult*dx/100;
	sp[sc.isphere].radius = sp[sc.isphere].radius + 0.05*dy;
        if(sp[sc.isphere].radius<0.0) 
               sp[sc.isphere].radius=0.0;
	break;
    case 18:	/* change x,y of sphere center using middle button*/
	sp[sc.isphere].cenpos[0] = sp[sc.isphere].cenpos[0] - transmult*dx/100;
	sp[sc.isphere].cenpos[1] = sp[sc.isphere].cenpos[1] - transmult*dy/100;
	break;
    }
    if( (mstate==2)||(mstate==3)||(mstate==4)||(mstate==1) ){
      glMultMatrixf(objmat);
      glGetFloatv(GL_MODELVIEW_MATRIX,objmat);
      glPopMatrix();
    }
    draw_screen(client_data);
}

void input(Widget w, XtPointer client_data, XtPointer call) //left window
{
   char buf[31];
   KeySym keysym;
   XEvent *event = ((GLwDrawingAreaCallbackStruct *) call)->event;
   static int mstate, omx, omy, mx, my;

   switch(event->type) {
   case KeyRelease:
      XLookupString(&event->xkey, buf, sizeof buf, &keysym, NULL);
      switch(keysym) {
      case XK_Escape:
	 exit(EXIT_SUCCESS);
	 break;
      default: break;
      }
      break;

    case KeyPress:  /* Use A key equivalent to button 1*/
      XLookupString(&event->xkey, buf, sizeof buf, &keysym, NULL);
      switch(keysym) { /* key a,s,d = rotate;q,w,e = translate */ 
      case XK_a:
 	    mstate = 2;
 	    break;
      case XK_s:
 	    mstate = 3;
 	    break;
      case XK_d:
 	    mstate = 4;
 	    break;
      case XK_q:  // change depth of field (sc.cut)
            mstate=5;
 	    break;
      case XK_e:  //rsize
 	    mstate = 7;
 	    break;
      default: break;
      }
      break; 
       
     case ButtonPress:
        if ( (event->xbutton.button == Button3)) {
	    if(mstate ==0){//If right button and no letter - rotate
 	       if(sc.movesphere ==1)// change z pos
                  mstate = 11;
               else // standard - change depth of field
                  mstate = 1;
	       mx = event->xbutton.x;
	       my = event->xbutton.y;
            }
            else 
	       mx = event->xbutton.x;
        }
        else if (event->xbutton.button == Button2) { // If push middle button and no letter pushed - translate
	    if(mstate ==0){
 	       if(sc.movesphere ==1)// change z pos
                  mstate = 18;
               else // standard - change depth of field
                  mstate = 8;
	       mx = event->xbutton.x;
	       my = event->xbutton.y;
            }
	}
        else if(event->xbutton.button == Button1) {//Start picking routine
           mstate =0;  // This keeps from going to update_view
	   mx =event->xbutton.x;
	   my = event->xbutton.y;
           pick(mx,my);
        }
	break;
    case ButtonRelease:
        if ( (event->xbutton.button == Button2)|| (event->xbutton.button == Button3))
	    mstate = 0;
	break;
    case MotionNotify:
	if (mstate) {
	    omx = mx;
	    omy = my;
	    mx =event->xbutton.x;
	    my = event->xbutton.y;
//printf("omx = %d mx = %d\n",omx,mx);
	    update_view(mstate, omx,mx,omy,my);
	}
	break;
   }
}




#define BUFSIZE 512

void pick(int mx,int my)
{
   GLuint selectBuf[BUFSIZE];// unsigned integer
   GLint hits,viewport[4];
   float eyepos,cutfront,cutback; 

   GLwDrawingAreaMakeCurrent(sc.drawwind1, sc.glxcontext);
   glGetIntegerv(GL_VIEWPORT,viewport);
   glSelectBuffer(BUFSIZE,selectBuf);   
   (void) glRenderMode(GL_SELECT); // Puts into select mode
   glInitNames();// clears the name stack 
//   glPushName(-1); // put -1 on the top of the buffer
   glPushName(0); // put -1 on the top of the buffer
      // drawing routines use either pushname (add name to buffer)or
      //  loadname  (replace top name by new name) 
      //  These calls ignored if not in selection mode

   glMatrixMode(GL_PROJECTION);
   glPushMatrix();
   glLoadIdentity();
   gluPickMatrix( (GLdouble)mx,(GLdouble)(viewport[3]-my),
         5.0,5.0,viewport);//create 5x5 pixel picking region near cursor location
// This section is copied directly from drawscreen
   cutfront = 0.001+sc.cmeyedist*(1-sc.cut);
   cutback = sc.cmeyedist+sc.cmeyedist*sc.cut; 
   gluPerspective(sc.degfovy, (float)sc.winwidth/sc.winheight,
         cutfront,cutback);// keep front and back planes equal distance from 0
   eyepos = sc.cmeyesep/2.0; /* left window -crosseye- viewpoint at +x pos */
    draw_scene(sc.drawwind1,eyepos);
    hits = glRenderMode(GL_RENDER); // Puts back into (default) feedback mode
    if(hits>0)
       processhits(hits,selectBuf);
}         

static void pickagain_pushed(Widget w,XtPointer  client_data,XtPointer  call_data)
// Quit fuse procedure
{
    geo.fuse[1] = geo.fuse[2]= -1;  // reset to defaault vale
    geo.fuse[0]=0;
    Widget shell = (Widget) client_data;
    XtDestroyWidget (shell);
}

static void fuse_cancel_pushed(Widget w,XtPointer  client_data,XtPointer  call_data)
// Quit fuse procedure
{
    geo.fuse[0] = geo.fuse[1] = geo.fuse[2]= -1;  // reset to defaault vale
    Widget shell = (Widget) client_data;
    XtDestroyWidget (shell);
}

void fuse_dialog(Widget w)  // Dialog for inputting fit #
/* DO NOT UNDERSTAND -  the second ("clear") button cannot be
    be given an arbitrary function, apparently, it must be
    a clear type button
*/
{
    Widget dialog, rc;
    XmString string;
    Arg args[5];
    int n = 0;
    static ActionAreaItem action_items[] = {
        { "OK",     fuse_ok_pushed,    NULL         },
        { "PICK AGAIN",    pickagain_pushed ,NULL         },
        { "QUIT", fuse_cancel_pushed,  NULL          },
        { "Help",   help,          "Help Button" },
    };
    /* Create a TemplateDialog that will contain the control area
     * and the action area buttons for the dialog
     */ 
    string = XmStringCreateLocalized ("Dialog Shell");
    XtSetArg (args[n], XmNdialogTitle, string); n++;
    XtSetArg (args[n], XmNautoUnmanage, False); n++;
#if 0 // tried different things to position window - NONE WORKED
//    XtSetArg (args[n],  XmNtopAttachment, 300);n++;
    XtSetArg (args[n],   XmNtopPosition,      300  ); n++;
//    XtSetArg (args[n],  XmNleftAttachment,   XmATTACH_WIDGET); n++;
#endif
    dialog = XmCreateTemplateDialog (toplevel, "dialog", args, n);
    XmStringFree (string);

    /* now that the dialog is created, set the Cancel button's
     * client data, so close_dialog() will know what to destroy.
     */
    action_items[0].data=action_items[1].data = action_items[2].data = (XtPointer) dialog;

    /* create the control area which contains a
     * Label gadget and a TextField widget.
     */
    rc = XtVaCreateWidget ("control_area", 
        xmRowColumnWidgetClass, dialog,
        NULL);
    string = XmStringCreateLocalized ("FUSE THE 2 PICKED FITS");
    XtVaCreateManagedWidget ("label", xmLabelGadgetClass, rc,
        XmNlabelString,    string,
        NULL);
    XmStringFree (string);
   /* RowColumn is full -- now manage */
    XtManageChild (rc);
  
    /* Create the action area. */
    CreateActionArea (dialog, action_items, XtNumber (action_items));

    XtManageChild (dialog);
    XtPopup (XtParent (dialog), XtGrabNone);
    glClearColor(backcolv[0],backcolv[1],backcolv[2],0);//background color :red, green, blue, and alpha = black
}

void assignclose_dialog(Widget w,XtPointer  client_data,XtPointer  call_data)
{
    geo.resnum = -1; // this will turn off assignseq
    Widget shell = (Widget) client_data;
    XtDestroyWidget (shell);
}

static void assign_ok_pushed(Widget w,XtPointer  client_data,XtPointer  call_data)
{
    char* text[5];
    text[0] =   XmTextFieldGetString (ss.w[0]);
    geo.resnum = atoi(text[0]);
    printf(" Input seq num = %d\n",geo.resnum);
    if(geo.resnum>geo.seqnum){
         errormessage(toplevel,"ERROR: sequence number greater than");
         geo.resnum = -1; // this turns off assign
    }
     else
        buildseq(map[1],geo,rot);
//    sc.drawgeo = 1;
    Widget shell = (Widget) client_data;
    XtDestroyWidget (shell);
}


void assign_dialog(Widget w)  // Dialog for inputting sequence aa #
{
    Widget dialog, rc;
    XmString string;
    Arg args[5];
    int n = 0;
    static ActionAreaItem action_items[] = {
        { "OK",     assign_ok_pushed,     NULL          },
        { "Clear",  clear_pushed,  NULL          },
        { "Cancel", assignclose_dialog,  NULL          },
        { "Help",   help,          "Help Button" },
    };
    /* Create a TemplateDialog that will contain the control area
     * and the action area buttons for the dialog
     */ 
    string = XmStringCreateLocalized ("Dialog Shell");
    XtSetArg (args[n], XmNdialogTitle, string); n++;
    XtSetArg (args[n], XmNautoUnmanage, False); n++;
    dialog = XmCreateTemplateDialog (toplevel, "dialog", args, n);
    XmStringFree (string);

    /* now that the dialog is created, set the Cancel button's
     * client data, so close_dialog() will know what to destroy.
     */
    action_items[0].data=action_items[2].data = (XtPointer) dialog;

    /* create the control area which contains a
     * Label gadget and a TextField widget.
     */
    rc = XtVaCreateWidget ("control_area", 
        xmRowColumnWidgetClass, dialog, NULL);
    ss.num = 1;
    string = XmStringCreateLocalized ("ENTER SEQUENCE # FOR THIS RESIDUE:");
    XtVaCreateManagedWidget ("label", xmLabelGadgetClass, rc,
        XmNlabelString,    string,
        NULL);
    XmStringFree (string);
    ss.w[0] = XtVaCreateManagedWidget ("text-field",xmTextFieldWidgetClass, rc, NULL);
   /* RowColumn is full -- now manage */
    XtManageChild (rc);
  
    /* Create the action area. */
    CreateActionArea (dialog, action_items, XtNumber (action_items));

    /* callback for Return in TextField.  Use dialog as client data */
    XtAddCallback (ss.w[0], XmNactivateCallback,activate_cb, dialog);

    XtManageChild (dialog);
    XtPopup (XtParent (dialog), XtGrabNone);
    glClearColor(backcolv[0],backcolv[1],backcolv[2],0);//background color :red, green, blue, and alpha = black
}

void removecancel_dialog(Widget w,XtPointer  client_data,XtPointer  call_data)
{
    geo.removefit = -1; // this will turn off remove
    Widget shell = (Widget) client_data;
    XtDestroyWidget (shell);
}

static void remove_ok_pushed(Widget w,XtPointer  client_data,XtPointer  call_data)
{
    char* text[5];

    removefit(fit,geo.pickfitnum);
    Widget shell = (Widget) client_data;
    XtDestroyWidget (shell);
}

void removefit_dialog(Widget w)  // This is overkill, only need yes no
{
    Widget dialog, rc;
    XmString string;
    Arg args[5];
    int n = 0;
    static ActionAreaItem action_items[] = {
        { "OK",     remove_ok_pushed,     NULL          },
        { "Clear",  clear_pushed,  NULL          },
        { "Cancel", removecancel_dialog,  NULL          },
        { "Help",   help,          "Help Button" },
    };
    /* Create a TemplateDialog that will contain the control area
     * and the action area buttons for the dialog
     */ 
    string = XmStringCreateLocalized ("Dialog Shell");
    XtSetArg (args[n], XmNdialogTitle, string); n++;
    XtSetArg (args[n], XmNautoUnmanage, False); n++;
    dialog = XmCreateTemplateDialog (toplevel, "dialog", args, n);
    XmStringFree (string);

    /* now that the dialog is created, set the Cancel button's
     * client data, so close_dialog() will know what to destroy.
     */
    action_items[0].data=action_items[2].data = (XtPointer) dialog;

    /* create the control area which contains a
     * Label gadget and a TextField widget.
     */
    rc = XtVaCreateWidget ("control_area", 
        xmRowColumnWidgetClass, dialog, NULL);
    ss.num = 1;
    string = XmStringCreateLocalized ("PRESS OK TO REMOVE THIS FIT, CANCEL TO QUIT:");
    XtVaCreateManagedWidget ("label", xmLabelGadgetClass, rc,
        XmNlabelString,    string,
        NULL);
    XmStringFree (string);
    ss.w[0] = XtVaCreateManagedWidget ("text-field",xmTextFieldWidgetClass, rc, NULL);
   /* RowColumn is full -- now manage */
    XtManageChild (rc);
  
    /* Create the action area. */
    CreateActionArea (dialog, action_items, XtNumber (action_items));

    /* callback for Return in TextField.  Use dialog as client data */
    XtAddCallback (ss.w[0], XmNactivateCallback,activate_cb, dialog);

    XtManageChild (dialog);
    XtPopup (XtParent (dialog), XtGrabNone);
    glClearColor(backcolv[0],backcolv[1],backcolv[2],0);//background color :red, green, blue, and alpha = black
}

void processhits(GLint hits,GLuint buffer[])
/*  Buffer consists of a set of data for each hit:
   1) number of names for this hit
   2),3) min and max planes for the hit
   4) 5)....  the names for this hit
I will use only the last hit:
If pick residue:
hits = number of hits
buffer[0]= # of names 
buffer[3] = pick type = 2 for residue, 4 for core maxima (X)
buffer[4] = pdb file number (1,2,3,...)
buffer[5] = segment number
buffer[6] = residue number in that segment
*/
{ 
XtPointer client_data;
   int i,pdbn,segn,resn,pnum,ivec[3];
   int bonen,cenpos;
   Widget widg,widget;
   float fvec[3];

   if(verbose>=2){
      printf("number of hits = %d\n",hits);
      if(hits>0){
          printf("number of names = %d\n",buffer[0]);
          for(i = 1;i<=buffer[0];i++)
             printf("name =%d\n",buffer[2+i]);
      }
   }
   if(hits>0){
      if(verbose>=1)
         printf("buffer3 = %d  buffer4 = %d buffer5 = %d\n",buffer[3],buffer[4],buffer[5]);
     if( (buffer[3]==6)&&(geo.removefit==1) ){ // if picked fit and want to remove fit
         geo.pickfitnum = buffer[4];
         removefit_dialog(toplevel);
     }
     if( (sc.centerres == 1)&&(buffer[3]==6) ){ // Center on this last picked fit
          for(i=0;i<=2;i++)
             sc.centerpos[i] = fit[ buffer[4] ].pos[ buffer[5] ][1][i];
          for(i=0;i<=2;i++)
             sc.trans[i]=0.0;
       if(sc.pdbwindow==1){ // Then window about centerpos
          findpdbwind();
          findbonewind();
          sc.windowselect =1;
       }
       else if (sc.pdbwindow !=1)
           sc.windowselect = -1;
       for(i=1;i<= maidfile.nden;i++){ // contour about picked fit
            if(map[i].contonoff ==1)
                contourmap(map[i],fit[ buffer[4] ].pos[ buffer[5] ][1]);
            sc.drawcontour = 1; // Turn on drawing of contour
       } 
     }
     if( (buffer[3]==6)&&(fit[ buffer[4] ].rankfit[0]>0) ) // if picked fit and fit part of rankset
        printf("Fit %d  Number %d\n",
            fit[ buffer[4] ].rankfit[0], fit[ buffer[4] ].rankfit[1]);
     if( (geo.fuse[0]>=0) && (geo.fuse[0]<2)&& (buffer[3]==6) ){ // picked fit and want to fuse to fits
       if( (buffer[4]!=geo.fuse[1])&&(buffer[4]!=geo.fuse[2]) ){// not equal to previous picked
         if( buffer[5] < fit[ buffer[4] ].nres/2.0 ){ // then nterm end
           if(geo.fuse[1]>=0){ // already assigned
              geo.fuse[0]= -1;
              geo.fuse[1]=geo.fuse[2]= -1;
              errormessage(toplevel,"ERROR: BOTH PICKS ARE AT N TERMINAL - QUITTING");
           }
           else{ // OK
              geo.fuse[0]=geo.fuse[0]+1;
              geo.fuse[1]= buffer[4]; // Number of nterm fit
           }
         }
         else{ // cterm
           if(geo.fuse[2]>=0){ // already assigned
              geo.fuse[0]= -1;
              geo.fuse[1]=geo.fuse[2]= -1;
              errormessage(toplevel,"ERROR: BOTH PICKS ARE AT C TERMINAL - QUITTING");
           }
           else{ // OK
              geo.fuse[0]=geo.fuse[0]+1;
              geo.fuse[2]= buffer[4]; // Number of nterm fit
           }
         }
         draw_screen(client_data);// Must redraw to draw picked in different color
         if(geo.fuse[0]==2){ 
           fuse_dialog(toplevel);
         } // geo.fuse[0]==2
       } // not equal
     } 
     else if( (geo.resnum== -1000)&& (buffer[3]==6) ){ // picked fit and want to assign sequence 
         geo.pickfitnum = buffer[4];
         geo.fnum = buffer[5];         
         assign_dialog(toplevel); // sets geo.resnum = sequence number for this fit res.
#if 0
printf("AT 2  geo.resnum = %d seqnum = %d\n",geo.resnum,geo.seqnum);
         if(geo.resnum>0){// Indicates did not cancel
              if(geo.resnum>geo.seqnum){
                errormessage(toplevel,"ERROR: sequence number greater than");
                geo.resnum = -1; // this turns off assign
              }
              else
                buildseq(map[1],geo,rot);
         } // resnum >0
#endif
     }
     else if( (sc.pickfit==1)&&(buffer[3]==6) ){ // picked fit
       geo.pickfitnum = buffer[4];
       copyfittogeo(geo,fit[ geo.pickfitnum ]);
       if(geo.badfit==0) // fit is  too long
          errormessage(toplevel,"ERROR: THIS FIT IS TOO LONG TO MAKE 'ACTIVE' - CUT FIT INTO 40 RESIDUE SECTIONS or increase nafit");
       else{
        if (widget = XtNameToWidget (fitmenu,"button_12")) //  Turn off check box for drawing picked fit (geo)
            XtVaSetValues (widget, XmNset,TRUE, NULL);
        sc.drawgeo=1;
#if 0 // run aveden
float avecb2,aveco,finalden,phi,psi;
float psimin,psimax,phimin,phimax;
int resbad[NAFIT],nbad,rtype;
       avedenreg(map[1],geo,geo.pos,&finalden,&avecb2,&nbad,
                  geo.nend[1]-1,geo.nend[1]); // Find average density for last cterm residue additions
//       aveden(map[1],geo,geo.pos,&finalden,&avecb2,&aveco,&nbad,resbad); // Find average density (ave) for initial orientation 
       for(i=geo.nend[0];i<=geo.nend[1]-2;i++){
          findphipsi(geo.pos,i,&phi,&psi);
          rtype = geo.rtype[ i+1 ];
          if(checkangrange(geo,phi,psi,rtype,&phimin,&phimax,&psimin,&psimax)==0)
              printf(" PHI= %6.3f/PSI=%6.3f GREATER THAN LIMIT (2PI = 360)\n",
                phi*180.0/PI,psi*180.0/PI);
       }
       printf("In picked fit:  finalden = %5.3f avecb2 = %5.3f avecb = %5.3f nbad = %d \n",
             finalden,avecb2,aveco,nbad);
#endif
#if 0
        sc.pickfit = -1;
        if (widg = XtNameToWidget (fitmenu,"button_10")) //  Turn off check box for picking fit
            XtVaSetValues (widg, XmNset,FALSE, NULL);
        if (widg = XtNameToWidget (fitmenu,"button_11")) //  Turn off check box for picking fit
            XtVaSetValues (widg, XmNset,TRUE, NULL);
#endif
      } // fit is not to long
     }  // picked fit
     else if(buffer[3]==2){ // picked residue
       pdbn= buffer[4];
       segn= buffer[5];
       resn= buffer[6];
       pnum=  ++maidfile.pdb[pdbn].pickres[0].num;
       if(pnum>PDBLABEL-5)
         errormessage(toplevel,"You are near allowed limit of pdb labels.. 'Delete Label' and continue");
       maidfile.pdb[pdbn].pickres[pnum].seg = segn;
       maidfile.pdb[pdbn].pickres[pnum].num = resn;
       if(sc.addsphere == 1){
          maidfile.nsphere = maidfile.nsphere+1;
          sc.isphere = maidfile.nsphere; // Make new sphere active
          for(i=0;i<=2;i++)
             sp[sc.isphere].cenpos[i] = maidfile.pdb[pdbn].mainch[segn][resn][1][i];
          sc.addsphere = -1; // turn off adding of sphere
          if (widg = XtNameToWidget (spheremenu,"button_1")) //DEFAULT:  Turn off check box for label on/off button
              XtVaSetValues (widg, XmNset,FALSE, NULL);
          sc.drawspheres = 1;// turn on drawing of spheres
          if (widg = XtNameToWidget (spheremenu,"button_0")) //DEFAULT:  Turn off check box for label on/off button
              XtVaSetValues (widg, XmNset,TRUE, NULL);
printf("isphere = %d nsphere= %d sphere center = %f %f %f\n",sc.isphere,maidfile.nsphere,sp[sc.isphere].cenpos[0],sp[sc.isphere].cenpos[1],sp[sc.isphere].cenpos[2]);
          goto endrespick;
       }
       if(sc.centerres == 1){ // Center on this last picked residue
          for(i=0;i<=2;i++)
             sc.centerpos[i] = maidfile.pdb[pdbn].mainch[segn][resn][1][i];
          for(i=0;i<=2;i++)
             sc.trans[i]=0.0;
       }
       if(sc.pdbwindow==1){ // Then window about centerpos
          findpdbwind();
          findbonewind();
          sc.windowselect =1;
       }
       else if (sc.pdbwindow !=1)
           sc.windowselect = -1;
       for(i=1;i<= maidfile.nden;i++){
            if(map[i].contonoff ==1)
                contourmap(map[i],maidfile.pdb[pdbn].mainch[segn][resn][1]);
            sc.drawcontour = 1; // Turn on drawing of contour
       } 
       if(sc.refinepdbside == 1){
          fitnum[0]=1;
          fitnum[1] = 1; // renumber fits so run consecuatively
          allocatefit(fit[1],10);
          refineonepdbside(0,map[1], geo,rot,maidfile.pdb[pdbn],segn,
                   resn); // 0 means to not write data to file
#if 0
          copyreglevtofit(fit[1],maidfile.pdb[pdbn],segn,resn-2,resn+2);          
          copyonefittolev(fit[1],maidfile.pdb[pdbn],2,segn,resn);
#endif          
       } 
printf(" total picked = %d, pdbn = %d segn = %d resn = %d\n",pnum,pdbn,segn+1,resn); // add 1 to segn because seg starts with 0
       endrespick:;   
     } // picked residue
     else if( (geo.fuse[0]<0) && (buffer[3]==3) ){ // picked geo
        geo.fnum = geo.nend[0]+buffer[4]; // geo number of picked residue
        printf(" Setting geo.num = %d, geo.nend[0,1]= %d %d\n",geo.fnum,geo.nend[0],geo.nend[1]);
        if(sc.centerres == 1){ // Center on this last picked residue
          for(i=0;i<=2;i++)
             sc.centerpos[i] = geo.pos[ geo.nend[0]+buffer[4] ][1][i];
          for(i=0;i<=2;i++)
             sc.trans[i]=0.0;
        }
        if(sc.pdbwindow==1){ // Then window about centerpos
          findpdbwind();
          findbonewind();
          sc.windowselect =1;
        }
        else if (sc.pdbwindow !=1)
           sc.windowselect = -1;
        for(i=1;i<= maidfile.nden;i++){
            if(map[i].contonoff ==1)
                contourmap(map[i],geo.pos[ geo.nend[0]+buffer[4] ][1]);
            sc.drawcontour = 1; // Turn on drawing of contour
        } 
      }
      else if(buffer[3]==4){ // Picked bone pt
       bonen = buffer[4];
       cenpos= buffer[5];
       d1tod3(bone[bonen].cmax,bone[bonen].ptv[cenpos].pos,ivec);
       imcv(map[bonen],ivec,fvec);// convert int grid ivec to float world co-ord fvec
       if(sc.addsphere == 1){
          maidfile.nsphere = maidfile.nsphere+1;
          sc.isphere = maidfile.nsphere; // Make new sphere active
          for(i=0;i<=2;i++)
             sp[sc.isphere].cenpos[i] = fvec[i];
          sc.addsphere = -1; // turn off adding of sphere
          if (widg = XtNameToWidget (spheremenu,"button_1")) //DEFAULT:  Turn off check box for label on/off button
              XtVaSetValues (widg, XmNset,FALSE, NULL);
          sc.drawspheres = 1;// turn on drawing of spheres
          if (widg = XtNameToWidget (spheremenu,"button_0")) //DEFAULT:  Turn off check box for label on/off button
              XtVaSetValues (widg, XmNset,TRUE, NULL);
printf("isphere = %d nsphere= %d sphere center = %f %f %f\n",sc.isphere,maidfile.nsphere,sp[sc.isphere].cenpos[0],sp[sc.isphere].cenpos[1],sp[sc.isphere].cenpos[2]);
          goto endbonepick;
       }
       if(sc.centerres == 1){ // Center on this last picked bone
          for(i=0;i<=2;i++)
             sc.centerpos[i] =fvec[i];
          for(i=0;i<=2;i++)
             sc.trans[i]=0.0;
       }
       if(sc.pdbwindow==1){ // Then window about centerpos
          findpdbwind();
          findbonewind();
          sc.windowselect =1;
       }
       else if (sc.pdbwindow !=1)
           sc.windowselect = -1;
       for(i=1;i<= maidfile.nden;i++){
            if(map[i].contonoff ==1)
                contourmap(map[i],fvec);
            sc.drawcontour = 1; // Turn on drawing of contour
       }  
       endbonepick:;
#if 1
int num,ivec[3];
num =  bone[bonen].ptv[cenpos].pos;
d1tod3(bone[bonen].cmax,num,ivec); // convert from 1d to 3d vec[i,j,k]
printf("bone num = %d  Center num = %d: (%d %d %d) num of con. = %d:",
      bonen,cenpos,ivec[0],ivec[1],ivec[2],bone[bonen].ptv[cenpos].num);
    for(i=0;i<bone[bonen].ptv[cenpos].num;i++)
       printf("  %d",bone[bonen].ptv[cenpos].con[i]);
printf("\n");
#endif
     }  // picked bone       
     else if(buffer[3]==5){ // Picked sphere
printf("buffer[4] = %d\n",buffer[4]);
        if(sc.selectsphere==1){
           sc.selectsphere = -1;
           sc.isphere = buffer[4];
           if (widg = XtNameToWidget (spheremenu,"button_3")) //DEFAULT:  Turn off check box for label on/off button
              XtVaSetValues (widg, XmNset,FALSE, NULL);
        }
        else if(sc.deletesphere==1){// delete picked sphere = buffer[4]
          for(i=buffer[4]+1;i<=maidfile.nsphere;i++){
#if 0
            for(k=0;k<=2;k++)
               sp[i].cenpos[k] = sp[i-1].cenpos[k];
            sp[i].radius= sp[i-1].radius;
#endif
            sp[i-1]=sp[i];
          }
          maidfile.nsphere= maidfile.nsphere-1;
          sc.deletesphere = -1;
          if (widg = XtNameToWidget (spheremenu,"button_4")) //DEFAULT:  Turn off check box for label on/off button
              XtVaSetValues (widg, XmNset,FALSE, NULL); 
        }              
     } // selected sphere
     draw_screen(client_data);// Must redraw to draw label
   } // hits>0
   quitrout:;
}

void findbone1(struct griddata& den1,struct bonedata& bone1)
{
   int totnum,i,k,ivec[3],tcount;
   float fvec[3],dist2;

   if(verbose>=2)
      printf("Finding bones in window  #\n");
   totnum = bone1.ptv.end()-bone1.ptv.begin();// total number bones 
   tcount = 0;  // total in window
   for(i=0;i<totnum;++i){
     d1tod3(bone1.cmax,bone1.ptv[i].pos,ivec);
     imcv(den1,ivec,fvec);// convert int grid ivec to float world co-ord fvec
     dist2 = 0.0;
     for(k = 0;k<=2;k++)
        dist2 = dist2 + sqr(fvec[k] -sc.centerpos[k]);
     if(dist2 < sc.boneradius2){
        tcount = tcount+1;
        if(tcount > NBONEWIN -2){
            errormessage(toplevel,"WARNING: Too many bone pts in window");
            tcount = tcount -1;
            goto endrout; // This clips # bone pts in window at max allowed
        }
        else 
          bone1.winpts[tcount] = i;
     } // dist2 < 
   }// i
   endrout:;
   bone1.winpts[0] = tcount;
} 

void findbonewind()
/*
   Go through list of possible bone file, it it exists,
     then bind list of pts in window
*/
{
    int ib;

   if(verbose>=2)
       printf("Finding bones in window \n");
    for(ib=1;ib<NMAPS;ib++){ // each potential bone file
      if(verbose>=2)
         printf("ib = %d  bone[ib].iexist = %d\n",ib,bone[ib].iexist);
      if( (bone[ib].iexist == 2)||(bone[ib].iexist == 4) )// =4 new minden picked
         findbone1(map[ib],bone[ib]); 
    }
}

void findpdbwind()
// For each residue - check if within sc.pdbradius of centerpos -
//    if it is - add it to maidfile.pdb[i].windmod[segnum][resnum]
{
   int ip,is,ires,i,ii;
   float dist2,xyz;

   for(ip=1;ip<=maidfile.npdb;ip++){ // each pdb file
      for(is=0;is<maidfile.pdb[ip].resdat[0][0][0];is++){ // all segments
	 ii=maidfile.pdb[ip].windmod[is][0]=0;
	 for(ires=1;ires<=maidfile.pdb[ip].resdat[is+1][0][0];ires++) { //all residues
	    dist2 = 0.0;
            for(i = 0;i<=2;i++){
	       xyz= maidfile.pdb[ip].mainch[is][ires][1][i]; // Ca positon ("1")
               dist2 = dist2 + sqr(xyz -sc.centerpos[i]);
            }
            if(dist2 < sc.pdbradius2){
	          ii = maidfile.pdb[ip].windmod[is][0] =maidfile.pdb[ip].windmod[is][0]+1;
	          maidfile.pdb[ip].windmod[is][ii]= ires;
            }
          } // ires
	  if(ii > NRESWIND-2){
	    errormessage(toplevel,"WARNING: Number of residues in Window > NRESWIND");
	    maidfile.pdb[ip].windmod[is][0]=NRESWIND -2; // Clip number residues in window
          }
       } // is
    } // ip
printf(" In findpdb:  # res in 1 = %d  in 2 = %d\n", maidfile.pdb[1].windmod[0][0],
                     maidfile.pdb[2].windmod[0][0]);
}



void input2(Widget w, XtPointer client_data, XtPointer call) //right window
// Input from right window - same as input exceept add screewidth
//    to x mouse position
{
   char buf[31];
   KeySym keysym;
   XEvent *event = ((GLwDrawingAreaCallbackStruct *) call)->event;
   static int mstate, omx, omy, mx, my;

   switch(event->type) {
   case KeyRelease:
      XLookupString(&event->xkey, buf, sizeof buf, &keysym, NULL);
      switch(keysym) {
      case XK_Escape:
	 exit(EXIT_SUCCESS);
	 break;
      case XK_a:
	  mstate = 0;
         break;
      case XK_s:
	  mstate = 0;
         break;
      case XK_d:
	  mstate = 0;
         break;
      case XK_q:
	  mstate = 0;
         break;
      case XK_w:
	  mstate = 0;
         break;
      case XK_e:
	  mstate = 0;
         break;
      default: break;
      }
      break;

    case KeyPress:  /* Use A key equivalent to button 1*/
      XLookupString(&event->xkey, buf, sizeof buf, &keysym, NULL);
      switch(keysym) { /* key a,s,d = rotate;q,w,e = translate */ 
      case XK_a:
 	    mstate = 2;
	    mx = event->xbutton.x;/* rotate by x displacement*/
            mx=mx+sc.winwidth;
 	    break;
      case XK_s:
 	    mstate = 3;
	    mx = event->xbutton.x;
            mx=mx+sc.winwidth;
 	    break;
      case XK_d:
 	    mstate = 4;
	    mx = event->xbutton.x;
            mx=mx+sc.winwidth;
 	    break;
      case XK_q:
 	    mstate = 5;
	    mx = event->xbutton.x;/* rotate by x displacement*/
            mx=mx+sc.winwidth;
 	    break;
      case XK_w:
 	    mstate = 6;
	    mx = event->xbutton.x;
            mx=mx+sc.winwidth;
 	    break;
      case XK_e:  //rsize
 	    mstate = 7;
	    mx = event->xbutton.x;
            mx=mx+sc.winwidth;
 	    break;
      default: break;
      }
      break; 
     case ButtonPress:
        if (event->xbutton.button == Button3) {//If right button and no letter - rotate
	    if(mstate ==0){
               mstate = 1;
	       mx = event->xbutton.x;
               mx=mx+sc.winwidth;
	       my = event->xbutton.y;
            }
        }
        else if (event->xbutton.button == Button2) { // If push middle button and no letter pushed - translate
	    if(mstate ==0){
               mstate = 8;
	       mx = event->xbutton.x;
               mx=mx+sc.winwidth;
	       my = event->xbutton.y;
            }
	}
	break;
       
    case ButtonRelease:
        if( (event->xbutton.button == Button2)|| (event->xbutton.button == Button3) )
	    mstate = 0;
	break;
    case MotionNotify:
	if (mstate) {
	    omx = mx;
	    omy = my;
	    mx =event->xbutton.x;
            mx=mx+sc.winwidth;
	    my = event->xbutton.y;
	    update_view(mstate, omx,mx,omy,my);
	}
	break;
   }
}


void drawline(float fvec[3],float fvec2[3])
{
   glBegin(GL_LINE_STRIP);
     glVertex3fv(fvec);
     glVertex3fv(fvec2);
   glEnd();
}   

void drawx(float xlen,float pos[3])
// Draws X at world position pos, of lenght xlen (angstroms)
// Assume color has already been assigned
{
   glBegin(GL_LINE_STRIP);
     glVertex3f(pos[0]-xlen,pos[1],pos[2]);
     glVertex3f(pos[0]+xlen,pos[1],pos[2]);
   glEnd();
   glBegin(GL_LINE_STRIP);
     glVertex3f(pos[0],pos[1]-xlen,pos[2]);
     glVertex3f(pos[0],pos[1]+xlen,pos[2]);
   glEnd();
   glBegin(GL_LINE_STRIP);
     glVertex3f(pos[0],pos[1],pos[2]-xlen);
     glVertex3f(pos[0],pos[1],pos[2]+xlen);
   glEnd();
}


void drawpeptide(struct geometry& geo)  
/*
 = drawgeo 
New version: eacg side chain has individual # of atoms
Peptide (=geo)  begins at Ca for residue 1 and ends at Ca for residue geo.nafit
   Sidechains not drawn for last residue 
Define  1=ca,2=c,3=o;4=n,5=cb; 6=cb2 in geo.pos[i]...
*/
{
   int i,j,i1,count,atomname,rottype;

   glLineWidth(geo.linewidth);
   for(i=geo.nend[0]+1;i<=geo.nend[1];i++){ /* nafid is number of serines in fit */
//printf("geo.pos[%d-1][1,2,3,4][0]= %5.3f %5.3f %5.3f %5.3f\n",
//    i,geo.pos[i-1][1][0],geo.pos[i-1][2][0],geo.pos[i-1][3][0],geo.pos[i-1][4][0]);
     glPushName(i-geo.nend[0]); /*put residue # on stack*/
     glColor3fv(geo.geomaincol);  //This the simpler, vector version
     glBegin(GL_LINE_STRIP);
      glVertex3fv(geo.pos[i-1][1]);   
      glVertex3fv(geo.pos[i-1][2]);   
      glVertex3fv(geo.pos[i][4]);   
      glVertex3fv(geo.pos[i][1]); 
    glEnd(); 
    if( i<geo.nend[1] ){ /*Draw sideshain  second condition is test if glycine*/
       glColor3fv(geo.geosidecol);  //This the simpler, vector version
       if(geo.iside[i]>=0){ // indicates side is just series of C atoms
         glBegin(GL_LINE_STRIP);
         glVertex3fv(geo.pos[i][1]);
         for(j=1;j<=geo.iside[i];j++)  // draw side chain atoms  
           glVertex3fv(geo.pos[i][4+j]);   
         glEnd(); 
       } // >=0 
#if 1
       else{ //nside <0 indicates side is residue type = rottype
         rottype = -geo.iside[i];
         count = 0; 
         for(i1=1;i1<=side[rottype][0][0]+1;i1++){ /* draw each side chain branch*/
            glBegin(GL_LINE_STRIP);
            for(j=0;j<side[rottype][i1][0];j++){
              atomname=side[rottype][i1][j+1];
//printf("atomname = %d\n",atomname);
              if(atomname<=5)             
                 glVertex3fv(geo.pos[i][atomname]);
              else {
                 count = count+1;
                 glVertex3fv(geo.pos[i][count+5]);
              }
            }//j   
            glEnd();
         } //i1
      } // else 
#endif
    }
    glPopName();
  } /* i*/
  for(i=geo.nend[0];i<geo.nend[1];i++){ /*draw C-O bond*/
     glPushName(i-geo.nend[0]); /*put res # on stack*/
     glColor3fv(geo.geococol);  //This the simpler, vector version
     glBegin(GL_LINE_STRIP);
       glVertex3fv(geo.pos[i][2]);   
       glVertex3fv(geo.pos[i][3]);   
     glEnd();  
     glPopName();
  }
}


void drawfit(struct geometry& geo,struct pepfit fit[NFIT],
        int fitnum[NFIT])
/* 
Almost identical to drawpeptide - except draws all
 saved fits
Define  1=ca,2=c,3=o;4=n,5=cb; 6=cb2 in geo.pos[i]...
*/
{
 int i,j,ifit,ii,i1,rottype,atomname,count;

 glLineWidth(geo.linewidth);
 float maincolor[NFIT][3];

 if(geo.fuse[0]>0){ // then must change color of picked main chain
   for(ifit=1;ifit<=fitnum[0];ifit++){ // assign main color
     ii = fitnum[ifit]; // the bin fit is stored in
     if( (ii == geo.fuse[1])||(ii==geo.fuse[2]) ){
         for(j=0;j<=2;j++)
            maincolor[ii][j]= greenv[j];
     }
     else{ 
         for(j=0;j<=2;j++)
            maincolor[ii][j]= geo.fitmaincol[j];
     }
   } // ifit
 } // geo.fuse > 0
 else{//geo.fuse <= 0
   for(ifit=1;ifit<=fitnum[0];ifit++){ // assign main color
     ii = fitnum[ifit]; // the bin fit is stored in
     for(j=0;j<=2;j++)
        maincolor[ii][j]= geo.fitmaincol[j];
   }
 } //geo.fuse = 0
 for(ifit=1;ifit<=fitnum[0];ifit++){ // all stored fits
  ii = fitnum[ifit]; // the bin fit is stored in
#if 0
printf("drawing fit: ifit = %d ii=%d  caend[x] = %5.3f\n",
     ifit,ii, fit[ii].pos[fit[ii].nres-1][1][0]); 
printf("   ii=%d  geo.pickfitnum=%d  sc.drawgeo= %d\n",ii,geo.pickfitnum,sc.drawgeo);
#endif
  if( (ii!=geo.pickfitnum) || (sc.drawgeo== -1) ){
   glPushName(ii); /*put pdbfile # on stack*/
   for(i=1;i<fit[ii].nres;i++){ /* nafid is number of serines in fit */     
     glPushName(i); /*put res # on stack*/
//     glColor3fv(geo.fitmaincol);  //This the simpler, vector version
     glColor3fv(maincolor[ii]);  //This the simpler, vector version
     glBegin(GL_LINE_STRIP);
      glVertex3fv(fit[ii].pos[i-1][1]); 
//printf(" fit[%d].pos[%d-1][1][0] %5.3f\n",ii,i,fit[ii].pos[i-1][1][0]); 
      glVertex3fv(fit[ii].pos[i-1][2]);   
//printf(" fit[%d].pos[%d-1][2][0] %5.3f\n",ii,i,fit[ii].pos[i-1][2][0]); 
      glVertex3fv(fit[ii].pos[i][4]);   
//printf(" fit[%d].pos[%d][4][0] %5.3f\n",ii,i,fit[ii].pos[i][4][0]); 
      glVertex3fv(fit[ii].pos[i][1]); 
//printf(" fit[%d].pos[%d][1][0] %5.3f\n",ii,i,fit[ii].pos[i][1][0]); 
    glEnd();

    if( i<fit[ii].nres-1 ){ /*Draw sideshain  second condition is test if glycine*/
       glColor3fv(geo.fitsidecol);  //This the simpler, vector version
       if(fit[ii].nside[i]>=0){ // indicates side is just series of C atoms
//printf("fit[%d].nside[%d] = %d\n",ii,i,fit[ii].nside[i]);
           glBegin(GL_LINE_STRIP);
           glVertex3fv(fit[ii].pos[i][1]);
           for(j=1;j<=fit[ii].nside[i];j++){  // draw side chain atoms  
             glVertex3fv(fit[ii].pos[i][4+j]);
//printf(" fit[%d].pos[%d][4+%d][0] %5.3f\n",ii,i,j,fit[ii].pos[i][4+j][0]); 
           }   
           glEnd(); 
       }
#if 1
       else{ //nside <0 indicates side is residue type = rottype
         rottype = -fit[ii].nside[i];
         count = 0; 
         for(i1=1;i1<=side[rottype][0][0]+1;i1++){ /* draw each side chain branch*/
            glBegin(GL_LINE_STRIP);
            for(j=0;j<side[rottype][i1][0];j++){
              atomname=side[rottype][i1][j+1];
              if(atomname<=5)             
                 glVertex3fv(fit[ii].pos[i][atomname]);
              else {
                 count = count+1;
                 glVertex3fv(fit[ii].pos[i][count+5]);
              }
            }//j   
            glEnd();
         } //i1
      } // else 
#endif
    }
    glPopName();
   } /* i*/
   for(i=0;i<fit[ii].nres-1;i++){ /*draw C-O bond*/
     glPushName(i+1); /*put res # on stack*/
     glColor3fv(geo.fitcocol);  //This the simpler, vector version
     glBegin(GL_LINE_STRIP);
       glVertex3fv(fit[ii].pos[i][2]);   
       glVertex3fv(fit[ii].pos[i][3]);   
     glEnd();  
     glPopName();
   }
   glPopName();
 } // ii!= pickfitnum
} // ifit
}

void draweditbone(struct extbones ext)
/* Draw the edited bones used to extend fit
*/
{
  int i,j,k;

  for(i=0;i<=2;i++){
      ext.mecolor[i]=pinkv[i];
      ext.secolor[i]=lightbluev[i];
  }
   glLineWidth(geo.linewidth); // thin lines all the time
//glLineWidth(2.0); // thin lines all the time
   if(ext.nmain == -150){// this indicates have converted main to onemain - no sides
     glColor3fv(ext.mecolor);  //This the simpler, vector version    
     glPushName(0); // push the number of the mainone bone- used for drawing X and connecting line
     glBegin(GL_LINE_STRIP);
     for(j=0;j<ext.mainone.num;j++) // draw the main pts
          glVertex3fv(ext.mainone.mainpt[j].pos);
     glEnd();
     glPopName();
   }
   else for(i=1;i<=ext.nmain;i++){
     glColor3fv(ext.mecolor);  //This the simpler, vector version    
     glPushName(i); // push the number of the bonept- used for drawing X and connecting line
     glBegin(GL_LINE_STRIP);
     for(j=0;j<ext.main[i].num;j++) // draw the main pts
          glVertex3fv(ext.main[i].mainpt[j].pos);
     glEnd();
     glColor3fv(ext.secolor);  //This the simpler, vector version 
     for(j=1;j<=ext.main[i].nside;j++){ // each side for this main
        glBegin(GL_LINE_STRIP);
        for(k=0;k<ext.main[i].side[j].num;k++)   
           glVertex3fv(ext.main[i].side[j].pos[k]);
        glEnd();
     }
     glPopName();
  }//i
}     

void drawonebone(struct griddata& den1,struct bonedata& bone1,int i)
{
   int ivec[3],ivec2[3],j,icon;
   float xlen = 0.2;
   float fvec[3],fvec2[3];

     d1tod3(bone1.cmax,bone1.ptv[i].pos,ivec);
     imcv(den1,ivec,fvec);// convert int grid ivec to float world co-ord fvec
     glPushName(i); // push the number of the bonept- used for drawing X and connecting line
     if(sc.pickbone==1){
        if(i==bone1.pickbone){
           glColor3fv(whitev);  //This the simpler, vector version    
//           xlen = 0.5;
           xlen = 3;
           drawx(xlen,fvec); // draw x at bone pt pos
           glColor3fv(bone1.mcolor);  //This the simpler, vector version    
           xlen = 0.2;
        }
        else
          drawx(xlen,fvec); // draw x at bone pt pos
     }
     else
        drawx(xlen,fvec); // draw x at bone pt pos
     for(j=0;j<bone1.ptv[i].num;j++){ // draw lines to higher nubmer pts
//printf("num = %d cen =%d con = %d\n",(*i).num, (*i).pos,(*i).con[j]);
       icon = bone1.ptv[i].con[j]; // the pt num of the connect point
       if( icon> i){ // draw lines from lower number to higher numbered points
         d1tod3(bone1.cmax,bone1.ptv[icon].pos,ivec2);
          imcv(den1,ivec2,fvec2);// convert int grid ivec to float world co-ord fvec
          drawline(fvec,fvec2); 
        }
     }// j 
}       


void drawbone(struct griddata& den1,struct bonedata& bone1,int windowselect)
/* New version - where bone.ptv.con is the bone number  of connected point
   Bone window  (bone1.winonoff =1)
     winpts[1,2, .... max]  = # of bonept in window
   where max = winpts[0],  max < NBONEWIN
*/
{
  int i,ii;
   glColor3fv(bone1.mcolor);  //This the simpler, vector version    
//   glLineWidth(geo.linewidth);
   glLineWidth(1.0); // thin lines all the time
//glLineWidth(2.0); // thin lines all the time
   int totnum = bone1.ptv.end()-bone1.ptv.begin();
//printf("windowselect = %d\n",windowselect);
   if(windowselect == -1){ // window off
    for(i=0;i<totnum;++i){
       drawonebone(den1,bone1,i);
     glPopName();
    } // i 
   } // window off
   else{ // window on
    for(ii=1;ii<=bone1.winpts[0];++ii){
     i= bone1.winpts[ii];
     drawonebone(den1,bone1,i);
     glPopName();
    } // ii 
   } // window on
}

void definepdbcolor(int nseg,float **maincolor,float **sidecolor)
// this was in conv_lev, but removed so no screen stuff
{
  int i,k;
   for(k=0;k<=2;k++) { /* first seg is yellow and blue */
     maincolor[i][k]=yellowv[k];
     sidecolor[i][k]=lightbluev[k];
   }
 for(i=1;i<=nseg;i++){
   for(k=0;k<=2;k++) { /* all seg after first are in cyan and pink */
     maincolor[i][k]=cyanv[k];
     sidecolor[i][k]=pinkv[k];
   }
 }
}

#if 1  // has screen stuff - moved to drawwindow.h
void savepick(struct geometry& geo)
/* 
New - 26may99 - If fit result of fusion, remove higher number of fused fits and renumber
New - 10may99 - First read fitfile into fittemp - then
   write out file  -substituting new picked fit
Similar to save fit - except there is a picked fit (fit #=geo.pickfitnum)
  that hasbeen converted to geo and has been altered
Saves all the fit files - and the one geo file in same number order
*/
{
   int i,numfits,irank0,ifuse;
   FILE *fp;
   int nii[2];
   struct pepfit fittemp[NFIT];
   char backup[50]; 

   strcpy(backup,maidfile.fitfile);
   strcat(backup,"_backup");
   ifuse = 0;
#if 0
int fitsavnum;
   if( (fit[geo.pickfitnum ].rankfit[0]>0) ) // if picked fit is part of rankset - save in place of orginal
        fitsavenum = fit[geo.pickfitnum].rankfit[1];
   else
       fitsavenum = geo.pickfitnum;
#endif
   printf("Saving fitfile %s and copy to backup %s    geo.fuse[1]= %d\n",
               maidfile.fitfile,backup,geo.fuse[1]);
   printf("    Replacing original fit %d with picked fit %d\n",
                          fit[ geo.pickfitnum ].rankfit[0],geo.pickfitnum);
 for(i=1;i<NFIT;i++) //Initial default allocation - always allocate because new each time routine called
     fittemp[i].nalloc = -10;   
 if( (fp=fopen(maidfile.fitfile,"r")) == NULL){
     printf("ERROR: There are no fit files\n" );
 }
 else{ //fitfile opened
   fscanf(fp,"%d",&numfits); // total  # of fits;
   for(i=1;i<=numfits;i++){ // each fit
     irank0=0;  // there are no rankings for files read in readfit
     readonefit(fp,fittemp,irank0,i);
   }// each fit 
   fclose(fp);  
 } // fitfile opend
 fp=fopen(backup,"w"); // Save original fitfile in backup
 fprintf(fp,"%d\n",numfits); // total  # of fits
 for(i=1;i<=numfits;i++){
       saveonefit(fp,fittemp,i,i);
 }
 fclose(fp);
 fp= fopen(maidfile.fitfile,"w");
 if(geo.fuse[1]>=0){ // then this fit is result of fusing of two fits
    ifuse = 1;
    if(geo.fuse[2]<geo.fuse[1]){ // order fits so fuse[1] is smaller 
      i= geo.fuse[1];
      geo.fuse[1]=geo.fuse[2];
      geo.fuse[2]= i;
    }
 }
 if(ifuse == 0) { // Not fused fits
   fprintf(fp,"%d\n",numfits); // total  # of fits
   for(i=1;i<=numfits;i++){
     if(i!=fit[ geo.pickfitnum ].rankfit[0]){// not the fit that is being extended
       saveonefit(fp,fittemp,i,i);
     }// not the picked fit
     else{ // the picked fit = geo 
        nii[0]=fittemp[i].val1; // use original values
        nii[1]=fittemp[i].val2;
        copygeotofit(geo,fittemp[i],nii);
        saveonefit(fp,fittemp,i,i);
     }//  the picked fit
   } // i each fit
   fclose(fp); 
 } // not fused
 else {  // fused fits
   fprintf(fp,"%d\n",numfits-1); // remove one fit -  total  # of fits
   for(i=1;i<geo.fuse[1];i++) // write out fits up to fuse[1]
      saveonefit(fp,fittemp,i,i);
   nii[0]=fittemp[i].val1; // use original values
   nii[1]=fittemp[i].val2;
   i=geo.fuse[1];
   copygeotofit(geo,fittemp[i],nii);
   saveonefit(fp,fittemp,i,i); // saved fuse fit in place of lower number fit
   for(i=geo.fuse[1]+1;i<geo.fuse[2];i++) // save intervening fits with original number 
      saveonefit(fp,fittemp,i,i);
   for(i=geo.fuse[2]+1;i<=numfits;i++) // save remaining fits with number reduced by 1
      saveonefit(fp,fittemp,i,i-1);
   geo.fuse[0] = geo.fuse[1] = geo.fuse[2]= -1;  // reset to defaault vale
   fclose(fp); 
 }// fused fits  
 endsave:; 
}
#endif
