/*
************************************************************************
*
*   bib.C -
*
*   Copyright (c) 1995
*
*   ETH Zuerich
*   Institut fuer Molekularbiologie und Biophysik
*   ETH-Hoenggerberg
*   CH-8093 Zuerich
*
*   All Rights Reserved
*
*   Date of last modification : 95/09/15
*   Pathname of SCCS file     : /export/home3/cb/garant-1.0/src/SCCS/s.bib.C
*   SCCS identification       : 1.4
*
************************************************************************
*/
/**************************************************************************/
/*                              bib.cc                                    */
/*                                                                        */
/* classes to access text files                                           */
/**************************************************************************/

#include <stdio.h>
#include <string.h>
#include <math.h>
#include "global.h"
#include "log.h"		/* general log file                 */
#include "wert.h"		/* general measured value           */
#include "str.h"		/* tokenize string                  */

#include "bib.h"

/* path for libraries */
char LIB_PATH[]=BASEDIR "/lib/";

		// number of atoms in a equivalence class
#define MAXEQUIVALENCE 10

		// number of basic fragments
#define MAXFRAG 100	

		// number of subfragments of a mollecule
#define MAXAA 1000

		// number of different atom types
#define MAXATOMTYPE 40

		// number of bomds per atom
#define MAXATBOND 4

		// number of different atom names
#define MAXATNAMES 200

		// number of different secondary structure types
#define MAXSECTYPE 10

class fragT;
class subfragT;
class atomNamesT;

/* Data structures that hold information read from library */

class secShiftT {		/* 2nd structure chemical shifts     */
public:
   Wert w[MAXSECTYPE];		/* shifts                            */
   char *names[MAXSECTYPE];     /* name of secondary structure type  */
   int n;
};

class equivAtNamesT {		/* equivalent atomnames              */
public:
   atomNamesT **equiv;		/* list of equiv atoms               */
   int max;			/* max. size of list                 */
   int nr;			/* nr. of equiv atoms                */
   void add(atomNamesT *);
   int find(char *);		
   int isEquiv(atomNamesT *);	/* is the given name in the list     */
};

class atomNamesT {		/* all the allowed atomnames         */
public:
   char name[MAXNAME];
   int typeNr;
   equivAtNamesT equivAtNam;
};

class atomtypeT {		/* atom types                       */
public:
   char type[MAXNAME];
};

class atomT {			/* atoms as defined in the ecepp.lib */
public:
   int nr;			/* number within fragment            */
   int bondto[MAXATBOND];	/* bound atoms                       */
   int nrBd;			/* nr of bound atoms                 */
   int ssAtom;			/* id for compatible atoms           */
   atomNamesT *name;		/* name                              */
   char *type;			/* type                              */
   float shift;                 /* statisticaly expected shift       */
   float delshift;		/* variance of expected shift        */
   float prob;			/* probability to observe the atom   */
   secShiftT secShift;		/* 2nd structure chem. shifts        */
};

#define MAGNETIC 1
#define CHEMICAL 2


class equivalenceT {		/* pseudo atoms as defined in ecepp.lib */
public:
   atomNamesT *name;	     	/* name of pseudo atom               */
   char *type;			/* type of pseudo atom               */
   int magnchem;		/* MAGNETIC or CHEMICAL equivalence  */
   int nrAt;			/* number of atoms                   */
   int at[MAXEQUIVALENCE];	/* atom numbers within fragment      */
   int nr;			/* atom number of pseudoatom         */
};

class equivFragsT {		/* equivalent fargment types        */
public:
   int max;			/* max. size of list                */
   int nr;			/* number of entries                */
   fragT **equiv;		/* list of equiv fragment types     */
   equivFragsT() {equiv = 0; max = 0; nr = 0; }
   void add(fragT *);		/* add frag. type to list           */
   int isEquiv(fragT *);	/* is the given fragment in the list*/
};

class fragT {			/* fragment as defined in ecepp.lib  */
public:
   char name[MAXNAME];		/* name of the fragment              */
   atomT *at;			/* atoms of the fragment             */
   int nrAt;			/* number of atoms defined           */
   int firstAt;			/* first atom belonging to the frag. */
   int lastAt;			/* last atom belonging to the frag.  */
   equivalenceT *eq;		/* definition of pseudoatoms         */
   int nrEq;			/* number of pseudo atoms            */
   equivFragsT equivFrags;	/* list of equivalent fragments      */
   int equivNr;
   fragT() {nrAt=0; nrEq=0; }
};

class subfragT {		/* subfragment of a sequence         */
public:
   int molnr;			/* number of the mollecule           */
   int subnr;			/* number of frag. within mollecule  */
   int externnr;		/* number for naming the fragment    */
   int type;			/* fragment belongs to ...           */
				/*    0 ... destination fragment     */
				/*    1 ... subfragments             */
   fragT *frag;			/* fragment type                     */
   int offAt;			/* offset to add to frag. atom nr.   */
};

class coT {			/* coherence file entries            */
public:
   int nr;			/* number of the coherence           */
   Wert shift;			/* chemical shift of the coherence   */
   float score;			/* local socre of the coherence      */
   char name[MAXNAME];		/* name of the coherence             */
   int nrExt;			/* number of the subfrgamnet         */
   subfragT *refFrag;		/* internal reference into frag. list*/
   char flag;			/* unused flag (EASY)                */
};

class Bib {
   atomtypeT atomTypes[MAXATOMTYPE]; 	/* all the atom types        */
   int nrAtTypes;
   atomNamesT atomNames[MAXATNAMES]; 	/* all the possible names    */
   int nrAtNam;
   fragT frag[MAXFRAG];		/* fragment types                    */
   int nrFrag;
   coT liCo[MAXCOHERENCE];	/* list of valid coherences          */
   int nrCo;			/* number of coherences              */
   subfragT subs[MAXAA];	/* all the subfragments              */
   int nrSub;
   int readlib(char *);		/* read ecepp.lib                    */
   int readEquivAt(char *);	/* read lib. of equiv atom names     */
   int readEquivFrags(char *);	/* read lib. of equiv fragment names */
   int readmol(char *,int);	/* read fragments (dest. or subfrgs) */
   int readcolist(char *, int); /* read co's (dest. or subfrags)     */
   int readSecShifts(char *);   /* read lib. with 2nd struc. shifts  */
public:
   // friend Request;
   Bib();
   friend class BibAccess;
   friend class AtomRequest;
   friend class FragRequest;
   friend class BondRequest;
   friend class PseudRequest;
   friend class FreqRequest;
   friend class PeakRequest;
   friend class CoordRequest;
   friend class AssRequest;
   friend class FragNrRequest;
   friend class FragNameRequest;
   friend class SecShiftRequest;
   friend class CoWriter;
   friend class PeakWriter;
   friend class AssWriter;
   int equiv(int,int);		/* return methyl pseudo atom             */
   atomNamesT *addAtName(char *);/* adds name to list; returns reference */
				/* into atom names list                  */
   int getAtNameNr(char *);     /* gets the ref. number for an atom name */
				/* -1 for unknown name                   */
   void equivAtNamStart(char *, char *, int &, int &);
   void equivAtNamNext(char *, int &, int &);
				/* Iteration over all atom names         */
				/* equivalent to the provided atom name  */
				/* uses the integers as iteration vars   */
   char *addAtType(char *);	/* adds type to list returns name        */
   fragT *getRefFrag(char *);	/* gets pointer to fragment corresponding*/
				/* given name                            */
   int nrCoGet(char *,int);	/* get the number of the coherence with  */
				/* given name and number                 */
   int refCoGet(int);
				/* get index of co with given number     */
};

static Bib bib;

/* Implementation */

Bib::Bib()
{
   nrAtTypes = 0;
   nrAtNam = 0;
   nrFrag = 0;
   nrCo = 0;
   nrSub = 0;
}

int Bib::readlib(char *fn)
{
   FILE *fp;
   char s1[MAXLINE],s2[MAXLINE];
   int n1,i,j,k;
   int res = 1;
   int res2 = 1;
   int read;
   int nf=-1;
   Log log;
   char str[MAXLINE];
   char name[MAXLINE];
   char type[MAXLINE];
   char magn[MAXLINE];

   fp=fopen(fn,"r");
   if(fp == NULL) {
      sprintf(str,"\n... can not open library '%s'",fn);
      log.w(str);
      res2 = 0;
   } else {
      while(fgets(s1,MAXLINE,fp)!=NULL) {
         if(strncmp(s1,"ATOMTYPES",7)==0) {
            res &= (sscanf(s1,"%s%d",s2,&n1) == 2);
            for(i=0;i<n1;i++)
               res &= (fscanf(fp,"%*d%s%*f%*d%*d\n",atomTypes[i].type)==1);
            nrAtTypes = n1;
         }
         if(strncmp(s1,"RESIDUE",7)==0 ||
            strncmp(s1,"STARTGRP",7)==0 ||
            strncmp(s1,"TAILRES",7)==0 ||
            strncmp(s1,"HEADRES",7)==0) {
            nf++;
            frag[nf].equivNr = 0;
            res &= (sscanf(s1,"%*10s%5s%5d%5d%5d%5d%*5d%5d\n",
                              frag[nf].name,
                              &n1, &frag[nf].nrAt, &frag[nf].firstAt,
                              &frag[nf].lastAt,
                              &frag[nf].nrEq) >= 5);
            for(i=0;i<n1;i++) res &= (fgets(s1,MAXLINE,fp)!=NULL);

            frag[nf].at = new atomT[frag[nf].nrAt];
            for(i=0;i<frag[nf].nrAt;i++) {
               res &= (fgets(s1,MAXLINE,fp)!=NULL);
               read = sscanf(s1,
                      "%*4d %5s%5s%*5d%*f%*f%*f%*f%d%d%d%d%*d%f%f%f%*d\n",
                              name,
                              type,
                              &frag[nf].at[i].bondto[0],
                              &frag[nf].at[i].bondto[1],
                              &frag[nf].at[i].bondto[2],
                              &frag[nf].at[i].bondto[3],
                              &frag[nf].at[i].shift,
                              &frag[nf].at[i].delshift,
                              &frag[nf].at[i].prob);
               frag[nf].at[i].name = addAtName(name);
               frag[nf].at[i].type = addAtType(type);
	       frag[nf].at[i].secShift.n=0;
               if(read<9) {
                  frag[nf].at[i].shift = BAD_PPM;
                  frag[nf].at[i].delshift = 0.0;
                  if(frag[nf].at[i].bondto[0] != 0)
                     frag[nf].at[i].prob = 0.0;
                  else
                     frag[nf].at[i].prob = 1.0;
               }
               res &= (read>=6);
               res &= (frag[nf].at[i].name != 0);
               res &= (frag[nf].at[i].type != 0);
            }

            if(frag[nf].nrEq>0) frag[nf].eq = new equivalenceT[frag[nf].nrEq];
            else frag[nf].eq=0;
            for(i=0;i<frag[nf].nrEq;i++) {
               res &= (fgets(s1,MAXLINE,fp)!=NULL);
               res &= (sscanf(s1,"%*4d %5s%5s%5d%10s",name,type,
                              &frag[nf].eq[i].nrAt,magn)==4);
               for(k=0;k<frag[nf].nrAt;k++)
                  if(strcmp(name,frag[nf].at[k].name->name)==0)
                        frag[nf].eq[i].nr=k;
               frag[nf].eq[i].name = addAtName(name);
               frag[nf].eq[i].type = addAtType(type);
               if(strcmp(magn,"MAGN")==0) frag[nf].eq[i].magnchem = MAGNETIC;
               else frag[nf].eq[i].magnchem = CHEMICAL;
               for(j=0;j<frag[nf].eq[i].nrAt;j++) {
                  res &= (sscanf(&s1[30+j*10],"%10s",s2)==1);
                  for(k=0;k<frag[nf].nrAt;k++)
                     if(strcmp(s2,frag[nf].at[k].name->name)==0)
                        frag[nf].eq[i].at[j]=k;
               }
               fscanf(fp,"\n");
            }
            if(res==0) {
               sprintf(str,
                    "Error in library at fragment %s.\n",frag[nf].name);
               log.w(str);
               res=1;res2=0;
            }
         }
      }
   }
   fclose(fp);
   nrFrag =  nf+1;
   return(res2);
}

int Bib::readSecShifts(char *fn)
{
   FILE *fp;
   int res;
   Log log;
   char str[MAXLINE];
   char s1[MAXLINE];
   char atnam[MAXLINE],resnam[MAXLINE],secnam[MAXLINE],t1[MAXLINE], t2[MAXLINE];
   int iat,iSub,iSh;
   atomT *atom;
   float x, v;

   fp=fopen(fn,"r");
   if(fp == NULL) {
      sprintf(str, "\n ... can not open library %s\n",fn);
      log.w(str);
      res = 0;
   } else {
      res = 1;
      while(!feof(fp)) {
         fgets(s1,MAXLINE,fp);
         strip(s1," \n");
         if(s1[0] != '#' && s1[0] != '\0') {
            getToken(s1,atnam," :");
            getToken(s1,resnam," :");
            atom=0;
            for(iSub=0;iSub<nrFrag && !atom;iSub++) {
               if(strcmp(resnam,frag[iSub].name)==0) {
                  for(iat=frag[iSub].firstAt-1;iat<frag[iSub].lastAt-1;iat++) {
                     if(strcmp(atnam,frag[iSub].at[iat].name->name)==0) {
		        atom=&(frag[iSub].at[iat]);
		     }
                  }
               }
            }
            if(atom == 0) {
               sprintf(str,
                  "\n... nWARNING: unknown atom '%s %s'",atnam, resnam);
               log.w(str);
            } else {
               getToken(s1,secnam," :");
               getToken(s1,t1," :");
               getToken(s1,t2," :");
               while(secnam[0] != '\0' && 
	          sscanf(t1, "%f", &x)==1 &&
		  sscanf(t2, "%f", &v)==1) {
	          iSh=atom->secShift.n;
	          if(iSh<MAXSECTYPE) {
	             atom->secShift.names[iSh] = new char[strlen(secnam)+2];
	             strcpy(atom->secShift.names[iSh], secnam);
	             atom->secShift.w[iSh].x=x;
	             atom->secShift.w[iSh].v=v;
	             atom->secShift.n++;
	          } else {
	             log.w("\nWARNING: cannot handle more than\n");
                     sprintf(str,"   %d different types of secondary structure\n",MAXSECTYPE);
                     log.w(str);		  
	          }
	          getToken(s1,secnam," :");
                  getToken(s1,t1," :");
                  getToken(s1,t2," :");
               }
	       if(secnam[0] != '\0') {
	          log.w("\nWARNING: could not read beyond\n");
                  sprintf(str,"   %4s %5s : %s",atnam, resnam, secnam);
                  log.w(str);		  		
               }
	    }
         }
      }
   }
   return res;
}

int Bib::readEquivAt(char *fn)
{
   FILE *fp;
   int res;
   Log log;
   char str[MAXLINE];
   char s1[MAXLINE];
   char token[MAXLINE];
   int first;			/* name referenced by first may be used */
   int equiv;			/* instead of names ref. by equiv       */

   fp=fopen(fn,"r");
   if(fp == NULL) {
      log.w("\n ... can not open library");
      res = 0;
   } else {
      res = 1;
      while(!feof(fp)) {
         fgets(s1,MAXLINE,fp);
         strip(s1," \n");
         if(s1[0] != '#' && s1[0] != '\0') {
            getToken(s1,token," :");
            first = getAtNameNr(token);
            if(first == -1) {
               (void)addAtName(token);
               first = getAtNameNr(token);
            }
            getToken(s1,token," :");
            while(token[0] != '\0' && first != -1) {
               equiv = getAtNameNr(token);
               if(equiv == -1) {
                  sprintf(str,
                      "\n... warning: unknown atom name '%s'",
                      token);
                  log.w(str);
               } else {
                  atomNames[first].equivAtNam.add(&atomNames[equiv]);
               }
               getToken(s1,token," :");
            }
         }
      }
   }
   return res;
}

int Bib::readEquivFrags(char *fn)
{
   FILE *fp;
   int res;
   Log log;
   char str[MAXLINE];
   char s1[MAXLINE];
   char token[MAXLINE];
   fragT *equiv[MAXFRAG];	/* list of equivalent fragments         */
   int nrEquiv;			/* number of equivalent fragments       */
   int equivNr;			/* number of the equivalence class      */
   int i,j;

   fp=fopen(fn,"r");
   if(fp == NULL) {
      sprintf(str,"\n ... can not open library");
      log.w(str);
      res = 0;
   } else {
      res = 1;
      equivNr=0;
      while(!feof(fp)) {
         fgets(s1,MAXLINE,fp);
         strip(s1," \n");
         if(s1[0] != '#' && s1[0] != '\0') {
            getToken(s1,token," :");
            nrEquiv = 0;
            equivNr++;
            while(token[0] != '\0') {
               equiv[nrEquiv] = getRefFrag(token);
               equiv[nrEquiv]->equivNr = equivNr;
               if(equiv[nrEquiv] == 0) {
                  sprintf(str,
                      "\n... warning: unknown fragment name '%s'",
                      token);
                  log.w(str);
               } else {
                  nrEquiv++;
               }
               getToken(s1,token," :");
            }
            for(i=0;i<nrEquiv;i++)
               for(j=i+1;j<nrEquiv;j++) {
               equiv[i]->equivFrags.add(equiv[j]);
               equiv[j]->equivFrags.add(equiv[i]);
            }
         }
      }
   }
   return res;
}

int Bib::readmol(char *fn,int type)
{
   FILE *fp;
   char s1[MAXLINE];
   char str[MAXLINE];
   static int nrSubFrags=-1;
   static int nrMol=0;
   int nrExt=0;
   int fragnr;
   int res = 1;
   int warned=0;
   static int nr=-1;
   static int offAtom = 0;
   int i;
   Log log;
   int status;		/* 0 ... res name expected */
			/* 1 ... res name read     */
			/* 2 ... frag nr. read     */
			/* 3 ... mapping read      */

   fp=fopen(fn,"r");
   nrMol++;
   status=0;
   if(fp == NULL) {
      sprintf(str,"\n... can not open sequence '%s'",fn);
      log.w(str);
      res = 0;
   } else {
      while(fscanf(fp,"%s",s1)==1) {
         if(s1[0] == '#') fgets(s1,MAXLINE,fp);
         else if(s1[0] == '|') {
            nrMol++;
            status=0;
         }
         else if(sscanf(s1,"%d",&nrExt) == 1 && nr >= 0) {
            if(status==1) {subs[nr].externnr = nrExt; status=2; }
            if(status==2) status=3;
         }
         else {
            status=1;
            for(i=0;i<nrFrag;i++) {
               if(strcmp(frag[i].name,s1)==0) fragnr=i;
            }
            if(fragnr<nrFrag) {
               if(!warned && (strcmp("CYS",s1)==0 || strcmp("CYSS",s1)==0)) {
                  log.w("\nWARNING: Expected CB shift depends on whether\n");
                  log.w("           CYS (reduced) or CYSS (oxidized) is used\n");
                  warned=1;
               }
               if(nr>=MAXAA) {
                  sprintf(str,
                     "\n... too many fragments defined - discarting %s %d",
                     s1,nrExt+1);
                  log.w(str);
               } else {
                  nr++;
                  subs[nr].frag = &frag[fragnr];
                  subs[nr].subnr = ++nrSubFrags;
                  subs[nr].externnr = ++nrExt;
                  subs[nr].molnr = nrMol;
                  subs[nr].offAt = offAtom;
                  subs[nr].type = type;
                  offAtom += frag[fragnr].lastAt - frag[fragnr].firstAt + 1;
               }
            } else {
               res=0;
            }
         }
      }
   }
   fclose(fp);
   nrSub = nr+1;
   sprintf(str,"\n... %d fragments with %d atoms defined",nrSub,offAtom);
   log.w(str);
   return(res);
}

int Bib::readcolist(char *name, int typeFrag)
{
   FILE *fp;
   char s1[MAXLINE];
   char str[MAXLINE];
   int error = 0;
   int fatalError = 0;
   int end = 0;
   int nr = 0;
   int i,count;
   Log log;

   fp=fopen(name,"r");
   if(!fp) fatalError =1;
   while(fp != NULL && !end && nr < MAXCOHERENCE) {
      error = 0;
      if(feof(fp)) end = 1;
      else {
         fgets(s1,MAXLINE,fp);
         fatalError = ferror(fp);
         if(!error) {
            if(s1[0] != '\0' && s1[0] != '#') {
               /* read new or old format */
               error = 
                  (sscanf(s1,"%d %f %f %s %d",
                         &liCo[nr].nr, &liCo[nr].shift.x, &liCo[nr].shift.v,
                         liCo[nr].name, &liCo[nr].nrExt) != 5 &&
                   sscanf(s1,"%d %f %f %e %s %d %c",
                         &liCo[nr].nr, &liCo[nr].shift.x, &liCo[nr].shift.v,
                         &liCo[nr].score, liCo[nr].name, &liCo[nr].nrExt,
                         &liCo[nr].flag) != 7
                   );
               count = 0;
               liCo[nr].refFrag = 0;
               for(i=0;i<nrSub;i++) {
                  if(subs[i].externnr == liCo[nr].nrExt &&
                     subs[i].type == typeFrag)
                  {
                     count++;
                     liCo[nr].refFrag = &subs[i];
                  }
               }
               if(count == 0) {
                  if(typeFrag == DESTFRAG) {
                     sprintf(str,
                     "\n... WARNING: no destination fragment with number %d",
                     liCo[nr].nrExt);
                     log.warn(10,str);
                  } else {
                     sprintf(str,
                     "\n... WARNING: no subfragment with number %d",
                     liCo[nr].nrExt);
                     log.warn(10,str);
                  }
               } else if(count > 1) {
                  if(typeFrag == DESTFRAG) {
                     sprintf(str,
                     "\n... WARNING: %d destination fragment with number %d",
                     count,liCo[nr].nrExt);
                     log.warn(10,str);
                  } else {
                     sprintf(str,
                     "\n... WARNING: %d subfragment with number %d",
                     count,liCo[nr].nrExt);
                     log.warn(10,str);
                  }
               }
               if(!error) nr++;
            }
         }
         if(nr>=MAXCOHERENCE) {
            sprintf(str,
                "\n... WARNING more than %d coherences, some are discarded",
                 MAXCOHERENCE);
            log.warn(-1,str);
            error = 1;
         }
      }
   }
   fclose(fp);
   nrCo = nr;
   return(!fatalError);
}

int Bib::equiv(int nr, int fnr)
{
   int i,j;
   int iEq= -1;

   for(i=0; i<subs[fnr].frag->nrEq && iEq<0; i++)
      for(j=0; j<subs[fnr].frag->eq[i].nrAt && iEq<0; j++)
         if(nr == subs[fnr].frag->eq[i].at[j]) iEq = i;
   return(iEq);
}

atomNamesT *Bib::addAtName(char *n)
{
   int i;
   atomNamesT *found=0;
   Log lo;
   static int written = 0;

   for(i=0;i<nrAtNam && !found;i++) {
      if(strncmp(n,atomNames[i].name,MAXNAME) == 0) {
         found = &atomNames[i];
      }
   }
   if(!found) {
      if(nrAtNam < MAXATNAMES) {
         strncpy(atomNames[nrAtNam].name,n,MAXNAME);
         found = &atomNames[nrAtNam];
         nrAtNam++;
      } else {
         if(!written) {
            lo.w("Error: too many different atom names in library\n");
            written = 1;
         }
      }
   }
   return found;
}

int Bib::getAtNameNr(char *n) 
{
   int i;
   int found=-1;

   for(i=0;i<nrAtNam && found==-1;i++) {
      if(strncmp(n,atomNames[i].name,MAXNAME) == 0) found = i;
   }
   return found;
}

void Bib::equivAtNamStart(char *templ, char *result, int &curAtNam,
                          int &curEquiv)
{
   curAtNam = getAtNameNr(templ);
   curEquiv = 0;
   strcpy(result,templ);
}

void Bib::equivAtNamNext(char *result, int &curAtNam, int &curEquiv)
{
   result[0] ='\0';
   if(curAtNam != -1) {
      if(curEquiv < atomNames[curAtNam].equivAtNam.nr)
      {
         strcpy(result,
                atomNames[curAtNam].equivAtNam.equiv[curEquiv]->name);
         curEquiv++;
      }
   }
}

fragT *Bib::getRefFrag(char *n)
{
   int i;
   fragT *found=0;

   for(i=0;i<nrFrag && found==0;i++) {
      if(strncmp(n,frag[i].name,MAXNAME) == 0) found = &frag[i];
   }
   return found;
}


char *Bib::addAtType(char *n)
{
   int i;
   char *found=0;
   Log lo;
   static int written = 0;

   for(i=0;i<nrAtTypes && !found;i++) {
      if(strncmp(n,atomTypes[i].type,MAXNAME) == 0) found = atomTypes[i].type;
   }
   if(!found) {
      if(nrAtTypes < MAXATOMTYPE) {
         strncpy(atomTypes[nrAtNam].type,n,MAXNAME);
         found = atomTypes[nrAtNam].type;
         nrAtTypes++;
      } else {
         if(!written) {
            lo.w("Error: too many different atom types in library\n");
            written = 1;
         }
      }
   }
   return found;
}

int Bib::nrCoGet(char *atnam,int nrRes)
{
   int i;
   int res = -1;

   if(nrRes != 0 || atnam[0] != '\0') {
      for(i=0;i<nrCo && res == -1;i++) {
         if(liCo[i].nrExt == nrRes && strcmp(atnam,liCo[i].name) == 0)
            res = liCo[i].nr;
      }
   }
   return res;
}

int Bib::refCoGet(int searchedCo)
{
   int i;
   int res = -1;

   for(i=0;i<nrCo && res == -1;i++) {
      if(liCo[i].nr == searchedCo) res=i;
   }
   return res;
}

void equivFragsT::add(fragT *toAdd)
{
   int i;
   fragT **newList;

   if(nr == max) {
      newList = new fragT *[max+10];
      for(i=0;i<max;i++) newList[i]= equiv[i];
      delete equiv;
      equiv = newList;
      max += 10;
   } 
   equiv[nr]=toAdd;
   nr++;
}

int equivFragsT::isEquiv(fragT *test)
{
   int i;
   int found = 0;

   for(i=0;i<nr && !found;i++) if(equiv[i]==test) found = 1;
   return found;
}


void equivAtNamesT::add(atomNamesT *toAdd)
{
   int i;
   atomNamesT **newList;

   if(nr == max) {
      newList = new atomNamesT *[max+10];
      for(i=0;i<max;i++) newList[i]= equiv[i];
      delete equiv;
      equiv = newList;
      max += 10;
   } 
   equiv[nr]=toAdd;
   nr++;
}

int equivAtNamesT::isEquiv(atomNamesT *namAt)
{
   int i;
   int res = 0;
   for(i=0;i<nr && !res;i++) if(equiv[i] == namAt) res=1;
   return res;
}

/* Access to the library */

int BibAccess::fragLoad(char *fn, int type)
{
   return bib.readmol(fn,type);
}

int BibAccess::coLoad(char *fn, int type)
{
   return bib.readcolist(fn,type);
}

int BibAccess::readlib(char *nameLibFrag)
{
   int res;
   Log log;
   char str[MAXLINE];

   sprintf(str,"\n... reading fragment definition library '%s'\n",
                     nameLibFrag);
   log.w(str);
   res = bib.readlib(nameLibFrag);
   if(!res) log.w(" - aborted\n");
   else log.w("\n");

   return res;
}

int BibAccess::readSecShifts(char *nameLib)
{
   int res;
   Log log;
   char str[MAXLINE];

   sprintf(str,"\n... reading library '%s' defining 2nd structure dependent shifts\n",
                     nameLib);
   log.w(str);
   res = bib.readSecShifts(nameLib);
   if(!res) log.w("    - aborted\n");
   else log.w("\n");

   return res;
}

int BibAccess::readEquivAt(char *nameLibAt)
{
   int res;
   Log log;
   char str[MAXLINE];

   sprintf(str,"\n... reading library '%s' defining equivalent atom names\n",
                     nameLibAt);
   log.w(str);
   res = bib.readEquivAt(nameLibAt);
   if(!res) log.w("    - aborted\n");
   else log.w("\n");
   return res;
}

int BibAccess::readEquivFrags(char *nameLibSs)
{
   int res;
   Log log;
   char str[MAXLINE];

   sprintf(str,
       "\n... reading library '%s' defining fragments with equivalent spin systems\n",
       nameLibSs);
   log.w(str);
   res = bib.readEquivFrags(nameLibSs);
   if(!res) log.w("    - aborted\n");
   else log.w("\n");

   return res;
}

int BibAccess::equivAtNam(char *t, char *n) 
{
   int nt;
   int res = 0;
   int i;

   nt = bib.getAtNameNr(t);
   if(nt != -1) {
      i=0;
      while(i<bib.atomNames[nt].equivAtNam.nr && !res) {
         if(strcmp(n,bib.atomNames[nt].equivAtNam.equiv[i]->name) == 0)
            res = 1;
         i++;
      }
   }
   return res;
}

int BibAccess::equivAt(int a1, int f1, int a2, int f2)
{
   int res = 1;

   a1 = a1 - bib.subs[f1].offAt - 2 + bib.subs[f1].frag->firstAt;
   a2 = a2 - bib.subs[f2].offAt - 2 + bib.subs[f2].frag->firstAt;

   if(bib.subs[f1].frag->at[a1].name != bib.subs[f2].frag->at[a2].name) {
      if(!bib.subs[f1].frag->at[a1].name->equivAtNam.isEquiv(
          bib.subs[f2].frag->at[a2].name)) res = 0;
   }

   if(res && f1 != f2) {
      if(bib.subs[f1].type == bib.subs[f2].type) {
         if(bib.subs[f1].frag != bib.subs[f2].frag &&
            !bib.subs[f1].frag->equivFrags.isEquiv(bib.subs[f2].frag))
            res=0;
      } else {
         res=0;
      }
   }

   return res;
}

int BibAccess::equivFrag(char *n1, char *n2)
{
   int res = 0;
   int i,j;

   if(strcmp(n1,n2) == 0) return 1;

   for(i=0;i<bib.nrFrag;i++) {
      if(strcmp(bib.frag[i].name,n1)==0) {
         for(j=0;j<bib.frag[i].equivFrags.nr;j++) {
            if(strcmp(bib.frag[i].equivFrags.equiv[j]->name,n2)==0) res=1;
         }
      }
   }

   return res;
}

int SecShiftRequest::get(int iSub,int nr,char *secName)
{
   int iAt, iSec, ifound;
   atomT *atom;

   iAt=nr-bib.subs[iSub].offAt+bib.subs[iSub].frag->firstAt-2;
   atom=&(bib.subs[iSub].frag->at[iAt]);
   ifound=0;
   for(iSec=0;iSec<atom->secShift.n;iSec++) {
      if(strcmp(atom->secShift.names[iSec], secName)==0){
	 ifound=1;
	 shift=atom->secShift.w[iSec];
      }
   }
   return ifound;
}

AtomRequest::AtomRequest() {eof=1;}

void AtomRequest::start() 
{
   iSub=0;
   eof=0;
   oldMolNr=-1;
   if(iSub>=bib.nrSub) eof = 1;
   else {
      iAt = bib.subs[iSub].frag->firstAt - 1;
      if(iAt<0) eof=1;
   }
   read();
}

void AtomRequest::read()
{
   int iEq;

   if(iSub>=bib.nrSub) eof=1;
   if(!eof) {
      iEq = bib.equiv(iAt,iSub);
      nr = bib.subs[iSub].offAt + iAt
           - bib.subs[iSub].frag->firstAt + 2;
      strcpy(type,bib.subs[iSub].frag->at[iAt].type);
      sprintf(name,"%s",bib.subs[iSub].frag->at[iAt].name->name);
      nrAa = bib.subs[iSub].externnr;
      molNr = bib.subs[iSub].molnr;
      shift = bib.subs[iSub].frag->at[iAt].shift;
      delshift = bib.subs[iSub].frag->at[iAt].delshift;
      if(iEq<0 || bib.subs[iSub].frag->eq[iEq].magnchem == CHEMICAL)
         prob = bib.subs[iSub].frag->at[iAt].prob;
      else
         prob = 1.0;
      typeFrag = bib.subs[iSub].type;
      nrFrag = iSub;
/* remove HN of aminoterminus *********************************************/
      if(strcmp(name,"HN") == 0 && oldMolNr != molNr) prob = 1.0;
      if(strcmp(name,"N")  == 0 && oldMolNr != molNr) prob = 1.0;
/**************************************************************************/

      iAt++;
      if(iAt>=bib.subs[iSub].frag->lastAt) {
         oldMolNr=molNr;
         iSub++;
         if(iSub<bib.nrSub) iAt = bib.subs[iSub].frag->firstAt - 1;
      }
   }
//   printf("debug: %5s %3d   %7.3f\n",name,nrAa,shift);
}

FragRequest::FragRequest() {eof=1;}

void FragRequest::start() 
{
   iSub=0;
   eof=0;
   if(iSub>=bib.nrSub) eof = 1;
   read();
}

void FragRequest::read()
{
   if(iSub>=bib.nrSub) eof=1;
   if(!eof) {
      name = bib.subs[iSub].frag->name;
      nrInt = iSub;
      nrExt = bib.subs[iSub].externnr;
      type = bib.subs[iSub].type;
      equivNr = bib.subs[iSub].frag->equivNr;
      iSub++;
   }
}

BondRequest::BondRequest() {eof=1;}

void BondRequest::start() 
{
   iSub=0;
   eof=0;
   if(iSub>=bib.nrSub) eof = 1;
   else {
      iAt = bib.subs[iSub].frag->firstAt - 1;
      if(iAt<0) eof=1;
   }
   iBd = 0;
   read();
}

void BondRequest::read()
{
   int iAt2;
   int same_mol=0;
   int iEq1,iEq2;

   do {
      if(iSub>=bib.nrSub) eof=1;
      if(!eof) {
         iEq1 = bib.equiv(iAt,iSub);
         if(iEq1<0 || bib.subs[iSub].frag->eq[iEq1].magnchem == CHEMICAL)
            nr1 = bib.subs[iSub].offAt + iAt
                  - bib.subs[iSub].frag->firstAt +2;
         else
            nr1 = bib.subs[iSub].offAt
                  + bib.subs[iSub].frag->eq[iEq1].nr
                  - bib.subs[iSub].frag->firstAt +2;

         iAt2 = bib.subs[iSub].frag->at[iAt].bondto[iBd] - 1;
         iEq2 = bib.equiv(iAt2,iSub);
         if(iEq2<0 || bib.subs[iSub].frag->eq[iEq2].magnchem == CHEMICAL)
            nr2 = bib.subs[iSub].offAt +iAt2
                 - bib.subs[iSub].frag->firstAt +2;
         else
            nr2 = bib.subs[iSub].offAt
                 + bib.subs[iSub].frag->eq[iEq2].nr
                 - bib.subs[iSub].frag->firstAt +2;
         molNr = bib.subs[iSub].molnr;

         same_mol = 1;
         if(iSub == 0 || bib.subs[iSub].molnr != bib.subs[iSub-1].molnr)
            if(iAt2 < bib.subs[iSub].frag->firstAt-1) same_mol = 0;
         if(iSub >= bib.nrSub ||
            bib.subs[iSub].molnr != bib.subs[iSub+1].molnr)
            if(iAt2 >= bib.subs[iSub].frag->lastAt) same_mol = 0;

         iBd++;
         if(iBd>=4 || nr2==0) {
            iBd = 0;
            iAt++;
            if(iAt>=bib.subs[iSub].frag->lastAt) {
               iSub++;
               if(iSub<bib.nrSub) iAt = bib.subs[iSub].frag->firstAt - 1;
            }
         }
      }
   } while((nr2==0 || !same_mol || nr1>=nr2) && !eof);
}

PseudRequest::PseudRequest() {eof=1;}

void PseudRequest::start() 
{
   iSub=0;
   eof=0;
   if(iSub>=bib.nrSub) eof = 1;
   iPseud = 0;
   iAt = 0;
   read();
}

void PseudRequest::read()
{
   atnr=-1;
   while (!eof && atnr==-1) {
      if(iSub>=bib.nrSub) { eof=1; }
      else if(bib.subs[iSub].frag->nrEq==0) {iSub++; iPseud=0;}
      else if(bib.subs[iSub].frag->eq[iPseud].magnchem == MAGNETIC) {
         iPseud++;
         if(iPseud>=bib.subs[iSub].frag->nrEq) {
            iPseud=0;
            iSub++;
         }
      } else {
         pseudnr = bib.subs[iSub].offAt
                   + bib.subs[iSub].frag->eq[iPseud].nr
                   - bib.subs[iSub].frag->firstAt +2;
         if(iAt < bib.subs[iSub].frag->eq[iPseud].nrAt)
           atnr = bib.subs[iSub].offAt
                  + bib.subs[iSub].frag->eq[iPseud].at[iAt]
                  - bib.subs[iSub].frag->firstAt +2;

         iAt++;
         if(iAt>=bib.subs[iSub].frag->eq[iPseud].nrAt) {
            iAt = 0;
            iPseud++;
            if(iPseud>=bib.subs[iSub].frag->nrEq) {
               iPseud=0;
               iSub++;
            }
         }
      }
   }
}

FreqRequest::FreqRequest() {eof=1; }

void FreqRequest::start()
{
   iCo =0;
   read();
}

void FreqRequest::read()
{
   while(iCo < bib.nrCo && bib.liCo[iCo].refFrag == 0) iCo++;
   if(iCo >= bib.nrCo) eof = 1;
   else {
      strcpy(name,bib.liCo[iCo].name);
      nrExt = bib.liCo[iCo].refFrag->externnr;
      shift = bib.liCo[iCo].shift;
      iCo++;
      eof = 0;
   }
}

PeakRequest::PeakRequest() {eof=1;}

void PeakRequest::start(char *name,int d,float noise,int *p)
{
   Log log;
   char str[MAXLINE];
   int dimFile;
   int i;

   dim = d;
   nl = noise;
   for(i=0;i<dim;i++) perm[i] = p[i];

   peaknr = 0;
   fp=fopen(name,"r");
   if(fp == 0) {
      eof = 1;
   } else {
      eof = 0;
   }

   if(!eof) {
      if(fscanf(fp,"# Number of dimensions %1d\n",&dimFile) != 1) {
         eof = 1;
         fclose(fp);
      }
   }

   if(!eof && dimFile != dim) {
      sprintf(str,"\n... ERROR: expected %d dimensional peaks",dim);
      log.w(str);
      eof = 1;
      fclose(fp);
   }

   if(!eof) {
      if(fscanf(fp,"%s",line)!=1) {
         eof = 1;
         fclose(fp);
      }
   }

   if(eof) {
      sprintf(str,"\n... ERROR: Could not read peaklist '%s'",name);
      log.w(str);
   } else {
      this->read();
   }
}

void PeakRequest::read()
{
   int res = 0;
   while(!eof && !res) {
      res = readline();
   }
}

int PeakRequest::readline()
{
   int res = 1;
   int i;
   Log log;
   char str[MAXLINE];

   if(feof(fp)) {
      eof = 1; res = 0;
   } else {
      if(line[0] == '#') {
         fgets(line,MAXLINE - 1,fp);
         fscanf(fp,"%s",line);
         res = 0;
      } else {
         ++peaknr;
// peak number read from file -->
         res &= (sscanf(line,"%d",&nr)==1);

// peak number is equal to peak count -->
//         nr = peaknr;

         for(i=0;i<dim;i++) res &= (
                 fscanf(fp,"%f",&w[perm[i]]) == 1);

         res &= (fscanf(fp,"%d %*s %e %e %c %*d",
                            &color,&vol.x,&vol.v,&intflag) == 4);

         for(i=0;i<dim;i++) res &= (
                 fscanf(fp,"%d",&nrCo[perm[i]]) == 1);

         res &= (fscanf(fp,"%*d") == 0);

         if(!res && !feof(fp)) {
            sprintf(str,
               "\n... WARNING: Difficulties reading %dth peak",
               peaknr);
            log.warn(10,str);
         }
         prob = 1.0;
         fscanf(fp,"%s",line);
         if(line[0] == '#') {
            fgets(str,MAXLINE - 1,fp);
            strncpy(comment,str,MAX_COMMENT - 2);
            comment[MAX_COMMENT-2]='\n';
            comment[MAX_COMMENT-1]='\0';
            fscanf(fp,"%s",line);
         } else {
            comment[0] = '\0';
         }
      }
   }
   return(res);
}

CoordRequest::CoordRequest() {eof=1;}

void CoordRequest::start(char *name)
{
   Log log;
   char str[MAXLINE];
   int i;

   fp=fopen(name,"r");
   if(fp == 0) {
      eof = 1;
   } else {
      eof = 0;
   }

   if(!eof) for(i=0;i<3;i++) fgets(line,MAXLINE,fp);

   if(eof) {
      sprintf(str,"\n... ERROR: Could not open coordinate file '%s'",name);
      log.w(str);
   } else {
      this->read();
   }
}

void CoordRequest::read()
{
   int res = 0;
   while(!eof && !res) {
      res = readline();
   }
}

int CoordRequest::readline()
{
   int res = 1;

   if(feof(fp)) {
      eof = 1; res = 0;
   } else {
      fgets(line,MAXLINE,fp);
      res = (sscanf(line,"%*d %s %d %*s %f %f %f",name,&nrExt,
                &coord[0], &coord[1], &coord[2]) == 5);
   }
   return(res);
}

void AssRequest::read()
{
   int res = 0;
   while(!eof && !res) {
      res = readline();
      if(res) genCoNames();
   }
}

void AssRequest::genCoNames()
{
   int i;
   int nt;
   char str[MAXLINE];
   Log log;
   int coRef;

   for(i=0;i<dim;i++) {
      coRef=bib.refCoGet(nrCo[i]);
      if(coRef != -1 && bib.liCo[coRef].refFrag !=0) {
         nt = bib.getAtNameNr(bib.liCo[coRef].name);
         if(nt != -1) {
            names[i] = &bib.atomNames[nt];
            nrFrag[i] = bib.liCo[coRef].refFrag->externnr;
         } else {
            sprintf(str,"\n... WARNING unknown atom name %s",
               bib.liCo[coRef].name);
            log.warn(10,str);
            names[i] = 0;
            nrFrag[i] = -1;
         }
      } else {
         nrFrag[i] = -1;
         names[i] = 0;
         if(nrCo[i]>0) {
            sprintf(str,"\n... WARNING atom number %d is not defined",nrCo[i]);
            log.warn(-1,str);
         }
      }
   }
}

int AssRequest::equivName(int d, char *nam)
{
   int i;
   int res;

   if(names[d] == 0) {res = 1; }
   else if(strcmp(names[d]->name,nam) == 0) {res = 1; }
   else {
      i=0; res=0;
      while(i<names[d]->equivAtNam.nr && !res) {
         if(strcmp(nam,names[d]->equivAtNam.equiv[i]->name) == 0)
            res = 1;
         i++;
      }
   }
   return res;
}

char *AssRequest::getName(int d)
{
   return names[d]->name;
}

void FragNrRequest::start()
{
   Log log;
   char str[MAXLINE];
   char number[MAXLINE];

   eof = 0;
   getToken(name,resname," ");
   getToken(name,number," ");
   if(name[0] == '\0') {
      nrDef = 1; nameDef = 1;
      if(number[0] == '\0') {
         if(sscanf(resname,"%d",&nrExt) != 1) {
            nrDef = 0;
         } else {
            nameDef = 0;
            strcpy(number,resname);
         }
      }
   }
   if(nrDef) {
      if(sscanf(number,"%d",&nrExt) != 1) {
         sprintf(str,"\n... WARNING: '%s' is not a residue number",name);
         log.warn(10,str);
         eof=-1;
      }
   }
   if(nameDef && resname[0] == '*') nameDef = 0;

   curFrag = 0;

   read();
}

void FragNrRequest::read()
{
   int res=0;

   while(!res && curFrag<bib.nrSub) {
      res = 1;
      if(nrDef && bib.subs[curFrag].externnr != nrExt) res = 0;
      if(nameDef && strcmp(bib.subs[curFrag].frag->name,resname) != 0) res = 0;
      if(bib.subs[curFrag].type != type && type != DESTSUB) res = 0;
      curFrag++;
   }

   if(!res) eof = 1;
   else nr = curFrag-1;
}

void FragNameRequest::start()
{
   if(nr>=0 && nr<bib.nrSub) {
      strcpy(resname,bib.subs[nr].frag->name);
      nrExt = bib.subs[nr].externnr;
      sprintf(name,"%s %d",resname,nrExt);
      eof = 0;
   } else {
     eof = 1;
   }
}

OptPosRequest::OptPosRequest() {eof=1;}

void OptPosRequest::start(char *name)
{
   fp=fopen(name,"r");
   if(fp == 0) {
      eof = 1;
   } else {
      eof = 0;
   }

   read();
}

void OptPosRequest::read()
{
   int ok = 0;

   while(!eof && !ok) {
      fgets(line,MAXLINE,fp);
      eof = feof(fp);
      if(line[0] != '#') {
         ok = (sscanf(line,"%d %f",&nrPeak,&prob) == 2);
      }
   }
   if(eof) fclose(fp);
}

void Writer::writeCom(char *com)
{
   if(!eof) {
       fprintf(fp,"# %s\n",com);
   }
}

void CoWriter::start(char *name)
{

   nrCo = 0;
   if(eof != 1) {eof = -1; }
   else {
      fp = fopen(name,"w");
      if(!fp) eof=-1;
      else eof = 0;
   }
}

int CoWriter::write()
{
   nrCo++;

   if(fprintf(fp, "%4d %7.3f %6.3f %-5s %4d \n",
   	      nrCo, shift.x, shift.v,
   	      name, nrAa) < 0) eof = -1;
   return(nrCo);
}

void PeakWriter::start(char *name)
{
   eof = 0;
   fp = fopen(name,"w");
   if(!fp) eof=-1;
   if(eof != -1 &&
      fprintf(fp,"# Number of dimensions %1d\n",dim)<0) {
      eof = -1;
   }
   if(eof != 0 && fp) fclose(fp);
}

void PeakWriter::write()
{
   char coords[MAXLINE], cos[MAXLINE], str[MAXLINE];
   int i;
   int nrCos[DIMENSION];
   Log log;

   for(i=0;i<dim;i++) {
      nrCos[i] = bib.nrCoGet(coAtNam[i],coResNum[i]);
      if(nrCos[i] == -1 && coAtNam[i][0] != '\0') {
         sprintf(str,
            "\n... WARNING: number of coherence '%s %d' is not defined",
            coAtNam[i],coResNum[i]);
         log.warn(50,str);
      }
      if(nrCos[i] == -1) nrCos[i] = 0;
   }

   strcpy(cos,"");
   strcpy(coords,"");
   for(i=0;i<dim;i++) {
      sprintf(str,"%7.3f ",w[perm[i]].x);
      strcat(coords,str);
      sprintf(str,"%4d ",nrCos[perm[i]]);
      strcat(cos,str);
   }

   curPeak++;
   if(fprintf(fp,"%4d %s%1d %-9s %10.3e %9.2e %c %3d %s%1d\n",
              nrPeak,coords,color,"?",v.x,v.v,(int)intflag, 0,
              cos, 0) < 0) eof = -1;
   if(comment[0] != '\0')
      if(fprintf(fp,"#%s\n",comment) < 0) eof = -1;
   if(eof != 0 && fp) fclose(fp);
}

void OptPosWriter::start(char *name)
{
   fp = fopen(name,"w");
   if(!fp) eof=-1;
   else eof = 0;
   if(eof != 0 && fp) fclose(fp);
}


void OptPosWriter::write()
{
   if(fprintf(fp,"%d %f\n",nrPeak,prob) < 0) eof = -1;
   if(eof != 0 && fp) fclose(fp);
}

void AssWriter::start(char *name)
{
   int res,i;

   eof = 0;
   curPeak = 0;
   fp = fopen(name,"w");
   if(!fp) eof=-1;
   if(eof != -1) {
      res = 1;
      res &= (fprintf(fp,"Assignment file\n") > 0);
      res &= (fprintf(fp,"   Corresponding peaklist: null\n") > 0);
      res &= (fprintf(fp,"   Corresponding coherencelist: %s\n",namCoLi) > 0);
      res &= (fprintf(fp,"   Number of dimensions: %1d\n",dim) > 0);
      res &= (fprintf(fp,"   Uncertainties: ") > 0);
      for(i=0;i<dim;i++) {
         res &= (fprintf(fp," %7.3f",tolerance[i]) > 0);
      }
      res &= (fprintf(fp,"\n\n\n") >= 0);

      if(!res) eof = -1;
   }
   if(eof != 0 && fp) fclose(fp);
}

void AssWriter::startPeak()
{
   curPeak++;
   if(fprintf(fp,"# %5d\n",nrPeak) < 0) eof = -1;
}

void AssWriter::writeAssPos()
{
   char cos[MAXLINE], str[MAXLINE];
   int i;
   int nrCos[DIMENSION];
   Log log;

   for(i=0;i<dim;i++) {
      nrCos[i] = bib.nrCoGet(coAtNam[i],coResNum[i]);
      if(nrCos[i] == -1 && coAtNam[i][0] != '\0') {
         sprintf(str,
            "\n... WARNING: number of coherence '%s %d' is not defined",
            coAtNam[i],coResNum[i]);
         log.warn(10,str);
      }
   }

   strcpy(cos,"");
   for(i=0;i<dim;i++) {
      sprintf(str,"%4d ",nrCos[perm[i]]);
      strcat(cos,str);
   }

   curPeak++;
   if(fprintf(fp,"\t%s\n",cos)<0) eof=-1;
   if(eof != 0 && fp) fclose(fp);
}