#include "drawwindow.h"  // the data that needs to be seen by other routines
#include "screen.h"
#include "pdbtolev.h"
#include "model.h"
#include "peptide.h"
#include "math.h"
#define EXTERN extern  /*when read def.h make all EXTERN variables extern*/
#include "global.h"  // global variables  color functions  (external variable)
extern struct maiddatafile maidfile;// Need global variable for X-menu calls


// Original old main variables - now allocated at run time

void findmodelphipsi(int ires,struct pdbfile& pdb1,int is,
       float *phi,float *psi)
/*
Modified from routine findphipsi using pos[][][] as input
Find the phi/psi angle for residue i - 
for i: phi is angle formed by Ci-Ni+1-Cai+1-Ci+1
       psi is angle formed by Ni+1-Cai+1-Ci+1-Ni+2
Returns the phi/psi for residue ires,(segment = is) which means
   i in above definition = ires-1
 mainch[0,1,2] = N,Ca,C
*/
{
   findchi(pdb1.mainch[is][ires-1][2],pdb1.mainch[is][ires][0],
       pdb1.mainch[is][ires][1],pdb1.mainch[is][ires][2],phi);   
   findchi(pdb1.mainch[is][ires][0],pdb1.mainch[is][ires][1],
       pdb1.mainch[is][ires][2],pdb1.mainch[is][ires+1][0],psi);   
}  

void drawrot(int restype,float rotpos[NATOMS][3],int side[AA][NBRN][NPTS])
/* 
  Draws the  rotomer for amino acid type restype, using
    the positions determined by buildrot routine
 */
{
    int ii,j,k,atomname;
    float vertx[3];

//printf("restype = %d\n",restype);
    glColor3fv(yellowv);  //This the simpler, vector version
    for(ii=1;ii<=side[restype][0][0]+1;ii++){ /* draw each side chain branch*/
       glBegin(GL_LINE_STRIP);
           for(j=0;j<side[restype][ii][0];j++){
              atomname=side[restype][ii][j+1];
              for(k=0;k<=2;k++)
                  vertx[k]=rotpos[atomname][k];
//printf("atomname = %d  branch = %d  j= %d  vertx = %5.3f %5.3f %5.3f\n",
//    atomname,ii,j,vertx[0],vertx[1],vertx[2]);
              glVertex3fv(vertx);
           }   
       glEnd();
    } //ii
}

void drawresidue(struct pdbfile& pdb1,int is,int i,
       int iprev,int windowselect)
/* draws one residue for segment is residue i
 mainch[0,1,2] = N,Ca,C
 sideonoff = +/- 1 (connect all mainchain atoms, + draw side)
           = +/- 2  Connect only Ca's
*/
{
   int ii,j,ig;
   
    glColor3fv(pdb1.maincolor[is]);  //This the simpler, vector version
   glBegin(GL_LINE_STRIP);
   for(ig=1;ig<=pdb1.gaps[is][0];ig++){
            if(pdb1.gaps[is][ig]==i){ /* implies gap, do not connect */
//printf("i=%d  agap at ig = %d\n",i,ig);
                goto agap;
            }
   }
   if( windowselect==1 ){
     if(i==iprev+1){ /* do not connect if not conseecuative*/
       if(abs(pdb1.sideonoff)==1) /* connect all mainch residues*/       
                glVertex3f(pdb1.mainch[is][i-1][2][0],pdb1.mainch[is][i-1][2][1],pdb1.mainch[is][i-1][2][2]);   
       else  // sideonoff = 0  indicates Ca trace
                glVertex3f(pdb1.mainch[is][i-1][1][0],pdb1.mainch[is][i-1][1][1],pdb1.mainch[is][i-1][1][2]);   
     }
   }
   else{ /* if windows not selected, always connect */
        if(abs(pdb1.sideonoff)==1){ /* connect all mainch residues*/ 
                glVertex3f(pdb1.mainch[is][i-1][2][0],pdb1.mainch[is][i-1][2][1],pdb1.mainch[is][i-1][2][2]);  
       } 
       else /*sideonfff =+/- 2 - indicates  Connect just the Ca's */
                glVertex3f(pdb1.mainch[is][i-1][1][0],pdb1.mainch[is][i-1][1][1],pdb1.mainch[is][i-1][1][2]);   
   }
   agap:;
   if(abs(pdb1.sideonoff)==1){ /* connect all mainch residues*/
      for(j=0;j<=2;j++) { /* complete redidue main chain (old program, toppt=2)*/
                glVertex3f(pdb1.mainch[is][i][j][0],pdb1.mainch[is][i][j][1],pdb1.mainch[is][i][j][2]);   
      }
   }
   else{ /* Connect just the Ca's */
                glVertex3f(pdb1.mainch[is][i][1][0],pdb1.mainch[is][i][1][1],pdb1.mainch[is][i][1][2]);   
   }
   glEnd();
   if(pdb1.sideonoff>0) { /*draw side chain*/
         for(ii=1;ii<=pdb1.resdat[is+1][i][0];ii++) {
            if(ii==pdb1.resdat[is+1][i][0])/* C-O bond*/
                 glColor3f(redv[0],redv[1],redv[2]);
            else 
              glColor3f(pdb1.sidecolor[is][0],pdb1.sidecolor[is][1],pdb1.sidecolor[is][2]);
            glBegin(GL_LINE_STRIP);
            for(j=0;j<pdb1.resdat[is+1][i][ii];j++) {
                glVertex3f(pdb1.sidech[is][i][ii][j][0],pdb1.sidech[is][i][ii][j][1],pdb1.sidech[is][i][ii][j][2]);   
            }
            glEnd();
         }
   } /*if sideonoff*/
} 

void allocatelev(struct pdbfile& pdb1) 
{  
   int iseg,i,k,maxres,nseg,nres[1000],nsubst,nwater; // Assume nseg < 1000
      findpdbdimensions(pdb1.name,&nseg,nres,&nsubst,&nwater);
      maxres = -1;
      for(i=0;i<nseg;i++)
        if(nres[i]>maxres-2)
          maxres= nres[i]+2;
printf("nseg = %d  nres[0] = %d nres[1] = %d maxres=%d nsubst= %d  nwater = %d\n",
               nseg,nres[0],nres[1],maxres,nsubst,nwater);
      imat2(pdb1.windmod,nseg,NRESWIND);
      imat2(pdb1.gaps,nseg,NGAPS);// Allocate memory
      imat3(pdb1.resdat,nseg+1,maxres,NBRN+1);
      fmat2(pdb1.maincolor,nseg,3);
      fmat2(pdb1.sidecolor,nseg,3);
      fmat4(pdb1.mainch,nseg,maxres,3,3);
      fmat5(pdb1.sidech,nseg,maxres,NBRN+1,NPTS,3);
      cmat3(pdb1.rname,nseg,maxres,4);
      cmat3(pdb1.resnumch,nseg,maxres,5);
//      if(nsubst>0){ // always need to allocate
          fmat4(pdb1.subch,nsubst+1,SUBBRN,SUBNPTS,3);
          imat2(pdb1.subdata,nsubst+1,SUBBRN);
//      }
      pdbtolev(nseg,nwater,maxres,nsubst,pdb1.subch,pdb1.rname,pdb1.subdata,pdb1.resnumch,
           pdb1.mainch,pdb1.sidech,pdb1.gaps,pdb1.resdat,
           pdb1.maincolor,pdb1.sidecolor,pdb1.name,pdb1.xmin,pdb1.xmax);
      iseg = 0;
      for(k=0;k<=2;k++) { /* first seg is yellow and blue */
         pdb1.maincolor[iseg][k]=pdb1maincol[k];
         pdb1.sidecolor[iseg][k]=pdb1sidecol[k];
      }
      if(iseg == nseg-1)
         goto skipcolor;
      iseg = 1;
      for(k=0;k<=2;k++) { 
         pdb1.maincolor[iseg][k]=cyanv[k];
         pdb1.sidecolor[iseg][k]=pinkv[k];
      }
      if(iseg == nseg-1)
         goto skipcolor;
      iseg = 2;
      for(k=0;k<=2;k++) {
         pdb1.maincolor[iseg][k]=bluev[k];
         pdb1.sidecolor[iseg][k]=redv[k];
      }
      if(iseg == nseg-1)
         goto skipcolor;
      iseg = 3;
      for(k=0;k<=2;k++) { 
         pdb1.maincolor[iseg][k]=redv[k];
         pdb1.sidecolor[iseg][k]=yellowv[k];
      }
      if(iseg == nseg-1)
         goto skipcolor;
      iseg = 4;
      for(k=0;k<=2;k++) { 
         pdb1.maincolor[iseg][k]=lightbluev[k];
         pdb1.sidecolor[iseg][k]=pinkv[k];
      }
      if(iseg == nseg-1)
         goto skipcolor;
      iseg = 5;
      for(k=0;k<=2;k++) { 
         pdb1.maincolor[iseg][k]=pinkv[k];
         pdb1.sidecolor[iseg][k]=redv[k];
      }
      if(iseg == nseg-1)
         goto skipcolor;
      iseg = 6;
      for(k=0;k<=2;k++) { 
         pdb1.maincolor[iseg][k]=magentav[k];
         pdb1.sidecolor[iseg][k]=yellowv[k];
      }
      if(iseg == nseg-1)
         goto skipcolor;
      iseg = 7;
      for(k=0;k<=2;k++) { 
          pdb1.maincolor[iseg][k]=greenv[k];
          pdb1.sidecolor[iseg][k]=lightbluev[k];
      }
      if(iseg == nseg-1)
         goto skipcolor;
      iseg = 8;
      for(k=0;k<=2;k++) { 
          pdb1.maincolor[iseg][k]=mycolor1[k];
          pdb1.sidecolor[iseg][k]=lightbluev[k];
      }
      if(iseg == nseg-1)
         goto skipcolor;
      iseg = 9;
      for(k=0;k<=2;k++) { 
          pdb1.maincolor[iseg][k]=mycolor2[k];
          pdb1.sidecolor[iseg][k]=lightbluev[k];
      }
      if(iseg == nseg-1)
         goto skipcolor;
      iseg = 10;
      for(k=0;k<=2;k++) { 
          pdb1.maincolor[iseg][k]=mycolor3[k];
          pdb1.sidecolor[iseg][k]=lightbluev[k];
      }
      if(iseg == nseg-1)
         goto skipcolor;
      iseg = 11;
      for(k=0;k<=2;k++) { 
          pdb1.maincolor[iseg][k]=mycolor4[k];
          pdb1.sidecolor[iseg][k]=lightbluev[k];
      }
      if(iseg == nseg-1)
         goto skipcolor;
      iseg = 12;
      for(k=0;k<=2;k++) { 
          pdb1.maincolor[iseg][k]=mycolor5[k];
          pdb1.sidecolor[iseg][k]=lightbluev[k];
      }
      if(iseg == nseg-1)
         goto skipcolor;
      iseg = 13;
      for(k=0;k<=2;k++) { 
          pdb1.maincolor[iseg][k]=mycolor6[k];
          pdb1.sidecolor[iseg][k]=lightbluev[k];
      }
      for(i=14;i<nseg;i++){ // For remaining segments - just rotate through first 1-7 colors
        iseg = i % 14;
        for(k=0;k<=2;k++) { 
          pdb1.maincolor[i][k]=pdb1.maincolor[iseg][k];
          pdb1.sidecolor[i][k]=pdb1.sidecolor[iseg][k];
        }
      }
   skipcolor:;
   for(i=0;i<nseg;i++){ 
        for(k=1;k<=nres[i];k++)  // This resdata[][]4] is used as indicator of whether residue has been changed and needs to be rewritten in writelevtopdb
            pdb1.resdat[i+1][k][4] = 0; //0 - default - not changed FIXED - increase NBRN to 5 and change 100 to 4 WARNING - OVERFLOW - NOT USED ANYMORE indicates has been modidied
   }
   quitallocate:;
}

void drawmodel(struct geometry& geo,struct screenvalue&  sc,struct maiddatafile& maidfile,
      int& windowselect)
/* 
New C++ version - Assume want to draw all segments in a pdbfile
   unless have picked pdb - they just draw that segment.
As default - draw only first pdb in file -
   Use this pdb to center and scale view
mainch[0,1,2] = N,Ca,C
NOTE:  if drawpickdbb =1, then draw entire picked pdb, not the 
   windowed pdb
New feature:  draw only the segments listed in segdraw
  segdraw[0]= # of segments to draw, segdraw[i]= segment #
*/
{  
 int is,i,ip,iprev,ires;

   for(ip=1;ip<=maidfile.npdb;ip++){ // each pdb file
     if(windowselect ==1)
          glLineWidth(2.0);
     else{ 
          glLineWidth(geo.linewidth);
     }
     glPushName(ip); /*put pdbfile # on stack*/
     if(maidfile.pdb[ip].onoff==1){ // if turned on
        if(maidfile.pdb[ip].start==1){
printf(" pdb[%d].start = %d maidfile.start = %d\n",ip,maidfile.pdb[ip].start,maidfile.start);
           allocatelev(maidfile.pdb[ip]);
#if 0 // find angle between N - Ca - C for all residues
 int ii,k,count;
float ang,dist1,dist2,cav[3],cgv[3],ave;
is = 0; // segment number - default assume only 1 segment
ip = 1; // peptide number 
  count =0;
  ave = 0.0;
  for(i=2;i<maidfile.pdb[ip].resdat[is+1][0][0];i++){ // i = each residue  
//     if( strcmp(maidfile.pdb[ip].rname[is][i],"TRP")==0 ){
         for(k=0;k<=2;k++){
            cav[k] =maidfile.pdb[ip].mainch[is][i][1][k]- maidfile.pdb[ip].mainch[is][i][0][k]; //Vector from Ca to N  
            cgv[k] =maidfile.pdb[ip].mainch[is][i][1][k]- maidfile.pdb[ip].mainch[is][i][2][k]; //Vector from Ca to C  
         }
         dist1 = sqrt( dot(cav,cav) );
         dist2 = sqrt( dot(cgv,cgv) );
         ang = acos( dot(cav,cgv)/(dist1*dist2) );
         ang = 180.0*ang/PI;
         count = count+1;
         ave = ave+ang;
         printf(" res # %d: angle = %5.3f\n",i,ang);
//     }
   }// i 
   ave = ave/count;
   printf(" Ave ang = %5.3f\n",ave);
#endif
#if 0 // find angle between Cb-Ca and Cb-Cg for ring compounds
 int ii,k,count;
float ang,dist1,dist2,cav[3],cgv[3],ave;
is = 0; // segment number - default assume only 1 segment
ip = 1; // peptide number 
  count =0;
  ave = 0.0;
  for(i=2;i<maidfile.pdb[ip].resdat[is+1][0][0];i++){ // i = each residue  
     if( strcmp(maidfile.pdb[ip].rname[is][i],"TYR")==0 ){
         for(k=0;k<=2;k++){
            cav[k] =maidfile.pdb[ip].sidech[is][i][1][0][k]- maidfile.pdb[ip].sidech[is][i][1][1][k]; //Vector from Cb to Ca   
            cgv[k] =maidfile.pdb[ip].sidech[is][i][1][2][k]- maidfile.pdb[ip].sidech[is][i][1][1][k]; //Vector from Cg to Ca  
         }
         dist1 = sqrt( dot(cav,cav) );
         dist2 = sqrt( dot(cgv,cgv) );
         ang = acos( dot(cav,cgv)/(dist1*dist2) );
         ang = 180.0*ang/PI;
         count = count+1;
         ave = ave+ang;
         printf(" For TYR, res # %d: angle = %5.3f\n",i,ang);
     }
   }// i 
   ave = ave/count;
   printf(" Ave ang = %5.3f\n",ave);
#endif
#if 0 // print out phi/psi angles
/* classifyphipsi returns: 
     nclass = number of constraint types for which angle falls in geo.sheet range.
       if nclass > 0; type[1,2,.. nclass] are the types.
        IF nclass = 0 , type[1,2] = type of closest and next closest region
        If nclass = 1;  type[1] = type for nclass, type[2] = next closest radiee
*/
float phi,psi,radius[10];
int type[10],nclass;
is = 0; // segment number - default assume only 1 segment
ip = 1; // peptide number 
printf(" This is a routine in model.c++  (drawmodel)that is normal commenet out \n");
printf("     See  // print out phi/psi angles\n");
 for(i=2;i<maidfile.pdb[ip].resdat[is+1][0][0];i++){ // i = each residue  
   findmodelphipsi(i,maidfile.pdb[ip],is,&phi,&psi);
   classifyphipsi(geo,phi,psi,&nclass,type,radius); // Uses geo.sheet parameters
   if(strcmp(maidfile.pdb[ip].rname[is][i],"PRO")==0 )//Print only if Proline
     printf("Res. %d %s: phi = %6.3f  psi = %6.3f  nclass = %d  type1 = %d type 2 = %d\n",
        i,maidfile.pdb[ip].rname[is][i],phi*180.0/PI,psi*180.0/PI,nclass,type[1],type[2]);
 }
#endif
#if 0 // find N[i]-O[i-2] and O[i]-O[i-1] dist
float dist,dist2,mino_n,mino_o;
int nbrn,nbrn2,imin,imin2;
is = 0; // segment number - default assume only 1 segment
ip = 1; // peptide number 
printf(" This is a routine in model.c++  (drawmodel)that is normal commenet out \n");
#if 0 //for debugging
i=1;
     nbrn= maidfile.pdb[ip].resdat[is+1][i][0]-1; // # of branches in residue
printf("nbrn = %d O[1][0] = %5.3f  n[1][0] = %5.3f\n",
        nbrn,maidfile.pdb[ip].sidech[is][i][nbrn+1][1][0], maidfile.pdb[ip].mainch[is][i][0][0]);
#endif
 mino_o = mino_n = 1000.0;
 for(i=3;i<maidfile.pdb[ip].resdat[is+1][0][0];i++){ // i = each residue  
     nbrn= maidfile.pdb[ip].resdat[is+1][i-2][0]-1; // # of branches in residue
     dist=distv(maidfile.pdb[ip].sidech[is][i-2][nbrn+1][1], maidfile.pdb[ip].mainch[is][i][0]); // O[i-2} to N[i] dist
     nbrn= maidfile.pdb[ip].resdat[is+1][i-1][0]-1; // # of branches in residue
     nbrn2= maidfile.pdb[ip].resdat[is+1][i][0]-1; // # of branches in residue
     dist2=distv(maidfile.pdb[ip].sidech[is][i-1][nbrn+1][1],maidfile.pdb[ip].sidech[is][i][nbrn2+1][1]);// O[i-1] to )[i] dist
//     printf("nbrn = %d   O[i-2] - N[i] dist = %5.3f   O[i-1] - O[i] dist = %5.3f\n",nbrn,dist,dist2);
     if(dist < mino_n){
        mino_n = dist;
        imin = i;
     }
     if(dist2 < mino_o){
        mino_o = dist2;
        imin2 = i;
     }
 }
 printf(" min o_n dist = %5.3f (i=%d)  mino_o dist = %5.3f (i=%d)\n",
       mino_n,imin,mino_o,imin2);
#endif

           maidfile.pdb[ip].start=0; // do not call allocatelev again for this pdb
           if(maidfile.start>=0){ // First time - Center about first pdb read 
               sc.centerpos[0]=(maidfile.pdb[ip].xmin[0]+maidfile.pdb[ip].xmax[0])/2.0; // center about first pdb
               sc.centerpos[1]=(maidfile.pdb[ip].xmin[1]+maidfile.pdb[ip].xmax[1])/2.0;
               sc.centerpos[2]=(maidfile.pdb[ip].xmin[2]+maidfile.pdb[ip].xmax[2])/2.0;
               sc.cmmaxh=2.0*(maidfile.pdb[ip].xmax[1]-maidfile.pdb[ip].xmin[1]);//  Maximum hight of visualized object
               sc.degfovy = sc.mag*(180.0/PI)*2.0*asin(sc.cmmaxh/(2.0*sc.cmeyedist));
               sc.run_wpid = XtAppAddWorkProc(sc.appctx,draw_screen,&sc.drawwind1);//This draw screen one time
               for(i=0;i<=2;i++){ // set rotation and translation back to 0
                  sc.rot[i] = sc.trans[i]=0;
               }
               maidfile.start = -1;
           }
        }
        for(is=0;is<maidfile.pdb[ip].resdat[0][0][0];is++){ // all segments
          glPushName(is); /*put segment # on stack*/
          if(windowselect !=1){
             for(i=1;i<=maidfile.pdb[ip].resdat[is+1][0][0];i++) { 
                glPushName(i); /*put res # on stack*/
                drawresidue(maidfile.pdb[ip],is,i,iprev,sc.windowselect);
                glPopName();
             }// each residue
          }
          else if(windowselect ==1){
             iprev= -3; /*value so do not connect first residue*/
             for(i=1;i<=maidfile.pdb[ip].windmod[is][0];i++){ // all residues in window
                ires = maidfile.pdb[ip].windmod[is][i];
                glPushName(ires); 
                drawresidue(maidfile.pdb[ip],is,ires,iprev,sc.windowselect);
                iprev=ires;
                glPopName();
             }
          }
          glPopName(); // is (segment)
        } // each segment
      } // is onoff = 1
      glPopName();// ip (pdb)
  }// each pdb file
}

void labelres() 
/* New version - prints the real residue number */
{
   char aname[10];
   float xpos[3];
   int segn,resn;
   glColor3fv(writecolv);
   for(int ip=1;ip<=maidfile.npdb;ip++){ // each pdb file
      for(int ii=1;ii<=maidfile.pdb[ip].pickres[0].num;ii++){ //loop through selected residues for this pdb
          segn =maidfile.pdb[ip].pickres[ii].seg;
          resn = maidfile.pdb[ip].pickres[ii].num; 
          for(int i=0;i<=2;i++) /*find position of CA  (="1")*/
              xpos[i]=maidfile.pdb[ip].mainch[segn][resn][1][i]; 
           glRasterPos3f(xpos[0],xpos[1],xpos[2]);
           strcpy(aname,maidfile.pdb[ip].rname[segn][resn]);
           strcat(aname,maidfile.pdb[ip].resnumch[segn][resn]);
    printString(aname);
    glFlush ();
      } // ii
   }// ip  
}

