#define EXTERN extern  /*when read def.h make all EXTERN variables extern*/
#include "global.h"  // global variables  color functions  (external variable)
//#include "math.h"

#define NANG 12  // Number of different angle types  1 (helix)  2 (sheet)
//Type 1 = helix, 2 = sheet, 3  = left helix 4= epsilon 5 = delta 6= gammar 7 = gammal 

#if 0 // old version - sets all arrays
struct pepfit
{
   int nres,nalloc;
   float ***pos;
   int *nside;
   int val1,val2,itype[NAFIT]; // = pt.val[1],pt.val[ni]
   int irank[NAFIT]; // irank[0] = # of previous rankings, irank[i] = ith ranking
   float den[NAFIT],qual; // density for this fit
   float fbad; // # of bad for this fit (float so that can be avebad)
   float fside; // ave # of side chains at c or n end of fit
   int cbden; //Value of cbden for last residue added to fit If = 1 then include Cb in den
};
#endif
#if 1 // old version -  arrays allocated by allocatefit
struct pepfit
{
   int nres;
   int nalloc; // initialized to -10, indicate not yet allocated
   int main[3]; // the main bone and point that fit is being extended on (NOTE - not saved to fit file)
   float ***pos;
   int *nside;
   int *irot;
   int val1,val2,*itype; // = pt.val[1],pt.val[ni]
   int *cis; // = 0 if standard trans proline, = 1 for cis proline
   int *irank; // irank[0] = # of previous rankings, irank[i] = ith ranking
   float *den,qual; // density for this fit
   float fbad; // # of bad for this fit (float so that can be avebad)
   float fside; // ave # of side chains at c or n end of fit
   int cbden; //Value of cbden for last residue added to fit If = 1 then include Cb in den
   int rankfit[2]; //if part of rankfit then rankfit[0] = fit# it belongs to, rankfit[1] is it order in that group
/*NOTE: MUST BE CAREFUL WITH LAST 2 - THEY ARE NOT ALLOCATED IN 
     allocatefit or freed in freepepfit or copied but are handled separately
     To be safe should allocate and free them, each time used.
  It is allocated by call to readresqual
    Be careful not to call this if already allocated
*/
   float **qualaa,**quallast,**qualfirst; // qual for each of 20 a.a. for each res.
   int  **ibad; // ibad for each of 20 a.a. for each res.
   int itemp[2]; // int variable used within routine and then can be discarded.
   float ***ftemp; // float variable for each residue and main atom, used and then discarded
   int ***itemp2; // int variable for each atom 
};
#endif

#if 1
struct pepset
/* A pepset is a set of rankedfitfiles plus the individual best fit files:
set.fit[0] = number of fit files (usually numbered 1,2,3...)
set.fit[i] = rankedfit filename corresponding to ith fit file
set.rank[0][0] = number of fits that have associated ranked fits
set.rank[i][0] = number of ranked fits in rankedfitset  i
set.rank[i][j] = fit number corresponding to jth ranked fit in rankset i
Exampe:  suppose read only two fits - fit 3 and fit 5 - asign to fit1 and fit2
  then set.fit[1] = 3;  set.fit[2] = 5
Save rankedfits for fit 3 in file rankfit3
*/
{
   int fit[NFIT];
   int rank[NFIT][NSAVEFIT];
};
#endif

struct savegeo // save geo.pos and side and type data
{
   float spos[NAFIT][NAA][3];
   int sside[NAFIT],stype[NAFIT];
};

struct geometry // Parameters and values used for making peptides
{
  int al[2],da[2],de[3][3],ep[3][3][3]; // constants used for rotation, interpolation
#if 1 // These constants are needed for totomer routine
  // following values are constant angle and bond lenght values
       // from M. Levitt (1973)  - assigned in assignlevconst
  float cnc,anc,bnb,ccn,acn,acc,bcc,ccs,ccv,oao,nao,nan,cao,can,aaa,aab,aav;
  float aag,nbn,cbn,abc,aba,bbn,bbc,gbc,bbg,csc,css,agn,agb,ggn,agg,bgg,bng,ccc;
  float bbb;//five membered ring aaa is for six membered ring
  float cn,an,en,bn,gn,ao,ae,cv,av,cc,ac,ec,bc,aa,ab,bb,ag,bg,gg,cs,ss,eo;
  float tyr,phe; // New - 8nov00 - Ca-Cb-Cg angles for tyrosine and phe = his.
  float rotcol[11][3];
#endif
  // following values are standard distances etc,  assigned in findconstants()
  float tt[8],coneang1,coneang2,coneang4,dda[15],dsqrt[15],hbonddist;
  float cadist,cardist; // cadist is distance from Ca to Ca, rdist is sum of dist along atoms
  // following values are positions of mainch atoms in peptine residues
  float cw[NAFIT][3],nw[NAFIT][3],caw[NAFIT][3],ow[NAFIT][3], 
             cbw[NAFIT][3]; 
  int nend[2],nafit,begend[2]; // nend is begin and end of geo, begend is ends of tor dyn region, dir = 0 fits in for N to C directon 
  float fitmaincol[3],fitsidecol[3],fitcocol[3];
  float geomaincol[3],geosidecol[3],geococol[3];
  float pos[NAFIT][NAA][3],stpos[3];
  float phihel,psihel;
  float phi[NANG],psi[NANG],phi0[NANG],phi1[NANG],psi0[NANG],psi1[NANG];
  float phipsi[NAFIT][2];  // phi (=phipsi[0]) and psi(=phipsi[1]) for each residue
  float plerr,frden,extendden,sideden,arange,thtot,phitot,trandel,rotdel,axdel;
  float minrigid,poorrigid,angweight,bondwt;
  int tstep,itype,nbadwhole,nbadend,nstartres,axstep;
  int rtype[NAFIT]; // the itype for each indiviual residue
  int cis[NAFIT]; // = 0 for standard trans pro, = 1 for cis proline
  int maxside;  // maximum # of side chain atoms
  int nside,iside[NAFIT];  // if side atoms if simple (iside) for each residue 
  int cbden[NAFIT]; // If = 1 then include Cb in den
  float stdist[3],curdist[3],errcurv[3],curlimdist[3]; // Distances along mainchain (1 = Helix, 2= Sheet)
  float covechel[3]; // Co-ordinates that relate Ca[i] to Ca[i+1] vector to average Co for helix
  int ncbang; // # of rotation steps to find best cb2
  int dir,endmove;// direction of dynamics, and number of end residues allowed to move
  int fixend; // 0 (move all), 1( fixed C term) 2( fixed N temp)
  float nbadwt; //subtract nbadwt*nbad/residue from avedensity 
  int savlen;// After extension, total length must be this long for fit to be saved
  int quit; // if point out of range in findforden, set quit = 1 and quit this trace
  float sheet[NCONSTRAINTS][8]; // parameter used to constrain phi/psi angles in tordyn 
  int setrad; // number of grid points about atom for set[][][] = -1
  int spherelim; //if = 1 In trace.c++, only use start points that are in spheres in spherfile
  int  usecb;  //if = 1 then use Cb for aveden and nbad, if = 2 use for just aveden, if = 0 do not use
  float minavecb2,minavecb;  // Minimum allowed value of ave Cb2 density to accempt final fit
  float minfrbad;  // Minimum allowed fraction of residues allowed to have bad atoms  
  float ca1_4min;// Minimum Ca[i] to Ca[i+4] dist for helix
  int ***set,nowset; //  keeps track of whether den pt has been already fit, nowset = 1 if set allocated
  int ntypes; // The number of constraint types that are tested
  int newfit; // if = 0, then start tracepep new; if =1 - continue on from fits in fitfile
  int oneside; // if = 0, then do standard dynamics, else just one side chain = oneside};
  int phenylden; // if = 1 then find density only over the last rigid ring of a.a. reisdue - e.ge. phenylaline
  int resside[NAFIT];// Top number of side chain atom - USED ONLY IF SIDE = RESIDUE
  float maxenddist; // maximum distance Ca ends allowed to move during "bestden" rfinemnt
  int fittype;  // =1 (helix) or 2(sheet) - type fit in checkbones
  float frpoor; // Allowed fraction of poor atoms in fit
  int obadn; // Number of bad O atoms in fit
  int res0[NAFIT]; // Number of atoms in each residue that have < 0 density
  int nisheet,nihel; // Cut off values of ni for sheet and helix
  int torend;  // Number of residues whose density is averaged in avectermden
  float linewidth;
  float aveatom[10];  // ave densith for atoms 1 to 6
  int npoormax; // Maximum allowed number of poor atom
  char setfile[10];  //Name of set describing set of ranked fits
  float fpos[2][7][3]; // pos of end residue - used when fuse two fits or run dyn. in center
  int fnum; // residue in geo that middle dynamics is run about 
  int userank,nmaxext; // in extend - which fit to extend, and how many to add
  int pickfitnum; // used to be sc.pickfitnum
  int fuse[3];  //fuse[0] = -1 (default); 0,1 or 2 (number of fits picked); fuse[1,2] = nterm, cterm fitnumber picked;  
  int setval,setside,setmain;  // val in set[][][] = setval indicated point is set
  int fitn; // A fit number that can be passed via geo
  float camaindist; // Max allowed distance of Ca from main bone
  float finalden,avecb2,avecb;
  int nbad,npoor;
  int main[3]; // the main bone that geo is being extended on (NOTE - not saved to fit file)
  int seqnum; // Number of residues in sequence file
  int *seq,maxseq,resnum; // geo.seq[i] is residue name of i residue, resnum = initial assigment
  int removefit;
  int badfit; // = 0 if  fit is too long to copytogeo ; = 1 if OK
  int extendfitnum; // the number of the fit that is currently being extended
  float maxfusedist,goodfusedist; // OK and good dist between end fits to be fused
  int checkmaindist; // if = 1, check distance from edited main in tordynbat
  int val1,val2;
  int irot[NAFIT]; // the best rotomer number for each side chain amino acid
  float lmax,tempfact;
  float denmult,dentrans; // scaling parameters for density
  int useset; // if = 0 do not use set in tracetobr
  float phiallow[4][2],psiallow[4][2]; // another set of phi/psi allowed angle range
  int movecb;  // if = 0 do not move Cb in Powell refinement of rotomer
  int ringang; //if = 0 only use the default ring angle, else try others
  float maxcbmove; // Max allowed cb movement in Powell refinement
  float valueweight[4]; // parameters used to weight value in rotfunccb,etc. valueweight[0,1] = wieght for nbad, badset
  int maxngap; // Maximum allowed baseline (den0.8) length of gap to be used for extending
  int gapnum; // gap used in connectnum() varies from cnum-geo.gapnum to cnum+geo.gapnum
  float rankmax,rankgap; // parameters for when to best gap ranking is good enougn
  float deldist;  // lenght of den1.delx vector - corresponds to resolution of map
  int gnmax;  // max dist to move out from current point when checking distance from bone
  float maxline;  // maximum allowed values for short fits of ave line dist in bonedistgeo
  int minnres; //Minimum number of resdues required to try and test gap in extendfit()
  int maxugap;  //Maximum gap allowed to connect unassigned fits
  int resn; // Just used to pass arbitray residue number
  float sidemaindist; // Minimum allowed distance for long sidechain atoms to mainchain 
  int miniden; // the minimum value of bone iden used when looking for connection
  int trmetfit[2];// = 1 if trace meets another fit - used as criteria for cutting end
  int addseed;  // add to ipt when set srand48
  float frcon; // fraction of atoms connected using bones as test
};


void maketestpep(struct geometry& geo,float pos[NAFIT][NAA][3],int numaa,
        float phi,float psi);
void maketestpep(struct geometry& geo,float pos[NAFIT][NAA][3],int numaa,
        float phi,float psi);
void makealanine2(struct geometry& geo,float pos[NAFIT][NAA][3],int ifind,int aanum,
          float phi,float psi); 
void makecisalanine2(struct geometry& geo,float pos[NAFIT][NAA][3],int ifind,int aanum,
          float phi,float psi); 
void makerevalanine2(struct geometry& geo,float pos[NAFIT][NAA][3],int ifind,int aanum,
          float phi,float psi); 
void makecisrevalanine2(struct geometry& geo,float pos[NAFIT][NAA][3],int ifind,int aanum,
          float phi,float psi); 
void findcb(struct geometry& geo,float n[3],float ca[3],
               float c[3],float cb[3]);
void addcb2(struct geometry& geo,float pos[NAFIT][NAA][3],int ifit);
void setup(struct geometry& geo); /* creates array delta and epsilon*/
void assignlevconst(struct geometry& geo);
void findconstants(struct geometry& geo);
void findallowedang(struct griddata& den1,struct geometry& geo);
void findphipsierr(struct griddata& den1,struct geometry& geo);
void cajunct(float m[3][3],float det,float ang,float u1[3],float v1[3],float p1[3],
     float u2[3],float v2[3],float p2[3],float w1[3],float w2[3]);
void cajunctn(int t[],int ix,float r[5][3],float m[3][3],float det,float ang,
     float u1[3],float v1[3],float p1[3],
     float u2[3],float v2[3],float p2[3],float w1[3],float w2[3]);
void findphipsi(float pos[NAFIT][NAA][3],int i,float *phi,float *psi);
void findphipsifit(struct pepfit& fit,int i,float *phi,float *psi);
void  printphipsi(struct geometry geo);
void addcbn(struct geometry& geo,float pos[NAFIT][NAA][3],int ifit,int iatom);
void classifyphipsi(struct geometry& geo,float phi,float psi,
   int *nclass,int type[],float radius[]);
void nclassifyphipsi(struct geometry& geo,float phi,float psi,
   int *type);
