#include "dist.h"

#define MAXPDBS 2000

/* ===========================================
   Parameters for Weeding
   c compare volumes w/ & w/o HETATMS
   s atoms with surface area above threshhold
   t atoms touching surface atoms
   H HETATMS
   A atoms touching HETATMs

   sthresh [n] is the surface area threshhold
=========================================== */

int cflag;
int sflag;
int tflag;
double sthresh;
int hflag;
int wflag;
int Hflag;
int Aflag;

int debug;  /* for regular debugging */
int debug1; /* to discover why my volumes so small */
int test;   /* to figure out why this is crashing */

#define nResStats 21

typedef struct {
  char resname[4];
  Statistic *V_restat;    /* statistic for residue volumes */
  Statistic *SA_restat;   /* statistic for residue exposed surface areas */

  int atnum;
  Statistic ** atstats;    /* array for pointers to atom volumes */
  Statistic ** SA_atstats; /* array for pointers to atom exposed surface areas */

} Residue_Stats;

Residue_Stats ResStats[nResStats];
/* =============================== */
int nTypFiles;
FILE **TypFiles;
/* =============================== */


void SetupStatsForResidue(char *resfile)
{
  int ii, jj, non;
  char pdbnm[10], type[5];
  FILE *fp;

  /* Setting up Big Statistics Array */
  if( (fp = fopen(resfile,"r")) == NULL )  {
    fprintf(stderr, "Could not find %s.", resfile);
    exit(1);
  }

  for(ii=0;ii<nResStats;ii++)  {
    if(fscanf(fp, "%s %d",ResStats[ii].resname , &ResStats[ii].atnum) == EOF) {
      fprintf(stderr, "residue file ended prematurely, %d\n", ii);
      exit(33);  
    }

    /* for residue statistics */
    ResStats[ii].V_restat = (Statistic *)malloc(sizeof(Statistic));
    ResStats[ii].SA_restat = (Statistic *)malloc(sizeof(Statistic));
    InitStats(ResStats[ii].V_restat); InitStats(ResStats[ii].SA_restat);

    /* For atom stats */
    ResStats[ii].atstats = (Statistic **)calloc(ResStats[ii].atnum, sizeof(Statistic *));
    ResStats[ii].SA_atstats = (Statistic **)calloc(ResStats[ii].atnum, sizeof(Statistic *));
    
    for(jj=0, non=0;jj<ResStats[ii].atnum;jj++)  {
      if(fscanf(fp, "%s %s",pdbnm, type) == EOF) {
	fprintf(stderr, "residue file ended prematurely, %d\n", ii);
	fprintf(stderr, "within residue %s\n", ResStats[ii].resname);
	exit(34);  
      }
      if(type[0] == '=') {
	non++;
	continue;
      }
      
      ResStats[ii].atstats[jj] = (Statistic *)malloc(sizeof(Statistic));
      InitStats(ResStats[ii].atstats[jj]);
      ResStats[ii].atstats[jj]->name = (char *)calloc(4, sizeof(char));
      sprintf(ResStats[ii].atstats[jj]->name,"%-3s", pdbnm);

      ResStats[ii].SA_atstats[jj] = (Statistic *)malloc(sizeof(Statistic));
      InitStats(ResStats[ii].SA_atstats[jj]);
      ResStats[ii].SA_atstats[jj]->name = (char *)calloc(4, sizeof(char));
      sprintf(ResStats[ii].SA_atstats[jj]->name,"%-3s", pdbnm);
    }
    ResStats[ii].atnum -= non; 
  }
  fclose(fp);

  if(debug)  {
    for(ii=0;ii<nResStats;ii++)  {
      fprintf(stderr,"%s %d\n", ResStats[ii].resname ,ResStats[ii].atnum);
      for(jj=0;jj<ResStats[ii].atnum;jj++)  {
	fprintf(stderr,"\t%s\n", ResStats[ii].atstats[jj]->name);
      }
    }
  }
}


/* This code was edited from full-dump-polyhedra.main.c ------------------------------*/
/* mark's code -----------------------------------------------------------------------*/
static int VertexOK(vertex * v) 
{
  if (v == NULL) 
    return 0;
  return 1 ;
}

static int NbrOk (neibour_descr *nd) {
  if (nd == NULL) 
    return -1;
  if (nd->atom == NULL) 
    return -1 ;
  return nd->atom->number ;
}


void Res_FullDumpPoly (atom_record *atom,neibour_descr *neibours,int neibour_count)
{
  int ii,jj,n_vertices_count,ll, done ;  /* last one is my counter */
  double volume = 0.0 ; 

  /* to weed out if atoms has surface area */
  if(sflag && atom->b > sthresh)   return;
  /* or is a HETATM */
  if(Hflag && atom->ident[0]=='H') return;

  for (ii=0 ; ii < neibour_count ; ii++) {
    neibour_descr * nd = &neibours[ii] ; 

    if (nd->flag) {
      vertex * n_vertices[MAX_VERTEX] ;

      /* weeds out if neighbor has surface area */
      if(tflag && nd->atom->b > sthresh) return;
      /* or is an HETATM */
      if(Aflag && nd->atom->ident[0]=='H') return;

      n_vertices_count= 0;
      for (jj = 0 ; jj < vertices_count; jj++) {
	vertex *ve = &vertices[jj];
	if (ve->a1 == nd || ve->a2 == nd || ve->a3 == nd ) {
	  ve->flag = 0 ;
	  n_vertices[n_vertices_count++] =ve ;
	}
      }
      
      if (n_vertices_count < 3) {
	STDERR("FullDumpPoly(): Less than 3 vertices skipping neighbor...");
	exit(1);
      } else  { 
	double FaceVolume = 0;
	vertex * FaceVertices[MAX_VERTEX];
	vertex *first_vertex = n_vertices[0] ;
	vertex *second_vertex, *third_vertex, *previous_second = first_vertex;
	neibour_descr *xnd = first_vertex->a1;
	
	if (xnd == nd)
	  xnd = first_vertex->a2 ;

	for(jj =1 ; jj< n_vertices_count ; jj++)
	  if (n_vertices[jj]->a1 == xnd ||
	      n_vertices[jj]->a2 == xnd ||
	      n_vertices[jj]->a3 == xnd ) {
	    second_vertex = n_vertices[jj] ;
	    break ;
	  }
	first_vertex->flag = 1;
	second_vertex->flag = 1 ;
	
	FaceVertices[0] = first_vertex;
	FaceVertices[1] = second_vertex ;

	for (ll = 2 ; ll < n_vertices_count ; ll++) {
	  neibour_descr *next_nd = second_vertex->a1 ;
	  
	  if (next_nd == nd)
	    if  (second_vertex->a2 == xnd )
	      next_nd = second_vertex->a3 ;
	    else next_nd = second_vertex->a2 ;
	  else if (next_nd == xnd)
	    if  (second_vertex->a2 == nd )
	      next_nd = second_vertex->a3 ;
	    else next_nd = second_vertex->a2 ;
	  
	  for (jj = 1  ; jj < n_vertices_count ; jj++) {
	    third_vertex = n_vertices[jj] ;
	    if (!third_vertex->flag &&
		(third_vertex->a1 == next_nd ||
		 third_vertex->a2 == next_nd ||
		 third_vertex->a3 == next_nd))
	      break ;
	  }
	  if (jj == n_vertices_count) {
	    warning(ARG_TYPE_ATOM,"Did not find third_vertex", atom) ;
	    return;
	  }
	  { 
	    double t_volume = 
	      tet_volume(atom,first_vertex, second_vertex, third_vertex);
	    FaceVolume += t_volume ;
	    
	  }
	  xnd = next_nd;
	  third_vertex->flag = 1;
	  FaceVertices[ll] = third_vertex ;
	  second_vertex = third_vertex ;
	} /* for ll */ 
	volume += FaceVolume ;
      } /* else */
    }
  }

  atom->volume=volume;

  if(debug) fprintf(stderr,"# individual atoms %s %s -", atom->resname, atom->name);
  for(ii=0;ii<nResStats;ii++)  {
    done=0;
    if(STREQ(atom->resname, ResStats[ii].resname))  {
      if(debug) fprintf(stderr," %s", ResStats[ii].resname);
      for(jj=0;jj<ResStats[ii].atnum;jj++)  {
	if(STREQ(atom->name, ResStats[ii].atstats[jj]->name))  {
	  if(debug) fprintf(stderr," %s %f", ResStats[ii].atstats[jj]->name, volume);
	  IncrementStats(ResStats[ii].atstats[jj], volume);	   
	  done=1; break;
	}
      }
      if(done) break;
      else {
	fprintf(stderr, "Didn't type the following atom: "); write_pdb_record(stderr, atom, 0);
      }
    }
  }
  if(debug) fprintf(stderr,"\n");
}

void SaveAtomVertices (int n,atom_record * a, vertex *v1, vertex *v2) 
{}

void LoopOverNeighbors (atom_record *atom,neibour_descr *neibours,int neibour_count) 
{}
/* end mark's code --------------------------------------------------------------------*/
     

void DisulfideBondedOrNot (file_records *f)
{
  int ii, jj;

  for(ii=0;ii<f->atomnum;ii++)  {
    atom_record *a = f->atoms[ii];

    /* if SG covalently bonded, make entire residue CSS */
    if(!strcmp("CYS", a->resname) && !strcmp("SG ", a->name) && a->flags & SULPHUR_BRIDGE) {
      fprintf(stderr,"------- changing %s %d to CSS because of disulfide --------\n", a->resname, a->resnum);
      for(jj=0;jj<f->atomnum;jj++)  {
        atom_record * b = f->atoms[jj];
	if(b->resnum == a->resnum && !strcmp("CYS", b->resname)) strcpy(b->resname, "CSS");
      }
    }
  }
}

/* figures out residue volume and increments stats */
void figure_res_V (file_records *f)
{
  int ii, jj, kk, atres=f->atoms[0]->resnum;
  double resV=0.0, resSA=0.0;

  for(kk=0;kk<f->atomnum;kk++)  {
    atom_record *a = f->atoms[kk];

    if(a->ident[0]=='H' || a->extraflag=='B') continue; /* takes out HETATMS */

    if(debug) fprintf(stderr,"# residue volume %s %s -", a->resname, a->name);

    if(atres != a->resnum)  {
      if(debug) printf("---- totals v-%6.2lf sa-%6.2lf ----\n", resV, resSA);
      if(!resSA && resV > 0.0) {  /* gets at resnumber change and checks if there is surface area */
	for(ii=0;ii<nResStats;ii++)  {
	  if(STREQ(f->atoms[kk-1]->resname, ResStats[ii].resname))  {
	    if(debug) fprintf(stderr," %s %f", ResStats[ii].resname, resV);
	    IncrementStats(ResStats[ii].V_restat, resV);
	    break;
	  }
	}
      }
      resV=0.0; resSA=0.0; atres = a->resnum;
    }
    else if(kk == f->atomnum-1) { /* gets the last residue */
      resV+=a->volume; resSA+=a->b;
      if(!resSA && resV > 0.0)  {
	for(ii=0;ii<nResStats;ii++)  {
	  if(STREQ(a->resname, ResStats[ii].resname))  {
	    if(debug) fprintf(stderr," %d %s %f", ii, ResStats[ii].resname, resV);
	    IncrementStats(ResStats[ii].V_restat, resV);
	    break;
	  }
	}
      }
    }

    if(debug) {
      fprintf(stdout, "%s %d %s %d %s v-%6.2lf sa-%6.2lf\n", 
	     f->name, a->resnum, a->resname, a->number, a->name, a->volume, a->b);
      fprintf(stderr, "\n");
    }
    if(debug1) {
      fprintf(stdout, "%s %d %s %d %s v-%6.3lf sa-%6.3lf\n", 
	     f->name, a->resnum, a->resname, a->number, a->name, a->volume, a->b);
    }
    if(a->volume > 0)  { resV += a->volume;  resSA += a->b; }
    else               { resSA += a->b + 1;                 } /* takes into account HETATMS or weeding */
  }
}

/* increments atom and residue stats for exposed surface area */
void figure_exposed_SA (file_records *f)
{
  int ii, jj, kk, atres=f->atoms[0]->resnum;
  double resSA=0.0;

  for(kk=0;kk<f->atomnum;kk++)  {
    atom_record *a = f->atoms[kk];

    if(a->ident[0]=='H' || a->extraflag=='B') continue; /* takes out HETATMS and b conformers */

#ifdef OFF
    if(debug1)
      fprintf(stdout, "%s %d %s %d %s v-%6.3lf sa-%6.3lf b-%6.3lf\n", 
	      f->name, a->resnum, a->resname, a->number, a->name, a->volume, a->surface, a->b);
#endif

    if(debug) fprintf(stderr,"# %d surface area  %s %d %s %d -", nResStats, a->resname, a->resnum, a->name, a->number);

    /* for atoms exposed SA */
    for(ii=0;ii<nResStats;ii++)  {
      int done = 0;
      if(a->surface > 0.0  && STREQ(a->resname, ResStats[ii].resname))  {
	for(jj=0;jj<ResStats[ii].atnum;jj++)  {
	  if(STREQ(a->name, ResStats[ii].SA_atstats[jj]->name))  {
	    if(debug)  fprintf(stderr, " %d %d - %s %e", ii, jj, ResStats[ii].SA_atstats[jj]->name, a->surface);
	    IncrementStats(ResStats[ii].SA_atstats[jj], a->surface);
	    done=1; break;
	  }
	}
	if(done) break;
      }
    }
    if(debug) fprintf(stderr, "\n");

    /* for residues */
    if(atres != a->resnum && resSA > 0)  {  /* gets at resnumber change */ 
      for(ii=0;ii<nResStats;ii++)  {
	if(STREQ(f->atoms[kk-1]->resname, ResStats[ii].resname))  {
	  if(debug)  {
	    fprintf(stderr, "tSA %d %s %f\n", ii, ResStats[ii].resname, resSA);;
	  }
	  IncrementStats(ResStats[ii].SA_restat, resSA);
	  break; 
	}
      }
      resSA=0.0;  atres=a->resnum; 
    }
    else if(kk == f->atomnum-1)  {  /* gets last residue */
      resSA += a->surface;
      if(resSA > 0)  {
	for(ii=0;ii<nResStats;ii++)  {
	  if(STREQ(a->resname, ResStats[ii].resname))  { 
	    if(debug)  fprintf(stderr, "tSA %d %s %f\n", ii, ResStats[ii].SA_atstats[jj]->name, resSA);
	    IncrementStats(ResStats[ii].SA_restat, resSA);
	    break; 
	  }
	}
      }
    }
    resSA += a->surface; 
  }
}


/* Since the NN_Calculate_Volumes rewrites the surface value */
/* need to put the surface values into another register, b   */
/* also, zeros the volume register                           */
void put_surface_into_b (file_records *f)
{
  int ii;
  
  for(ii=0;ii<f->atomnum;ii++)  {
    atom_record *a = f->atoms[ii];
    a->b = a->surface; a->surface = 0.0;
    a->volume = 0.0;
  }
}



void SelectionOfAtoms(char *weed_level, double surf_thresh)
{
  int ii, slength;

  sflag = tflag = Hflag = hflag = wflag = Aflag = 0; sthresh = surf_thresh;

  if(strcmp("none", weed_level))  {
    slength = strlen(weed_level) +1;
    for(ii=0;ii<slength;ii++)  {
      fprintf(stderr, "Weed -- ");
      switch (weed_level[ii])  {
      case 'c':
	cflag = 1;
	fprintf(stderr, "Doing Statistics Based on change in volume w/ and w/o HETATMS\n");
	fprintf(stderr, "Weed -- ============= THIS IS DEFAULT =================\n");
	break;
      case 's':
	sflag = 1; 
	fprintf(stderr, "Protein Surface Atoms above %lf.\n", sthresh);
	break;
      case 't':
	tflag = 1;
	fprintf(stderr, "Protein Atoms touching surface protein atoms above %lf.\n", sthresh);
	break;
      case 'H':
	Hflag = 1;
	fprintf(stderr, "HETATMs.\n");
	break;
      case 'A':
	Aflag = 1;
	fprintf(stderr, "Atoms touching HETATMs.\n");
	break;
      case 'h':
	hflag = 1;
	fprintf(stderr, "taking HETATMs out of calculation\n");
	break;
      case 'w':
	wflag = 1;
	fprintf(stderr, "taking waters out of calculation\n");
	break;
      default:
	fprintf(stderr, "Don't have entry for %c.\n", weed_level[ii]);
	break;
      }
    }
  }
  if(!sflag && surf_thresh)  {
    sflag = 1; 
    fprintf(stderr, "Weed -- Protein Surface Atoms above %lf.\n", surf_thresh);
  }
}


/* list of pdb files to run over, */
/* allows ability to comment pdb out with a # */
/* automatically tags the directory of pdbs/ */
int get_pdbsfrom_stdin(char ***thepdbs, char *pdb_dir)  
{
  int ii, jj, mxSlength=0, nthepdbs=0, lpdb_dir;
  char line[100], tmp[MAXPDBS][50];

  lpdb_dir = strlen(pdb_dir);

  while (gets(line) != NULL) {
    sscanf(line,"%s", tmp[nthepdbs]);
    if(tmp[nthepdbs][0] == '#' || tmp[nthepdbs][0] == ' ') continue;
    ii = strlen(tmp[nthepdbs]); if(mxSlength<ii+1) mxSlength=ii+1;
    nthepdbs++;
    if(nthepdbs >= MAXPDBS)  {
      fprintf(stderr, "too many pdbs from stdin!\n");
      fprintf(stderr, "can only handle %d. exiting ......\n", MAXPDBS);
      exit(20);
    }
  }

  if(nthepdbs == 0)  {
    fprintf(stderr, "no files read in from stdin, exiting .....");
    exit(21);
  }

  fprintf(stderr, "will analyze the %d pdb files (longest filename %d characters)\n", 
	  nthepdbs, mxSlength);
 
  for(ii=0;ii<nthepdbs;ii++)  { 
    fprintf(stderr, "%s ", tmp[ii]);
    if((ii+1)%5==0) fprintf(stderr, "\n");
  }
  if(ii%5!=0) fprintf(stderr, "\n");
   
  *thepdbs = (char **)calloc(nthepdbs, sizeof(char *));
  mxSlength += lpdb_dir;
  for(jj=0;jj<nthepdbs;jj++) {
    (*thepdbs)[jj]=(char *)calloc(mxSlength, sizeof(char));
    sprintf((*thepdbs)[jj], "%s%s", pdb_dir, tmp[jj]);
  }

  return nthepdbs;
}


void New176Output (char *res_file, char *outname)
{
  int ii, jj;
  double SDavg;
  double Tnum=0.0, Tmean=0.0, Tsd, Ssd=0.0, Tsdp=0.0;
  char outfile[100];
  FILE *OF;

  /* outputting volumes of atoms */
  sprintf(outfile, "%s.a-vols", outname);
  if(NULL == (OF = fopen(outfile, "w"))) { fprintf(stderr, "Could not open %s.\n", outfile); exit(44); }
  fprintf(OF, "# atom's volumes :: for residue file %s\n", res_file);
  fprintf(OF, "   atom  Count  Mean volume    SD       SD(\%)    Minimum   Maximum\n");
  for(ii=0;ii<nResStats;ii++)  {
    fprintf(OF, "%s\n", ResStats[ii].resname);
    for(jj=0;jj<ResStats[ii].atnum;jj++)  {
      if(ResStats[ii].atstats[jj]->num == 0.0 || ResStats[ii].atstats[jj]->mean == 0.0)  {
	ResStats[ii].atstats[jj]->max = 0.0; ResStats[ii].atstats[jj]->min = 0.0; SDavg = 0.0;
      }
      else  {
	SDavg = (ResStats[ii].atstats[jj]->dev/ResStats[ii].atstats[jj]->mean)*100;

     }
      fprintf(OF, "   %s   %5.0lf   %7.3lf    %6.3lf    %6.3lf     %6.3lf     %6.3lf \n", 
	     ResStats[ii].atstats[jj]->name,
	     ResStats[ii].atstats[jj]->num, ResStats[ii].atstats[jj]->mean,
	     ResStats[ii].atstats[jj]->dev, SDavg,
	     ResStats[ii].atstats[jj]->min, ResStats[ii].atstats[jj]->max);
      
      Tnum  += ResStats[ii].atstats[jj]->num; 
      Tmean += ResStats[ii].atstats[jj]->mean;
      Ssd   += (ResStats[ii].atstats[jj]->dev*ResStats[ii].atstats[jj]->dev);
    }

    if(Tnum == 0.0 || Tmean == 0.0)  { Tsd = 0.0; Tsdp = 0.0; }
    else { Tsd = sqrt(Ssd); Tsdp = 100*(Tsd/Tmean); }
    fprintf(OF, "\n Total   %5.0lf   %7.3lf    %6.3lf    %6.3lf\n\n", Tnum, Tmean, Tsd, Tsdp);
    Tnum=0.0; Tmean=0.0; Ssd=0.0;
  }
  fclose(OF);

  /* outputting exposed surface area's of atoms */
  Tnum=0.0, Tmean=0.0, Tsd, Ssd=0.0, Tsdp=0.0;
  sprintf(outfile, "%s.a-surf", outname);
  if(NULL == (OF = fopen(outfile, "w"))) { fprintf(stderr, "Could not open %s.\n", outfile); exit(45); }
  fprintf(OF, "# atom's exposed surface area :: for residue file %s\n", res_file);
  fprintf(OF, "   atom  Count  Mean volume    SD       SD(\%)    Minimum   Maximum\n");
  for(ii=0;ii<nResStats;ii++)  {
    fprintf(OF, "%s\n", ResStats[ii].resname);
   for(jj=0;jj<ResStats[ii].atnum;jj++)  {
      if(ResStats[ii].SA_atstats[jj]->num == 0.0 || ResStats[ii].SA_atstats[jj]->mean == 0.0)  {
	ResStats[ii].SA_atstats[jj]->max = 0.0; ResStats[ii].SA_atstats[jj]->min = 0.0; SDavg = 0.0;
      }
      else  {
	SDavg = (ResStats[ii].SA_atstats[jj]->dev/ResStats[ii].SA_atstats[jj]->mean)*100;
	
     }
      fprintf(OF, "   %s   %5.0lf   %7.3lf    %6.3lf    %6.3lf     %6.3lf     %6.3lf \n", 
	     ResStats[ii].SA_atstats[jj]->name,
	     ResStats[ii].SA_atstats[jj]->num, ResStats[ii].SA_atstats[jj]->mean,
	     ResStats[ii].SA_atstats[jj]->dev, SDavg,
	     ResStats[ii].SA_atstats[jj]->min, ResStats[ii].SA_atstats[jj]->max);
      
      Tnum  += ResStats[ii].SA_atstats[jj]->num; 
      Tmean += ResStats[ii].SA_atstats[jj]->mean;
      Ssd   += (ResStats[ii].SA_atstats[jj]->dev*ResStats[ii].SA_atstats[jj]->dev);
    }

    if(Tnum == 0.0 || Tmean == 0.0)  { Tsd = 0.0; Tsdp = 0.0; }
    else { Tsd = sqrt(Ssd); Tsdp = 100*(Tsd/Tmean); }
    fprintf(OF, "\n Total   %5.0lf   %7.3lf    %6.3lf    %6.3lf\n\n", Tnum, Tmean, Tsd, Tsdp);
    Tnum=0.0; Tmean=0.0; Ssd=0.0;
  }
  fclose(OF);

  /* now outputting for the residues' volumes  */
  sprintf(outfile, "%s.res-vols", outname);
  if(NULL == (OF = fopen(outfile, "w"))) { fprintf(stderr, "Could not open %s.\n", outfile); exit(46); }
  fprintf(OF, "# buried residue volumes  :: for residue file %s\n", res_file);
  fprintf(OF, "   atom  Count  Mean volume    SD     Minimum     Maximum\n");
  for(ii=0;ii<nResStats;ii++)  {
    if(ResStats[ii].V_restat->num == 0.0 || ResStats[ii].V_restat->mean == 0.0)  {
      ResStats[ii].V_restat->max = 0.0; ResStats[ii].V_restat->min = 0.0; 
    }
    fprintf(OF, "   %s   %5.0lf   %7.3lf    %6.3lf    %7.3lf     %7.3lf \n", 
	    ResStats[ii].resname,
	    ResStats[ii].V_restat->num, ResStats[ii].V_restat->mean,
	    ResStats[ii].V_restat->dev,
	    ResStats[ii].V_restat->min, ResStats[ii].V_restat->max);
  }
  fclose(OF);

  /* now outputting for the residues' exposed surface area */
  sprintf(outfile, "%s.res-surf", outname);
  if(NULL == (OF = fopen(outfile, "w"))) { fprintf(stderr, "Could not open %s.\n", outfile); exit(47); }
  fprintf(OF, "# exposed residue surface areas :: for residue file %s\n", res_file);
  fprintf(OF, "   res   Count  Mean SurfAr    SD     Minimum     Maximum\n");
  for(ii=0;ii<nResStats;ii++)  {
    if(ResStats[ii].SA_restat->num == 0.0 || ResStats[ii].SA_restat->mean == 0.0)  {
      ResStats[ii].SA_restat->max = 0.0; ResStats[ii].SA_restat->min = 0.0; 
    }
    fprintf(OF, "   %s   %5.0lf   %7.3lf    %6.3lf    %7.3lf     %7.3lf \n", 
	    ResStats[ii].resname,
	    ResStats[ii].SA_restat->num, ResStats[ii].SA_restat->mean,
	    ResStats[ii].SA_restat->dev,
	    ResStats[ii].SA_restat->min, ResStats[ii].SA_restat->max);
  }
  fclose(OF);
}


void Old176Output (char *res_file)
{
  int ii, jj;
  double SDavg;

  printf("# for residue file %s\n", res_file);
  printf("# Res    atom  Count  Mean volume(A3) StdDev(A3)   StdDev(%%)   Maximum   Minimum\n\n");
  for(ii=0;ii<nResStats;ii++)  {
    for(jj=0;jj<ResStats[ii].atnum;jj++)  {
      if(ResStats[ii].atstats[jj]->num == 0.0)  {
	ResStats[ii].atstats[jj]->max = 0.0;
	ResStats[ii].atstats[jj]->min = 0.0;
	SDavg = 0.0;
      }
      else  {
	SDavg = (ResStats[ii].atstats[jj]->dev/ResStats[ii].atstats[jj]->mean)*100;
      }

      printf("  %s     %s  %5.1lf     %6.3lf        %6.3lf       %7.3lf     %6.3lf    %6.3lf \n", 
	     ResStats[ii].resname, ResStats[ii].atstats[jj]->name,
	     ResStats[ii].atstats[jj]->num, ResStats[ii].atstats[jj]->mean,
	     ResStats[ii].atstats[jj]->dev, SDavg,
	     ResStats[ii].atstats[jj]->max, ResStats[ii].atstats[jj]->min);
    }
    printf("\n");
  }
}



/* FOR ALL SELECTION ROUTINES, ALSO DESELECTING B CONFORMATION */
void FlagAllAtoms (file_records *ff)
{
  int ii;

  for(ii=0;ii<ff->atomnum;ii++)  {
    atom_record *b = ff->atoms[ii];
    if(b->extraflag == 'B') {
      b->b=-1.0;b->volume = 0.0;
      SWITCH_ON(b,FLAG_1_DESELECT|FLAG_2_DESELECT);
    }
    else  { 
      SWITCH_OFF(b,FLAG_1_DESELECT|FLAG_2_DESELECT);
    }
  }
}

void UnFlagHETATMs (file_records *ff)
{
  int ii;

  for(ii=0;ii<ff->atomnum;ii++)  {
    atom_record *b = ff->atoms[ii];
    if (b->ident[0] == 'H' || b->extraflag == 'B') {
      SWITCH_ON(b,FLAG_1_DESELECT|FLAG_2_DESELECT);
    } 
    else {
      SWITCH_OFF(b,FLAG_1_DESELECT|FLAG_2_DESELECT);
    }
  }
}

void UnFlagOnlyHOHatoms(file_records *ff)  
{
  int ii;

  for(ii=0;ii<ff->atomnum;ii++)  {
    atom_record *b = ff->atoms[ii];
    if (!strcmp("HOH", b->resname) || b->extraflag == 'B') {
      b->b = -1.0; b->volume = 0.00;
      SWITCH_ON(b,FLAG_1_DESELECT|FLAG_2_DESELECT);
    } else {
      SWITCH_OFF(b,FLAG_1_DESELECT|FLAG_2_DESELECT);
    }
  }
}

void put_volume_into_b (file_records *f)
{
  int ii;
  
  for(ii=0;ii<f->atomnum;ii++)  {
    atom_record *a = f->atoms[ii];
    a->b = a->volume; a->volume=0.0;
  }
}


void Compare_and_Statistics(file_records *f)
{
  int ii, jj, kk;

  for(kk=0;kk<f->atomnum;kk++)  {
    atom_record *atom = f->atoms[kk];

    if(atom->ident[0]=='H') continue; /* takes out HETATMS */

    if(debug1) {
      fprintf(stdout, "%s %d %s %d %s v_w/o_HET- %6.3lf v_w/_HET- %6.3lf sa- %6.3lf\n", 
	      f->name, atom->resnum, atom->resname, atom->number, atom->name, 
	      atom->volume, atom->b, atom->surface);
    }

    if(atom->b > 0.0 && atom->b  == atom->volume && atom->surface <= sthresh )  {
      if(debug) fprintf(stderr,"# %s %s %lf %lf -", atom->resname, atom->name, atom->b, atom->volume);
      for(ii=0;ii<nResStats;ii++)  {
	int done = 0;
	if(STREQ(atom->resname, ResStats[ii].resname))  {
	  if(debug) fprintf(stderr," %s", ResStats[ii].resname);
	  for(jj=0;jj<ResStats[ii].atnum;jj++)  {
	    if(STREQ(atom->name, ResStats[ii].atstats[jj]->name))  {
	      if(debug) fprintf(stderr," %s %f", ResStats[ii].atstats[jj]->name, atom->volume);
	      IncrementStats(ResStats[ii].atstats[jj], atom->volume);	   
	      done=1; break;
	    }
	  }
	  if(done) break;
	  else {
	    fprintf(stderr, "Didn't type the following atom: "); write_pdb_record(stderr, atom, 0);
	  }
	}
      }
      if(debug) fprintf(stderr,"\n");
    }
  }
}


void volume_stats(double surf_thresh, char *weed_level, int out_form, int method, 
		  char *types_file, char *res_file, char *pdb_dir, char *outname)
{
  int ii, nthepdbs;
  char ** thepdbs;
  
  nthepdbs = get_pdbsfrom_stdin(&thepdbs, pdb_dir);
  SetupStatsForResidue(res_file);
  SelectionOfAtoms(weed_level, surf_thresh);

  for(ii=0;ii<nthepdbs;ii++)  {
    file_records *f = open_and_read_pdb(thepdbs[ii],0);
    if(f==NULL)  {
      fprintf(stderr, "#### couldn't find %s %d\n", thepdbs[ii], ii);
      continue;
    }
    fprintf(stderr, "############========== now at %s %d ===========#############\n", thepdbs[ii], ii);
    DisulfideBondedOrNot (f);

    FlagAllAtoms(f);
    if (wflag) UnFlagOnlyHOHatoms(f);

    if(cflag)  {
      NN_Calculate_Volumes(f, method, 64.0, 8.0);
      put_volume_into_b(f);
      UnFlagHETATMs(f);
      NN_Calculate_Volumes(f, method, 64.0, 8.0);
      UnFlagOnlyHOHatoms(f);
      calculate_surface(f,1.4);
      Compare_and_Statistics(f);
    }
    else {
      /* calculate surface first, so that can weed them in res_fulldumppoly */
      if(hflag) UnFlagHETATMs(f);
      calculate_surface(f,1.4); 
      figure_exposed_SA (f);
      put_surface_into_b(f);
      CalcVHookFcn = Res_FullDumpPoly; 
      NN_Calculate_Volumes(f, method, 64.0, 8.0);
      figure_res_V(f);
    }
    free_file_records(f);
  }
 
  switch (out_form)  {
  case '1':
    Old176Output(res_file);
  default:
    New176Output(res_file, outname);
  }
}


main(int argc,char **argv)
{
  SAY_HI_FIRST ;

  int out_form = (int)DoubleAfterOption(argc, argv, "output_format", 0);
  int method = (int)DoubleAfterOption(argc, argv, "method", 2);
  double surf_thresh = DoubleAfterOption(argc, argv, "surf_thresh", 0.0);
  char * weed_level = StringAfterOption(argc, argv, "weed_level", "none");
  char * types_file = STRG_AFTER_OPT_W_ENV("types_file","$HOME/lib/default-atom-defs.dat") ;
  char * res_file = STRG_AFTER_OPT_W_ENV("res_file","$HOME/lib/default-residue-defs.dat") ;
  char * pdb_dir = StringAfterOption(argc, argv, "pdb_dir", "pdbs/");

  /* these next two get output name */
  char * tf = StringAfterOption(argc, argv, "types_file", "nVradii_out");
  char * outname = StringAfterOption(argc, argv, "outname", tf);

  debug  = ArgvMatch(argc, argv, "debug");
  debug1 = ArgvMatch(argc, argv, "debug1");
  test   = ArgvMatch(argc, argv, "test");

  initial_residues(types_file,res_file,0);

  if(debug) WriteAtomTypes();

  volume_stats(surf_thresh, weed_level, out_form, method, types_file, res_file, pdb_dir, outname);

  fprintf(stderr,"I'm done!\n");
}