#include <stdio.h>
#include <stdlib.h>
#include "libsx.h"
#include "globals.h"
#include "grid.h"
#include "graph.h"
#include "kraftwerk.h"
#include "post.h"
#include "module.h"
#include "geodata.h"

#define CANT_FIND 0
#define FOUND_IT 1

typedef struct {
  int seg;
  int n;
  int *nodes;
} segmentsOnNodes;

segmentsOnNodes *segnod;
int nseg;
float   xmin_orig;
float   xmax_orig;
float   ymin_orig;
float   ymax_orig;

int positionFilePointer(FILE *stream, char *tag);
void readEntry(char *file, char *tag, void (*readingFunction)() );
void  addSolutionMenus();

int thisType;
int eqDataItem;

void readMesh(FILE *stream)
{
  int i;

  if(stream != NULL){

    currentSolution = 0;
    SetToggleState(solButt[currentSolution], TRUE);
    NoOfMeshes = 1;

    fscanf(stream, "%i%i", &mySol[currentSolution].nnode, &mySol[currentSolution].elnum);
    if(mySol[currentSolution].elnum == 0)
      errorHandler("Invalid data file!");

      if(allocateMemory(currentSolution, mySol[currentSolution].elnum,
		     mySol[currentSolution].nnode) == 0){
	

    for( i = 0; i < mySol[currentSolution].nnode; i++)
      fscanf(stream,"%f%f", &mySol[currentSolution].coord[i][0] ,&mySol[currentSolution].coord[i][1]);
    for( i = 0; i < mySol[currentSolution].elnum; i++)
      fscanf(stream,"%i%i%i", &mySol[currentSolution].nodpek[i][0],&mySol[currentSolution].nodpek[i][1],&mySol[currentSolution].nodpek[i][2]) ;


/*
 * Renumber node numbers so they start from zero
 *
 */

    for(i = 0; i < mySol[currentSolution].elnum; i++)
      {
	mySol[currentSolution].nodpek[i][0] -= 1;
	mySol[currentSolution].nodpek[i][1] -= 1;
	mySol[currentSolution].nodpek[i][2] -= 1;
      }
  }
  }
  else
    errorHandler("Unable to open file");

}

void new_readPoints(FILE *stream)
{
  int i, j, n;
  float f1, f2;
  PTS_NODE *p;
  fscanf(stream, "%i", &n);

  for(i = 0; i < n; i++){
    fscanf(stream, "%i %f %f", &j, &f1, &f2);
    p = add_pnt(f1, f2);
    p->id = j;
  }

}
void new_readSegments(FILE *stream)
{
  PTS_NODE *pts1, *pts2;
  GEO_NODE *b;
  BND_NODE *s;
  int n, i, id;
  int p1, p2;

  fscanf(stream, "%i ", &n);
  for(i = 0; i < n; i++){
    fscanf(stream, "%i %i %i ", &id, &p1, &p2);
/* Now find the pointers associated with p1,p2 */
    pts1 = get_pts_pointer(p1);
    pts2 = get_pts_pointer(p2);
    FOR_ALL_BOUNDARIES(b){
      FOR_ALL_SEGMENTS(b, s){
	if(s->seg->id == id){
	  s->seg->pts1 = pts1;
	  s->seg->pts2 = pts2;
	  break;
	}
      }
    }
  }
    
}
void new_readBoundary(FILE *stream)
{
  int  j, id, n, dummy;
  GEO_NODE *g;

  fscanf(stream, "%i", &dummy);
  fscanf(stream, "%i %i", &j, &n);

  g = new_bnd();
  g->type = thisType;
  for(j = 0; j < n; j++){
    fscanf(stream, "%i ", &id);
    add_seg2(id, g);
  }
  
}


void readDomain(FILE *stream)
{
  fscanf(stream,"%f %f %f %f", &xmin, &xmax, &ymin, &ymax);
  xmin_orig = xmin;
  xmax_orig = xmax;
  ymin_orig = ymin;
  ymax_orig = ymax;

}
void readState(FILE *stream)
{
  fscanf(stream,"%i", &state);
}
void readModuleID(FILE *stream)
{
  int id;
  fscanf(stream,"%i", &id);
  module[id].state = state;
  setModule(id);
}


void readSolution(FILE *stream)
{
  int el, imesh, i;
  int slask;

  currentSolution = 0;
  SetToggleState(solButt[currentSolution], TRUE);

/* free old pointers */
  for( imesh = 0; imesh < NoOfMeshes; imesh++){
    free(mySol[imesh].nodpek);
    free(mySol[imesh].coord);
    free(mySol[imesh].u);
    free(mySol[imesh].f);
    free(mySol[imesh].vecX);
    free(mySol[imesh].vecY);
  }

  if(stream != NULL){
    fscanf(stream, "%i", &NoOfMeshes);
    for(imesh = 0; imesh < NoOfMeshes; imesh++){
      fscanf(stream, "%i", &mySol[imesh].elnum);
      fscanf(stream, "%i", &mySol[imesh].nnode);

      if(allocateMemory(imesh, mySol[imesh].elnum, mySol[imesh].nnode) == 0){

      for(el = 0; el < mySol[imesh].elnum; el++)
	fscanf(stream, "%i%i%i", &mySol[imesh].nodpek[el][0], &mySol[imesh].nodpek[el][1],
	                         &mySol[imesh].nodpek[el][2]);
      

      for(el = 0; el < mySol[imesh].elnum; el++)
	fscanf(stream, "%i%i%i", &slask, &slask, &slask);
      

      for(i = 0; i < mySol[imesh].nnode; i++){
	fscanf(stream, "%f", &mySol[imesh].coord[i][0]);
	fscanf(stream, "%f", &mySol[imesh].coord[i][1]);
      }

      for(i = 0; i < mySol[imesh].nnode; i++)
	fscanf(stream, "%f", &mySol[imesh].u[i]);
      for(i = 0; i < mySol[imesh].nnode; i++)
	fscanf(stream, "%f%f", &mySol[imesh].vecX[i], &mySol[imesh].vecY[i]);

      for(el = 0; el < mySol[imesh].elnum; el++){
	mySol[imesh].nodpek[el][0] -= 1;
	mySol[imesh].nodpek[el][1] -= 1;
	mySol[imesh].nodpek[el][2] -= 1;
      }
    } /* Next mesh */
    }
  }
  else
    errorHandler("Unable to open the solution file");
}
void readSolution2(FILE *stream)
{
  int el, imesh, i, j;
  int slask, ioff, ntot;
  char dummy;

  currentSolution = 0;
  SetToggleState(solButt[currentSolution], TRUE);
/* free old pointers */
  for( imesh = 0; imesh < NoOfMeshes; imesh++){
    free(mySol[imesh].nodpek);
    free(mySol[imesh].coord);
    free(mySol[imesh].u);
    free(mySol[imesh].f);
    free(mySol[imesh].vecX);
    free(mySol[imesh].vecY);
  }

  if(stream != NULL){
    fscanf(stream, "%i", &NoOfMeshes);
    for(imesh = 0; imesh < NoOfMeshes; imesh++){
      fscanf(stream, "%i", &mySol[imesh].elnum);
      fscanf(stream, "%i", &mySol[imesh].nnode);
      fscanf(stream, "%i %i ", &mySol[imesh].n, &ntot);

      if(allocateMemory(imesh, mySol[imesh].elnum, mySol[imesh].nnode) == 0){

      for(el = 0; el < mySol[imesh].elnum; el++)
	fscanf(stream, "%i%i%i", &mySol[imesh].nodpek[el][0], &mySol[imesh].nodpek[el][1],
	                         &mySol[imesh].nodpek[el][2]);
      

      for(el = 0; el < mySol[imesh].elnum; el++)
	fscanf(stream, "%i%i%i", &slask, &slask, &slask);
      

      for(i = 0; i < mySol[imesh].nnode; i++){
	fscanf(stream, "%f ", &mySol[imesh].coord[i][0]);
	fscanf(stream, "%f ", &mySol[imesh].coord[i][1]);
      }
      mySol[imesh].f  = malloc(ntot * mySol[imesh].nnode * sizeof(float));

      ioff = 0;
      for(j = 0; j < mySol[imesh].n; j++){
	mySol[imesh].ptr[j] = ioff;
	fscanf(stream, "%[^;]", mySol[imesh].name[j]);
	fscanf(stream, "%c", &dummy);
	fscanf(stream, " %i ", &mySol[imesh].type[j]);
	printf("         %s\n", mySol[imesh].name[j]);
	for(i = 0; i < mySol[imesh].nnode; i++){
	  if(mySol[imesh].type[j] == 1){
	    fscanf(stream, "%f ", &mySol[imesh].f[i + ioff]);
	  }
	  else{
	    fscanf(stream, "%f %f ", &mySol[imesh].f[i + ioff], 
                       &mySol[imesh].f[i + ioff + mySol[imesh].nnode]);
	  }
	}
	if(mySol[imesh].type[j] == 1)
	  ioff +=  mySol[imesh].nnode;
	else
	  ioff +=  2*mySol[imesh].nnode;
      }
      for(el = 0; el < mySol[imesh].elnum; el++){
	mySol[imesh].nodpek[el][0] -= 1;
	mySol[imesh].nodpek[el][1] -= 1;
	mySol[imesh].nodpek[el][2] -= 1;
      }
    } /* Next mesh */
    }
  }
  else
    errorHandler("Unable to open the solution file");
}
void readEqData(FILE *stream)
{
  int i, j;
  float x, y;
  int item;
  char dummy;

  item = eqDataItem;

  fscanf(stream, "%i ", &module[currentModule].data[item].n);
  fscanf(stream, "%i ", &module[currentModule].dataGroup[item].type);

  if(module[currentModule].dataGroup[item].type > 0){ 
    for(j = 0; j < module[currentModule].data[item].n; j++){
      if(module[currentModule].dataGroup[item].type == 2) /* elementwise */
	fscanf(stream, "%i ", &module[currentModule].data[item].value[j].el);
      else{                                               /* pointwise */
	fscanf(stream, "%f %f ", &x, &y);
	phys2pix(&module[currentModule].data[item].value[j].x,
		 &module[currentModule].data[item].value[j].y, x, y);
      }
      for(i = 0; i < module[currentModule].dataGroup[item].n; i++){
	fscanf(stream, "%f ", &module[currentModule].data[item].value[j].f[i]);
      }
    }
  }
  else{   /* analytical */
    for(i = 0; i < module[currentModule].dataGroup[item].n; i++){
      fscanf(stream, "%[^;]", module[currentModule].data[item].expr[i].text);
      fscanf(stream, "%c", &dummy);
    }
  }
}

void readNodesOnSegments(FILE *stream)
{
  int i, j, n, node;
  fscanf(stream, "%i", &nseg);
  free(segnod);
  segnod = malloc( nseg * sizeof(segmentsOnNodes));
  for(i = 0; i < nseg; i++){
    fscanf(stream, "%i", &n);
    segnod[i].seg = n;
    fscanf(stream, "%i", &n);
    segnod[i].n = n;
    segnod[i].nodes = malloc(n * sizeof(int));
    for(j = 0; j < n; j++){
      fscanf(stream, "%i", &node);
      segnod[i].nodes[j] = node;
    }
  }
}


void readBCData(FILE *stream)
{
  int ibc;
  int k;

  for(nvar = 0; nvar < module[currentModule].nvar; nvar++){
    fscanf(stream, "%i ", &module[currentModule].bc[nvar].n);
    for(ibc = 0; ibc < module[currentModule].bc[nvar].n; ibc++){
      fscanf(stream, "%i %c ", &module[currentModule].bc[nvar].segment[ibc],
	                       &module[currentModule].bc[nvar].type[ibc]);
      fscanf(stream, "%f ",    &module[currentModule].bc[nvar].value[ibc]);
      for(k = 0; k < module[currentModule].nbc; k++){
	if(module[currentModule].bc[nvar].type[ibc] == 
           module[currentModule].bc[k].tag[nvar]){
	   module[currentModule].bc[nvar].id[ibc] = k;
	 }
      }
    }
  }
}

void readNvar(FILE *stream)
{
  fscanf(stream, "%i", &nvar);
}

void loadFromFile(char *file)
{
  int i, j, n, x, y;
  int active[20], oldstate;
  float xp, yp;
  PTS_NODE *p;

  oldstate = state;
  nvar = 0;
  state = HAVE_NOTHING;
  clear_geometry(); 
  readEntry(file, "DOMAIN", readDomain);
  readEntry(file, "STATE", readState);
  readEntry(file, "MODULE ID", readModuleID);
  readEntry(file, "POINTS", new_readPoints);
  readEntry(file, "MESH", readMesh);
  readEntry(file, "NODES ON SEGMENTS", readNodesOnSegments);
  thisType = INNER;
  readEntry(file, "INNER BOUNDARY", new_readBoundary);
  thisType = OUTER;
  readEntry(file, "OUTER BOUNDARY", new_readBoundary);
  thisType = INTERNAL;
  readEntry(file, "INTERNAL BOUNDARY", new_readBoundary);
  readEntry(file, "SEGMENTS", new_readSegments);

  nvar = 0;
  readEntry(file, "BOUNDARY CONDITIONS", readBCData);
  for(i = 0; i < 20; i++)
    active[i] = FALSE;
  if(oldstate == HAVE_SOLUTION){
    for(i = 0; i < mySol[currentSolution].n + 1; i++)
      active[i] = GetMenuItemChecked(plotMenu.item[i+5]);
  }
  else
    active[0] = TRUE;
  rmSolutionMenus();
  readEntry(file, "SOLUTION", readSolution2);
  addSolutionMenus();
  for(i = 0; i < mySol[currentSolution].n + 1; i++)
      SetMenuItemChecked(plotMenu.item[i+5], active[i]);

  getDrawingScale();

  for( eqDataItem = 0; eqDataItem < module[currentModule].ngroups; eqDataItem++)
    readEntry(file, module[currentModule].dataGroup[eqDataItem].tag, readEqData);

  FOR_ALL_POINTS(p){
    xp = p->physicalX;
    yp = p->physicalY;
    phys2pix(&x, &y, xp, yp);
    p->x = x;
    p->y = y;
  }

}

/*
 *    Below are functions needed to save things to file
 *
 */
void writeNodesOnSegments(FILE *stream)
{
  int i, j;

  if( nseg > 0){
    fprintf(stream, "%s\n", "NODES ON SEGMENTS");
    fprintf(stream, "%i\n", nseg);
    for(i = 0; i < nseg; i++){
      fprintf(stream, "%i %i\n", segnod[i].seg, segnod[i].n);
      for(j = 0; j < segnod[i].n; j++){
	fprintf(stream, "%i\n", segnod[i].nodes[j]);
      }
    }
  }
}

void new_writeMeshData(FILE *stream)
{
/*
 * Dumps the boundary to the scratch file needed for meshgeneration
 *
 */
  int i, j, seg, ibc;
  int el, imesh;
  float x, y;
  BND_NODE *s;
  GEO_NODE *g;
  PTS_NODE *p;


/* Convert the X11-integer coordinates to physical coordinates */
  if(state <= HAVE_MESHSIZE)
    getDrawingScale();

/* Put to file */

  if(stream != NULL){
/*    Put a number on each point (needed for the fortran prog) */
/*
    printf("RENUMBERING ALL GEOMETRY STUFF ***************\n");
    i = 1;
    FOR_ALL_POINTS(p){
      p->id = i;
      i++;
    }
    i = 1;
    FOR_ALL_BOUNDARIES( g ){
      FOR_ALL_SEGMENTS(g, s){
	s->seg->id = i;
	i++;
      }
    }
*/
    fprintf(stream, "%s\n","POINTS");
    fprintf(stream, "%i\n", number_of_points()); 

    FOR_ALL_POINTS(p)
      fprintf(stream, "%i %f %f\n", p->id , p->physicalX, p->physicalY); 

    i = 1;
    FOR_ALL_BOUNDARIES( g ){
      fprintf(stream,"%s\n","SEGMENTS");
      fprintf(stream,"%i\n", number_of_segments(g));
      FOR_ALL_SEGMENTS(g, s)
	fprintf(stream,"%i %i %i \n", s->seg->id, s->seg->pts1->id, s->seg->pts2->id);
      if(g->type == INNER)
	fprintf(stream,"%s\n","INNER BOUNDARY");
      else if(g->type == OUTER)
	fprintf(stream,"%s\n","OUTER BOUNDARY");
      else
	fprintf(stream,"%s\n","INTERNAL BOUNDARY");
      fprintf(stream,"%i\n", 1);
      fprintf(stream,"%i %i \n", i, number_of_segments(g));
      FOR_ALL_SEGMENTS(g, s)
	fprintf(stream, "%i \n", s->seg->id);
      i++;
    }

    if(ih > 0){
      fprintf(stream,"%s\n", "SIZE");
      fprintf(stream,"%i\n", ih);
      for(i = 0; i < ih; i++){
	pix2phys(hx[i], hy[i], &x, &y);
	fprintf(stream,"%f %f %f\n", x, y, h[i]);
      }
    }
    if(state >= HAVE_MESH){
      imesh = 0;
      fprintf(stream,"%s\n", "MESH DATA");
      fprintf(stream, "%i %i\n", mySol[imesh].nnode, mySol[imesh].elnum);
      for(i = 0; i < mySol[imesh].nnode; i++){
	fprintf(stream, "%f\n", mySol[imesh].coord[i][0]);
	fprintf(stream, "%f\n", mySol[imesh].coord[i][1]);
      }
      for(el = 0; el < mySol[imesh].elnum; el++)
	fprintf(stream, "%i %i %i\n", 1 + mySol[imesh].nodpek[el][0],
		                      1 + mySol[imesh].nodpek[el][1], 
		                      1 + mySol[imesh].nodpek[el][2]);

      writeNodesOnSegments(stream);
    }
    
  }
  else
    errorHandler("Unable to create area.scratch");
}
void writeSolverData(FILE *stream)
{
  fprintf(stream, "%s\n", "REFINEMENT DATA");
  fprintf(stream, "%f %i '%s'\n", solverTol, solverLevels, solverNorm);
}


void writeEqData(FILE *stream)
{
  int item, i, j, ibc;
  float x, y;

  for(item = 0; item < module[currentModule].ngroups; item++){
    fprintf(stream, "%s\n", module[currentModule].dataGroup[item].tag);
    fprintf(stream, "%i\n", module[currentModule].data[item].n);
    fprintf(stream, "%i\n", module[currentModule].dataGroup[item].type);

    if(module[currentModule].dataGroup[item].type == 0){
	for(i = 0; i < module[currentModule].dataGroup[item].n; i++)
	  fprintf(stream, "%s;", module[currentModule].data[item].expr[i].text);
	fprintf(stream, "\n");
    }
    else{
      for(j = 0; j < module[currentModule].data[item].n; j++){
	if(module[currentModule].dataGroup[item].type == 1){
	  pix2phys(module[currentModule].data[item].value[j].x, 
		   module[currentModule].data[item].value[j].y , &x, &y);
	  fprintf(stream, "%f %f ", x, y);
	}
	else{
	  fprintf(stream, "%i ", module[currentModule].data[item].value[j].el);
	}
	for(i = 0; i < module[currentModule].dataGroup[item].n; i++)
	  fprintf(stream, "%f ", module[currentModule].data[item].value[j].f[i]);
	
	fprintf(stream, "\n");
      }
    }
  }
/*
  for( item = 0; item < numberOfEqData; item++){
    fprintf(stream, "%s\n", cddata[item].label);
    fprintf(stream, "%i\n", cddata[item].n);

    for(j = 0; j < cddata[item].n; j++){
      pix2phys(cddata[item].value[j].x, cddata[item].value[j].y, &x, &y);
      fprintf(stream, "%f %f ", x, y);

      for(i = 0; i < cddata[item].nvalues; i++){
	fprintf(stream, "%f ", cddata[item].value[j].f[i]);
      }
      fprintf(stream, "%s\n"," ");
    }
  }
*/
  fprintf(stream, "%s\n", "NUMBER OF VARIABLES");
  fprintf(stream, "%i\n", module[currentModule].nvar);

  fprintf(stream, "%s\n", "BOUNDARY CONDITIONS");
  for(j = 0; j < module[currentModule].nvar; j++){
    fprintf(stream, "%i\n", module[currentModule].bc[j].n);
    for(ibc = 0; ibc < module[currentModule].bc[j].n; ibc++){
      fprintf(stream, "%i %c\n", module[currentModule].bc[j].segment[ibc],
	                         module[currentModule].bc[j].type[ibc]);
      fprintf(stream, "%f\n",    module[currentModule].bc[j].value[ibc]);
    }
  }
}


void writeSolution(FILE *stream)
{
  int el, imesh, i, j, ioff, ntot;

  ntot = 0;
  for(i = 0; i < NoOfMeshes; i++){
    for(j = 0; j < mySol[i].n; j++)
      ntot += mySol[i].type[j];
  }

  if(stream != NULL){
    fprintf(stream, "%s\n", "SOLUTION");
    fprintf(stream, "%i\n", NoOfMeshes);
    for(imesh = 0; imesh < NoOfMeshes; imesh++){
      fprintf(stream, "%i\n", mySol[imesh].elnum);
      fprintf(stream, "%i\n", mySol[imesh].nnode);
      fprintf(stream, "%i %i\n ", mySol[imesh].n, ntot);
      for(el = 0; el < mySol[imesh].elnum; el++)
	fprintf(stream, "%i %i %i\n", 1 + mySol[imesh].nodpek[el][0],
		                      1 + mySol[imesh].nodpek[el][1], 
		                      1 + mySol[imesh].nodpek[el][2]);

      for(el = 0; el < mySol[imesh].elnum; el++)
	fprintf(stream, "%i %i %i\n", 1, 1, 1);

      for(i = 0; i < mySol[imesh].nnode; i++){
	fprintf(stream, "%f\n", mySol[imesh].coord[i][0]);
	fprintf(stream, "%f\n", mySol[imesh].coord[i][1]);
      }

      ioff = 0;
      for(j = 0; j < mySol[imesh].n; j++){
	mySol[imesh].ptr[j] = ioff;
	fprintf(stream, "%s;\n", mySol[imesh].name[j]);
	fprintf(stream, " %i\n ",mySol[imesh].type[j]);
	for(i = 0; i < mySol[imesh].nnode; i++){
	  if(mySol[imesh].type[j] == 1){
	    fprintf(stream, "%f\n", mySol[imesh].f[i + ioff]);
	  }
	  else{
	    fprintf(stream, "%f %f \n", mySol[imesh].f[i + ioff], 
                       mySol[imesh].f[i + ioff + mySol[imesh].nnode]);
	  }
	}
	if(mySol[imesh].type[j] == 1)
	  ioff +=  mySol[imesh].nnode;
	else
	  ioff +=  2*mySol[imesh].nnode;
      }

    } /* Next mesh */
  }
  else
    errorHandler("Unable to open the solution file");
}


void saveToFile(char *file)
{
  FILE *stream;

  getDrawingScale();
  stream = fopen(file, "w");

  if(stream != NULL){
    fprintf(stream, "#---------------------------------------------- \n");
    fprintf(stream, "# Data saved by FEMLAB\n");
    fprintf(stream, "# FEMLAB is brought to you by:\n");
    fprintf(stream, "#     Dep. of Mathematics\n");
    fprintf(stream, "#     Chalmers University of Technology\n");
    fprintf(stream, "#     Gothenburg Sweden\n");
    fprintf(stream, "#     (femlab@math.chalmers.se)\n");
    fprintf(stream, "# \n");
    fprintf(stream, "#     FEMLAB, the choice of generation X\n");
    fprintf(stream, "#---------------------------------------------- \n");
    fprintf(stream, "%s\n", "MODULE ID");
    fprintf(stream, "%i\n", module[currentModule].id);
    fprintf(stream, "%s\n", module[currentModule].name);
    fprintf(stream, "%s\n", "STATE");
    fprintf(stream, "%i\n", state);
    fprintf(stream, "%s\n", "DOMAIN");
    fprintf(stream,"%f %f %f %f\n", xmin, xmax, ymin, ymax);
    if(state >= HAVE_BOUNDARY){
      new_writeMeshData(stream);
    }
    if(state >= HAVE_EQDATA)
      writeEqData(stream);
    if(state >= HAVE_SOLVDATA)
      writeSolverData(stream);
    if(state >= HAVE_SOLUTION)
      writeSolution(stream);

    fprintf(stream, "%s\n","END OF DATA");
    fclose(stream);
  }
  else
    errorHandler("Unable to create file");

}

void readEntry(char *file, char *tag, void (*readingFunction)() )
{

  FILE *stream;

  stream = fopen(file, "r");

  if(stream != NULL){
    while( positionFilePointer(stream, tag) != CANT_FIND){
      printf("reading: %s\n", tag);
      readingFunction(stream);
    }
    fclose(stream);
  }
}

int positionFilePointer(FILE *stream, char *tag)
{
  char line[80];

  fgets(line, 80, stream);

  while( strstr(line, tag) == NULL && strstr(line, "END") == NULL){
    fgets(line, 80, stream);
  }

  if(strstr(line, "END") != NULL ){
    return CANT_FIND;
  }
  else{
/*
   We have a problem with keywords SEGMENTS and NODES ON SEGMENTS 
   It is not possible to detect which of these tags we have found

   The following quick hack tries to resolve this bug
*/
    if(strstr(line,"NODES ON SEGMENTS") != NULL && 
       strstr(tag, "NODES ON SEGMENTS")==NULL)
      return CANT_FIND;
    else
      return FOUND_IT;
  }
}