// landschaften.h  Berechnung von Landschaft
#pragma once
#ifndef Vec3
#define Vec3 glm::vec3
#endif
//#define MEERWEITE 50.0f  //in main.cpp definiert
const float wurzel3halbe = sqrtf(3.0f)*0.5;

// Koordinatensystem von STL: z-Achse nach oben, x-Achse nach rechts, y-Achse nach hinten 
// Koordinatensystem von OpenGL: y-Achse nach oben, x-Achse nach rechts, z-Achse auf uns zu

struct VecN {
 float a,b;
};

struct VecS { //Vektor fuer Sigmoid-Parameter
 float a,b,k1,k2;
};

struct Wave {
 float alfa, period;
};

struct Teichparameter {
 uint32 meshid,xnmax,znmax;
 float wellenhoehe, wasserstand;
 Vec3 gmin, gmax; //Grenzen des Wassers, TODO:auswerten fuer Test ob Link im Wasser
 std::vector<Wave> waves;
};

class Huegel {
 int numGauss;
 int numSigmoid;
 int numRaster;
 float ymini,ymaxi;
 float y0=0.0f;
 
 //Parameter fuer Teiche:
 //Teichparameter teichpara[8]; //TODO: vorlaeufig beschraenkt auf 8 Gewaesser
 std::vector<Teichparameter> teichpara;
 int numTeiche=0;

 //Parameter fuer Gausshuegel:
 std::vector<Vec3> position; //x-z-Position der Huegelmitte, y=Huegelhoehe
 std::vector<VecN> param; //Parameter: Huegellaenge, Huegelbreite
 // Gauss-Funktion: y(r) = k1*exp(-k2*r^2)  //r ist Abstand von Mitte
 
 //Parameter fuer Sigmoid-Huegel:
 std::vector<Vec3> sposition; //x-z-Position der Huegelmitte, y=Huegelhoehe
 std::vector<VecS> sparam; //Parameter: Huegellaenge, Huegelbreite,k1,k2
 // Sigmoid-Funktion: s(x) = 1.0 / (1+exp(-k*x))  //x kann auch negativ sein
 // ax=xmitte-laenge/2; bx=xmitte+laenge/2;
 // y(x) = (x<xmitte) ? 1/(1+exp(-k1*(x-ax))) : 1/(1+exp(-k1*(bx-x)))
 // az=xmitte-breite/2; bz=xmitte+breite/2;
 // y(z) = (z<zmitte) ? 1/(1+exp(-k2*(z-az))) : 1/(1+exp(-k2*(bz-z)))
 // y(x,z) = y(x)*y(z)

 //Parameter fuer Rasterdatei:
/*
#ifdef TEST_N4
 //#define N100 400
 #define N100 TEST_N4
#else
 #define N100 100
#endif
*/
//int N100=100; //in main.cpp als globale Variable
#define N200 (2*N100)
#define N201 (N200+1)
 //float rasterpunkt[401][401]; //TODO: bisher auf diese Groesse beschraenkt
 //float rasterhfeld[401][401]; //TODO: bisher auf diese Groesse beschraenkt
 BigMatrix2 rasterpunkt; //TODO: initialisieren [N201][N201];
 BigMatrix2 rasterhfeld; //TODO: initialisieren [N201][N201];

public:
 void set_y0(float y) {y0=y;}
 void setNumraster(int n) {numRaster=n;}
 void setRasterpunkt(int i,int j,float ywert) {rasterpunkt.set(i,j,ywert);}
 void setRasterhfeld(int i,int j,float ywert) {rasterhfeld.set(i,j,ywert);}
 float getRast(int i,int j) {return rasterpunkt.get(i,j);}
 float getRasth(int i,int j) {return rasterhfeld.get(i,j);}
 float xmini= -100.0f,xmaxi=100.0f,zmini= -100.0f,zmaxi=100.0f; //Grenzen der Landschaft
 float meerweite=MEERWEITE;
#define MEERGR (100.0f+MEERWEITE)
 float xmini2= -MEERGR, xmaxi2=MEERGR, zmini2= -MEERGR, zmaxi2=MEERGR; //Maximale Grenzen fuer Meer
 void setGrenzen(float x1,float z1,float x2,float z2) {
  xmini=x1; zmini=z1; xmaxi=x2; zmaxi=z2;
  xmini2=x1-meerweite; zmini2=z1-meerweite; xmaxi2=x2+meerweite; zmaxi2=z2+meerweite;
  //printf("setGrenzen() xmini2=%f ... zmaxi2=%f\n",xmini2,zmaxi2);//test
#ifdef _DEBUG
  printf("setGrenzen() xmini=%f xmaxi=%f zmini=%f zmaxi=%f\n",xmini,xmaxi,zmini,zmaxi);//test
  printf("     xmini2=%f xmaxi2=%f zmini2=%f zmaxi2=%f\n",xmini2,xmaxi2,zmini2,zmaxi2);//test
#endif
 }
 void setMeerweite(float meerweit) {meerweite=meerweit; setGrenzen(xmini,zmini,xmaxi,zmaxi);}
 bool ist_in_wasser(Vec3 linkpos,float tiefe) {
  for(int i=0;i<numTeiche-1;i++)
   {
    if(linkpos.x > teichpara[i].gmin.x && linkpos.x < teichpara[i].gmax.x &&
       linkpos.z > teichpara[i].gmin.z && linkpos.z < teichpara[i].gmax.z && linkpos.y<=teichpara[i].gmin.y+tiefe)
     {
      return true;
     }
   }
  return false;
 }
 
 Huegel() {ymini=1e6f; ymaxi= -ymini; numGauss=numSigmoid=0; numRaster=0;}
 void add_gausshuegel(Vec3 pos, VecN para) {
  position.push_back(pos); param.push_back(para);
  numGauss++;
 }
 void add_sigmoidhuegel(Vec3 pos, VecS para) {
  sposition.push_back(pos); sparam.push_back(para);
  numSigmoid++;
 }
 float yf(float x,float z); //Funktion zur Berechnung der Kombination aller Huegel
 Vec3 normalenvektor_berechnen(Vec3 v1,float dx);
 //void minmaxtest();//test
 void teich_minmax(float wellhoehe, Vec3 mitte, float dx, Vec3* vmin, Vec3* vmax);
 void set(uint32 id,uint32 xnma,uint32 znma,float wellho,float wastand,int teichnr,Vec3 vmin,Vec3 vmax) {
   //if(teichnr>=8) {printf("Fehler: zu viele Teiche\n"); teichnr=7;}//test
   if(teichnr >= numTeiche)
    {
     Teichparameter para;
     teichpara.push_back(para);
     numTeiche++;
    }
   if(teichnr >= numTeiche) {printf("Fehler: teichnr=%d falsch\n",teichnr); teichnr=numTeiche-1;}//test
   Teichparameter *p = &teichpara[teichnr];
   p->meshid=id; p->xnmax=xnma; p->znmax=znma; p->wellenhoehe=wellho; p->wasserstand=wastand;
   p->gmin=vmin; p->gmax=vmax;
  }
 void get(uint32 *id,uint32 *xnma,uint32 *znma,float *wellho,float *wastand,int teichnr) {
   if(teichnr>=8) {printf("Fehler: zu viele Teiche\n"); teichnr=7;}//test
   Teichparameter *p= &teichpara[teichnr];
   *id=p->meshid; *xnma=p->xnmax; *znma=p->znmax; *wellho=p->wellenhoehe; *wastand=p->wasserstand;
  }
 void addWave(Wave w,int teichnr) {teichpara[teichnr].waves.push_back(w);}
 float getalfa(uint32 j,float delta,int teichnr) { //TODO
  Teichparameter *p= &teichpara[teichnr];
  p->waves[j].alfa += delta/p->waves[j].period;
  return p->waves[j].alfa;
 }
 void setRaster(int num,const char *dateiname,int smooth);
};

int abrunden(double x)
{
 if(x>=0) return int(x);
 return int(x-0.999);
}

int aufrunden(double x)
{
 if(x>=0) return int(x+0.999);
 return int(x);
}

double sq(double x) {return x*x;}

float Huegel::yf(float x,float z) //Funktion zur Berechnung der Kombination aller Huegel
{
 float y=y0, rquadrat, x1,z1,k1,k2,k3;
 for(int i=0;i<numGauss;i++)
  {
   k1=position[i].y; k2=0.5/param[i].a; k3=0.5/param[i].b;
   x1 = x-position[i].x; z1 = (z-position[i].z)*k3/k2;
   rquadrat = x1*x1+z1*z1;
   y += k1*exp(-k2*rquadrat);
  }
 float yx,yz,ax,bx,az,bz;
 for(int i=0;i<numSigmoid;i++)
  {
   float breitehalbe = sparam[i].a/2;
   float laengehalbe = sparam[i].b/2;
   ax=sposition[i].x-laengehalbe; bx=sposition[i].x+laengehalbe;
   az=sposition[i].z-breitehalbe; bz=sposition[i].z+breitehalbe;
   k1=sparam[i].k1; k2=sparam[i].k2;
   yx = (x<sposition[i].x) ? 1/(1+exp(-k1*(x-ax))) : 1/(1+exp(-k1*(bx-x)));
   yz = (z<sposition[i].z) ? 1/(1+exp(-k2*(z-az))) : 1/(1+exp(-k2*(bz-z)));
   y += yx*yz*sposition[i].y;
  }
 for(int i=0;i<numRaster;i++)
  {
   if(x>xmini && x<xmaxi && z>zmini && z<zmaxi)
    {
     float yplus=0;
     int nn=N100/100;
#ifdef HUEGEL_INTERPOLATION_EINFACHER
     int ix0=int((x-xmini)*nn), jz0=int((z-zmini)*nn);
     int ix1=ix0+1,             jz1=jz0+1;
     //interpolieren zwischen den Rasterpunkten (Variante1):
     if(N100==100)
      {
       Vec3 p1(ix0+xmini, rasterpunkt.get(ix0,jz0), jz0+zmini);
       Vec3 p2(ix1+xmini, rasterpunkt.get(ix1,jz0), jz0+zmini);
       Vec3 p3(ix0+xmini, rasterpunkt.get(ix0,jz1), jz1+zmini);
       Vec3 p4(ix1+xmini, rasterpunkt.get(ix1,jz1), jz1+zmini);
       float ya=(x-p1.x)*(p2.y-p1.y)+p1.y; //Mittelwert zwischen p1 und p2
       float yb=(x-p3.x)*(p4.y-p3.y)+p3.y; //Mittelwert zwischen p3 und p4
       yplus = (z-p1.z)*(yb-ya)+ya; //Mittelwert zwischen ya und yb
      }
     else
      {//TODO
       float xmi=xmini*nn, zmi=zmini*nn;
       Vec3 p1(ix0+xmi, rasterpunkt.get(ix0,jz0), jz0+zmi);
       //Vec3 p2(ix1+xmi, rasterpunkt.get(ix1,jz0), jz0+zmi);
       float p2y=rasterpunkt.get(ix1,jz0);
       //Vec3 p3(ix0+xmi, rasterpunkt.get(ix0,jz1), jz1+zmi);
       float p3x=ix0+xmi, p3y=rasterpunkt.get(ix0,jz1);
       //Vec3 p4(ix1+xmi, rasterpunkt.get(ix1,jz1), jz1+zmi);
       float p4y=rasterpunkt.get(ix1,jz1);
       float ya=(x*nn-p1.x)*(p2y-p1.y)+p1.y; //Mittelwert zwischen p1 und p2
       float yb=(x*nn-p3x)*(p4y-p3y)+p3y; //Mittelwert zwischen p3 und p4
       yplus = (z*nn-p1.z)*(yb-ya)+ya; //Mittelwert zwischen ya und yb
      }
#else
     int ix[2],jz[2]; //4 Rasterpunkte innerhalb derer der Punkt x,z liegt
     //ix[0]=int(x-xmini); jz[0]=int(z-zmini);
     ix[0]=int((x-xmini)*nn); jz[0]=int((z-zmini)*nn);
     ix[1]=ix[0]+1;      jz[1]=jz[0]+1;
     double abst[4]; //Abstand zu jedem der 4 Rasterpunkte
     abst[0]=sqrt(sq(x-xmini-ix[0])+sq(z-zmini-jz[0]));
     abst[1]=sqrt(sq(x-xmini-ix[1])+sq(z-zmini-jz[0]));
     abst[2]=sqrt(sq(x-xmini-ix[0])+sq(z-zmini-jz[1]));
     abst[3]=sqrt(sq(x-xmini-ix[1])+sq(z-zmini-jz[1]));
     //wenn Punkt fast exakt auf Rasterpunkt liegt, dann diesen verwenden:
     const float dd=0.001;
     if(abst[0]<dd) yplus=rasterpunkt.get(ix[0],jz[0]);
     else if(abst[1]<dd) yplus=rasterpunkt.get(ix[1],jz[0]);
     else if(abst[2]<dd) yplus=rasterpunkt.get(ix[0],jz[1]);
     else if(abst[3]<dd) yplus=rasterpunkt.get(ix[1],jz[1]);
     else if(mittelwert_methode==0) //TODO
      {//sonst interpolieren zwischen den Rasterpunkten (Variante0):
       Vec3 p1(ix[0]+xmini, rasterpunkt.get(ix[0],jz[0]), jz[0]+zmini);
       Vec3 p2(ix[1]+xmini, rasterpunkt.get(ix[1],jz[0]), jz[0]+zmini);
       Vec3 p3(ix[0]+xmini, rasterpunkt.get(ix[0],jz[1]), jz[1]+zmini);
       Vec3 p4(ix[1]+xmini, rasterpunkt.get(ix[1],jz[1]), jz[1]+zmini);
       float ya=0;
       //ya=(p1.y+p2.y+p3.y+p4.y)/4; //Mittelwert zwischen allen 4 Punkten
       //if(ya!=0) printf("ya = (%f+%f+%f+%f)/4 = %f\n",p1.y,p2.y,p3.y,p4.y,ya);//test
       if(abst[0]<abst[1] && abst[0]<abst[2] && abst[0]<abst[3]) ya=p1.y; //naechstgelegener Punkt verwenden
       else if(abst[1]<abst[2] && abst[1]<abst[3]) ya=p2.y;
       else if(abst[2]<abst[3]) ya=p3.y;
       else ya=p4.y;
       yplus = ya;
      }
     else if(mittelwert_methode==1)
      {//sonst interpolieren zwischen den Rasterpunkten (Variante1):
       Vec3 p1(ix[0]+xmini, rasterpunkt.get(ix[0],jz[0]), jz[0]+zmini);
       Vec3 p2(ix[1]+xmini, rasterpunkt.get(ix[1],jz[0]), jz[0]+zmini);
       Vec3 p3(ix[0]+xmini, rasterpunkt.get(ix[0],jz[1]), jz[1]+zmini);
       Vec3 p4(ix[1]+xmini, rasterpunkt.get(ix[1],jz[1]), jz[1]+zmini);
       float ya=(x-p1.x)*(p2.y-p1.y)+p1.y; //Mittelwert zwischen p1 und p2
       float yb=(x-p3.x)*(p4.y-p3.y)+p3.y; //Mittelwert zwischen p3 und p4
       yplus = (z-p1.z)*(yb-ya)+ya; //Mittelwert zwischen ya und yb
      }
     else if(mittelwert_methode==2)
      {//sonst interpolieren zwischen den Rasterpunkten (Variante2):
       double gewichtsum=0.0, eins=1.0;
       for(int j=0;j<4;j++) gewichtsum += eins/abst[j];
       yplus = rasterpunkt.get(ix[0],jz[0])/(gewichtsum*abst[0]);
       yplus += rasterpunkt.get(ix[1],jz[0])/(gewichtsum*abst[1]);
       yplus += rasterpunkt.get(ix[0],jz[1])/(gewichtsum*abst[2]);
       yplus += rasterpunkt.get(ix[1],jz[1])/(gewichtsum*abst[3]);
      }
     else if(mittelwert_methode>=3) //TODO: eventuell noch andere Mittelwertmethode
      {//testvariante:
       yplus = rasterpunkt.get(ix[0],jz[0]);
      }
#endif //HUEGEL_INTERPOLATION_EINFACHER
     y += yplus;
    }
   if(numRaster>1)
    {printf("TODO: bisher auf eine einzige Rasterdatei beschraenkt\n"); numRaster=1;}//test
  }
 if(y<ymini) ymini=y;
 if(y>ymaxi) ymaxi=y;
 return y;
}

Vec3 Huegel::normalenvektor_berechnen(Vec3 v1,float dx)
{
 float dd=dx*0.01;
 Vec3 v2,v3,nv1,nv2;
 v2.x=v1.x; v2.z=v1.z+dd;    v2.y=yf(v2.x,v2.z);
 v3.x=v1.x+dd; v3.z=v1.z+dd; v3.y=yf(v3.x,v3.z);
 nv1 = v2 - v1;
 nv2 = v3 - v2;
 return glm::normalize(glm::cross(nv1,nv2));
}

Huegel landschaft1;

void landschaft1_erstellen(int nhueg,float *ktab,float *postabzx,
			   int nsighueg,float *sigmktab,float *sigmpostabzx)
{
#ifdef _DEBUG
 printf("test: landschaft1_erstellen(nhueg=%d,ktab,postabzx, nsighuegel=%d,..)\n",nhueg,nsighueg);//test
#endif
 Vec3 position; //x-z-Position der Huegelmitte, y=Huegelhoehe
 VecN param; //Parameter: Huegellaenge, Huegelbreite
 VecS sparam; //Parameter fuer Sigmoid-Huegel
 for(int i=0;i<nhueg;i++)
  {
   int j=i*2, k=i*3;
   position.y = ktab[k];
   param.a = ktab[k+1]; param.b = ktab[k+2];
   position.x = postabzx[j+1];  position.z = postabzx[j];
   landschaft1.add_gausshuegel(position,param);
  }
 for(int i=0;i<nsighueg;i++)
  {
   int j=i*2, k=i*5;
   position.y = sigmktab[k];
   sparam.a = sigmktab[k+1]; sparam.b = sigmktab[k+2];
   sparam.k1 = sigmktab[k+3]; sparam.k2 = sigmktab[k+4];
   position.x = sigmpostabzx[j+1];  position.z = sigmpostabzx[j];
   landschaft1.add_sigmoidhuegel(position,sparam);
  }
 landschaft1.setRaster(1,rasterdateiname,rastersmooth);
}

void Huegel::setRaster(int num,const char *dateiname,int smooth)
{
 FILE *fp=fopen(dateiname,"r");
 if(fp==NULL) {fprintf(stderr,"Error: \"%s\" nicht gefunden\n",dateiname); return;}
 char zeile[200];
 int xi0=0,xi1=0,yi0=0,yi1=0;
 float maxhoehe=0;
 float dx=1.0, dy=1.0;
 while(getline(fp,zeile,200))
  {
   if(strncmp(zeile,"xminxmax:",9)==0) {sscanf(&zeile[9],"%d %d",&xi0,&xi1);}
   else if(strncmp(zeile,"yminymax:",9)==0) {sscanf(&zeile[9],"%d %d",&yi0,&yi1);}
   else if(strncmp(zeile,"dxdy:",5)==0) {sscanf(&zeile[5],"%f %f",&dx,&dy);}
   else if(strncmp(zeile,"maxhoehe:",9)==0) {sscanf(&zeile[9],"%f",&maxhoehe); break;}
   else {printf("Fehler: fehlende Paramter in \"%s\"\n",dateiname); return;}
  }
 if(maxhoehe==0) {printf("Fehler: fehlerhafte Paramter in \"%s\"\n",dateiname); return;}
 if(xi0!= -100 || xi1!=100 || yi0!= -100 || yi1!=100) //TODO
  {printf("Fehler: xy-Werte in Rasterdatei muss immer im Bereich -100 bis 100 sein.\n"); return;}//test
 float ywert;
 int np=0;//test
 int jmax=(int)((yi1-yi0)/dy+0.5);
 int imax=(int)((xi1-xi0)/dx+0.5);
 rasterpunkt.init(N201,N201);
 rasterhfeld.init(N201,N201);
 for(int j=0; j<=jmax; j++)
  for(int i=0; i<=imax; i++)
   {
    int n=fscanf(fp,"%f",&ywert); //Hoehe lesen
    if(n!=1) {printf("Fehler: zu wenig Daten in \"%s\"\n",dateiname); return;}
    setRasterpunkt(i,j,ywert);
    np++;//test
   }
 fclose(fp);
 //printf("Huegel::setraster()  %d Rasterpunkte gespeichert\n",np);//test
 if(smooth>=1)
  {
   //printf("Testpunkt: Huegel::setraster() smooth=%d\n",smooth);//test
   int nmal=1;
   if(smooth>3) nmal=smooth-2;
   if(smooth>3) smooth=3;//TODO: glaetten ueber mehr als 3 Punkte
   for(;nmal>0;--nmal)
    {
     int i0=(smooth==2)?0:1, j0=(smooth>=2)?1:0;
     //int i2=(xi1-xi0)-i0,    j2=(yi1-yi0)-j0;
     int i2=imax-i0,    j2=jmax-j0;
     float y;
     int n1=(smooth==3)?4:2;
     for(int j=j0; j<=j2; j++)
      for(int i=i0; i<=i2; i++)
       {
	int n=n1;
	y=getRast(i,j)*n;
	if(smooth&1) {y += getRast(i-1,j); y += getRast(i+1,j); n+=2;}
	if(smooth&2) {y += getRast(i,j-1); y += getRast(i,j+1); n+=2;}
	setRasterhfeld(i,j,y/n);
       }
     for(int j=0; j<=jmax; j++)
      for(int i=0; i<=imax; i++)
       {
	y=getRasth(i,j);
	setRasterpunkt(i,j,y);
       }
    }
  }
 setNumraster(1);
}

/*
void Huegel::minmaxtest() //test
{
 float x,y,z;
 float ymi=1e6f, yma= -ymi, xmi=0, xma=0, zmi=0, zma=0;
 int ix,iz;
 float dx=0.125f, dz=dx*wurzel3halbe;
 int xnmax= (int)((xmaxi-xmini)/dx+1.5);
 int znmax= (int)((zmaxi-zmini)/dz+1.5);
 printf("xmini=%f xmaxi=%f\n",xmini,xmaxi);//test
 printf("zmini=%f zmaxi=%f\n",zmini,zmaxi);//test
 for(ix=0,x=xmini; ix<=xnmax; ix++,x+=dx)
  for(iz=0,z=zmini; iz<=znmax;iz++,z+=dz)
   {
    y=yf(x,z);
    if(y<ymi) {ymi=y; xmi=x; zmi=z;}
    if(y>yma) {yma=y; xma=x; zma=z;}
   }
 if(ymini!=ymi || ymaxi!=yma) printf("ymini = %f  ymaxi = %f\n",ymini,ymaxi);
 printf("Tiefster  Punkt: x: %f y: %f z: %f\n",xmi,ymi,zmi);
 printf("Hoechster Punkt: x: %f y: %f z: %f\n",xma,yma,zma);
}
*/

void Huegel::teich_minmax(float wellenhoehe, Vec3 mitte, float dx, Vec3* vmin, Vec3* vmax)
{
 float dz=dx*wurzel3halbe;
 float wasserstand=mitte.y+wellenhoehe;
 float x1=0,z1=0,x2=0,z2=0;
 float y1,y2; y1=y2=wasserstand;
 Vec3 v;
#ifdef _DEBUG
 printf("Huegel::teich_minmax() wasserstand=%f dx=%f dz=%f\n",wasserstand,dx,dz);//test
#endif
 for(y1=y2=wasserstand; y1<=wasserstand || y2<=wasserstand || x1<=xmini2 || x2>=xmaxi2;)
  {
   v=mitte; v.x += x1;
   y1=yf(v.x,v.z); if(y1<=wasserstand) x1-=dx;
   v=mitte; v.x += x2;
   y2=yf(v.x,v.z); if(y2<=wasserstand) x2+=dx;
  }
 for(y1=y2=wasserstand; y1<=wasserstand || y2<=wasserstand || z1<=zmini2 || z2>=zmaxi2;)
  {
   v=mitte; v.z += z1;
   y1=yf(v.x,v.z); if(y1<=wasserstand) z1-=dz;
   v=mitte; v.z += z2;
   y2=yf(v.x,v.z); if(y2<=wasserstand) z2+=dz;
  }
 
 //testen ob durch x1,z1,x2,z2 definiertes Viereck vollstaendig an Land liegt:
#ifdef _DEBUG
 printf("Grenzen pruefen: %f %f %f %f\n",x1,z1,x2,z2);//test
#endif
 for(int flag=4,timeout=0; flag!=0 && ++timeout<1000000;)
  {
   float x,y,z;
   float ymin=1e6f, xmin=0.0f, zmin=0.0f;
   flag=0;
   for(x=x1,z=z1; x<=x2; x+=dx)
    {
     if((y=yf(mitte.x+x,mitte.z+z)) <= wasserstand)
      {
       flag=1; if(y<ymin) {ymin=y; xmin=x;}
      }
    }
   if(flag==1)
    {
     while(z1>zmini2 && yf(mitte.x+xmin,mitte.z+z1) <= wasserstand) z1 -= dz;
    }
   ymin=1e6f; 
   for(x=x1,z=z2; x<=x2; x+=dx)
    {
     if((y=yf(mitte.x+x,mitte.z+z)) <= wasserstand)
      {
       flag=2; if(y<ymin) {ymin=y; xmin=x;}
      }
    }
   if(flag==2)
    {
     while(z2<zmaxi2 && yf(mitte.x+xmin,mitte.z+z2) <= wasserstand) z2 += dz;
    }
   ymin=1e6f; 
   for(z=z1,x=x1; z<=z2; z+=dz)
    {
     if((y=yf(mitte.x+x,mitte.z+z)) <= wasserstand)
      {
       flag=3;
       if(y<ymin) {ymin=y; zmin=z;}
      }
    }
   if(flag==3)
    {
     while(yf(mitte.x+x1,mitte.z+zmin) <= wasserstand) x1 -= dx;
    }
   ymin=1e6f; 
   for(z=z1,x=x2; z<=z2; z+=dz)
    {
     if((y=yf(mitte.x+x,mitte.z+z)) <= wasserstand)
      {
       flag=4;
       if(y<ymin) {ymin=y; zmin=z;}
      }
    }
   if(flag==4)
    {
     while(x2<xmaxi2 && yf(mitte.x+x2,mitte.z+zmin) <= wasserstand) x2 += dx;
    }
  }
#ifdef _DEBUG
 printf("neue Grenzen: %f %f %f %f\n",x1,z1,x2,z2);//test
#endif

 vmin->x = x1; vmin->z = z1;
 vmax->x = x2; vmax->z = z2;
}
