/* gpxread.h   letzte Aenderung: 21.7.2022 */

//#include <stdio.h>
//#include <stdlib.h>
//#include <string.h>
//#include <math.h>

//const float PI = atanf(1.0f)*4;
//const float GRAD = 2*PI/360;

/************************ Globale Variablen ***************************/
const int MINXY= -100, MAXXY=100; //TODO: sollte von User setzbar sein
//const float MAXZ=50.0; //von User setzbar mit flag -a (hmax)
//#define MAXB (MAXXY-MINXY)
#define MAXB ((MAXXY-MINXY)*N100/100)
const float FMIN= -1e20f, FMAX=1e20f; //minimale und maximale Flisszahl bei minmax-Suche
float Gxmin=0, Gxmax=0, Gymin=0, Gymax=0;
bool swissflag=false;
//int bflag=0;//TODO schon in main.cpp

std::vector<glm::vec3> gpx_wegpunkte; //Wegpunkte um km-Faehnchen zu setzen
int iwpkt=0; //Anzahl Eintraege in gpx_wegpunkte

glm::vec3 gpx_gipfelpos={0,0,0};//hoechster Punkt von einem GPX-Track, oder auch xyz-Datei
glm::vec3 gpx_tiefpunkt={0,0,0};//tiefster Punkt
glm::vec3 gpx_startpos={0,0,0};//Start von GPX-Track
glm::vec3 gpx_zielpos={0,0,0};//Ziel von GPX-Track

float globalScale=1.0;//TODO

void denkmal_scalefaktor_setzen(float x) //x=Gxmax-Gxmin, Denkmale Tuerme ... Skalierfaktor anpassen
{
 //printf("denkmal_scalefaktor_setzen(x=%f)\n",x);//test
 if(x<=0.05) globalScale=1.0; //TODO Schwellwerte anpassen
 else if(x<=0.10) globalScale=0.5;
 else if(x<=0.20) globalScale=0.25;
 else if(x<=0.40) globalScale=0.15;
 else globalScale=0.1;//TODO: staerkste Verkleinerung bei grosser Karte
}

/********************** Koordinaten umrechnen *************************/
class Koordinatenumrechnung
{
 public:
 float xmitte,ymitte,znullpunkt,hscal;
 //float xscalfaktor,yscalfaktor;//schon in main.cpp als globale Variablen definiert
 void init(float xm,float ym,float z0,float zscal) {xmitte=xm; ymitte=ym; znullpunkt=z0; hscal=zscal;}
 float xuser2openglx(float x) {return (x-xmitte)*xscalfaktor;}
 float yuser2openglz(float y) {return (y-ymitte)*yscalfaktor;}
 float zuser2opengly(float z) {return (z-znullpunkt)*hscal;}
 float xopengl2userx(float x) {return x/xscalfaktor+xmitte;}
 float zopengl2usery(float z) {return z/yscalfaktor+ymitte;}
 float yopengl2userz(float y) {return (hscal==0 ? 0 : y/hscal) + znullpunkt;}
};
Koordinatenumrechnung koord;

//class Vector //TODO: direkt glm::vec3 benutzen besser?
class Vector : public glm::vec3 //TODO: direkt glm::vec3 benutzen besser?
{
public:
 //float x=0,y=0,z=0; //schon in glm::vec3
 void set(float a,float b,float c) {x=a; y=b; z=c;}
 void set(glm::vec3 v) {x=v.x; y=v.y; z=v.z;}
 void scal()
  {x = koord.xuser2openglx(x);
   float zh=z;
   z = koord.yuser2openglz(y);
   y = koord.zuser2opengly(zh);
  }
};

void vec3scal(glm::vec3& v)
{
 v.x = koord.xuser2openglx(v.x);
 float z=v.z;
 v.z = koord.yuser2openglz(v.y);
 v.y = koord.zuser2opengly(z);
}

double distanz(double lat0,double lon0,double lat1,double lon1)
{
 //Distanz zwischen zwei GPS-Koordinaten:
 double latb=lat1*GRAD,lata=lat0*GRAD,lonb=lon1*GRAD,lona=lon0*GRAD;
 const double r=6.38e6; //Erdradius in Meter
 return acos(sin(latb)*sin(lata)+cos(latb)*cos(lata)*cos(lonb-lona))*r;
}

double meter_nach_gpx_umrechnen(double meter)
{
 const double r=6.38e6; //Erdradius in Meter
 return 360.0*meter/(ZWEIPI*r);
}

/************************* Vordeklarationen ***************************/
//von gpxread.cc uebernommen und angepasst:
int gpxtrack_in_xyz_umwandeln(FILE *fp1,FILE *fp2);
void rasterdatei_erstellen(int npunkte,const char *quelle,const char *ziel,float minaquadrat,bool geoflag);
float naechster_punkt_suchen(int ix,int iy,Vector punkte[],int npunkte,float minaquadrat);
float swissalti3d(double x,double y);

//vom Hauptprogramm gebrauchte Funktionen:
void gpxread(const char *ziel,const char *quelle); //ziel=Rasterdaei, quelle=GPX-Datei
void rasterdatei_erstellen(const char *ziel,const char *quelle); //aus xyz-Datei erstellen

//nur in gpxread.h gebrauchte Funktionen:
int karten_zusammensetzen(int j,int nr,char *kartenname); //TODO Rueckgabewwert=nr, j muss von 0 bis 10 laufen

/************** Klassen und Hauptteil des Programms *******************/
Vec3 linien_color_tabelle[8] = {
  {0.0, 0.0, 1.0}, //[0]=blau
  {1.0, 0.0, 1.0}, //violet
  {1.0, 1.0, 0.0}, //gelb
  {0.0, 1.0, 1.0}, //blaugruen
  {0.5, 0.5, 1.0}, //hellblau
  {1.0, 0.5, 1.0}, //hellviolet
  {0.6, 1.0, 0.6}, //hellgruen
  {0.6, 0.6, 0.0} //[7]=dunkelgelb, braun
 };//TODO: ev. groessere Farbtabelle
 
class Linienzug
{
public:
 uint32 numPunkte;
 uint32 typ;//0=Werte in GPX-Koordinaten, 1=OpenGL-Koordinaten, ... //TODO
 std::vector<Vec3> punkte;
 double strecke=0;
 Material material;
 void setcolor(float r,float g,float b) {setcolor(Vec3{r,g,b});} //TODO
 void setcolor(Vec3 v) //TODO
  {//Vec3 v(r,g,b);
   material.ambient=material.diffuse=material.specular=v;
   material.emissive={0,0,0};
   material.shininess=4.0;
  }
 void setcolor(int i) {setcolor(linien_color_tabelle[i&7]);}//TODO: ev. groessere Farbtabelle
 Linienzug() {numPunkte=0; typ=0; setcolor(0,0,1); strecke=0;}
 void init() {strecke=0; typ=0; if(numPunkte!=0) {punkte.clear(); numPunkte=0;}}
 void add(Vec3 punkt) {punkte.push_back(punkt); numPunkte++;}
 void add(float x,float y,float z) {add(Vec3{x,y,z});}
};
std::vector<Linienzug> linienzuege; //mit linienzuege.size() erhaelt man Anzahl

float ele_read(const char *s)
{
 float z;
 sscanf(s,"%f",&z);
 return z;
}

int latlon_read(const char *s,double *px,double *py)
{
 int n=0;
 while(n<2 && *s!=0)
  {
   if(strncmp(s,"lat=\"",5)==0)      {sscanf(&s[5],"%lf",py); n++;}
   else if(strncmp(s,"lon=\"",5)==0) {sscanf(&s[5],"%lf",px); n++;}
   while(*s!=' ' && *s!=0) {s++;}
   while(*s==' ') {s++;}
  }
 if(n!=2) fprintf(stderr,"Fehler in latlon_read(): keine 2 Werte gefunden\n");//test
 return n;
}

void gpxread(const char *ziel,const char *quelle)
{
 const char *tmpname="tmpgpx.xyz";
 FILE *fp1=myfopen(quelle,"r");
 FILE *fp2=myfopen(tmpname,"w"); //TODO: eigentlich unnoetige Zwischendatei
 if(!fp1 || !fp2) return;
 xscalfaktor=yscalfaktor=0.0; //Skalierfaktoren bei GPX immer automatisch ausrechnen
 hscal = (gpxhscal!=0) ? gpxhscal : 0.0;
 //hscal = (gpxhscal!=0) ? gpxhscal : 1.0;//test
 logflags=0; //bei GPX niemals logarithmische Skalen
 int npkt=gpxtrack_in_xyz_umwandeln(fp1,fp2);
 fclose(fp1); fclose(fp2);
#ifdef _DEBUG
 printf("%s eingelesen: %d Punkte\n",quelle,npkt);//test
#endif
 if(tracksmooth>1)
  {
   printf("Testpunkt: tracksmooth=%d\n",tracksmooth);//test
   Vector *v1=new Vector[npkt];
   fp1=myfopen(tmpname,"r"); if(!fp1) return;
   for(int i=0;i<npkt;i++)
    fscanf(fp1,"%f %f %f",&v1[i].x,&v1[i].y,&v1[i].z);
   fclose(fp1);
   Vector *v2=new Vector[npkt];
   int k1 = -(tracksmooth-1)/2, k2=k1+tracksmooth;
   int nz=linienzuege.size(); //Anzahl Linienzuege (Tracks)
   int i1=0;//index fuer v1[]
   for(int nzi=0;nzi<nz;nzi++)
    {
     int np=linienzuege[nzi].numPunkte;
     for(int i=0;i<np;i++,i1++)
      {
       int n=1;
       //float summe=v1[i1].z; //nur Hoehe mitteln
       float summex=v1[i1].x, summey=v1[i1].y, summez=v1[i1].z; //alles mitteln
       for(int k=k1+i,i2=k1+i1; k<k2+i; k++,i2++)
	{
	 if(k>=0 && k<np)
	  //{summe += v1[i2].z; n++;}//nur Hoehe mitteln
	 {summex += v1[i2].x; summey += v1[i2].y; summez += v1[i2].z; n++;}
	}
       //v2[i1].x = v1[i1].x;  v2[i1].y = v1[i1].y;  v2[i1].z = summe/n;// nur Hoehe mittln
       v2[i1].set(summex/n, summey/n, summez/n); //alles mitteln
      }
    }
   
   fp2=myfopen(tmpname,"w"); if(!fp2) return;
   for(int i=0;i<npkt;i++)
    fprintf(fp2,"%f %f %f\n", v2[i].x, v2[i].y, v2[i].z);
   fclose(fp2);
   delete[] v1; delete[] v2;
  }
 float minaquadrat = (gpxtrackwidth==0) ? 16 : gpxtrackwidth*gpxtrackwidth; //TODO: Defaultwert?
 printf("1.Aufruf von rasterdatei_erstellen()\n");//test
 rasterdatei_erstellen(npkt,tmpname,ziel,minaquadrat,true);
 unlink(tmpname);
}

int gpxtrack_in_xyz_umwandeln(FILE *fp1,FILE *fp2)
{
 char zeile[200], *str;
 int nr=0;
 bool trkpt_flag=false;
 //float x=0,y=0,z=0;//falsch
 double x=0,y=0,z=0;
 double lat0=0,lon0=0,lat1=0,lon1=0, dist=0,dist1km=1000;//fuer Distanz-Berechnungen
 double aufstieg=0, zvorher=0;//TODO
 bool firstpoint=true;
 Linienzug linienzug;
 while(getline(fp1,zeile,200))
  {
   for(str=zeile;*str==' ';str++) {}//Leerstellen ueberlesen
   if(trkpt_flag)
    {
     if(strncmp(str,"</trkpt>",8)==0) {trkpt_flag=false;}
     else if(strncmp(str,"<ele>",5)==0)
      {
       z = ele_read(&str[5]);
       if(bflag) {wgs2swiss(x,y);}
       fprintf(fp2,"%f %f %f\n",x,y,z);
       linienzug.add(x,y,z); //TODO: anpassen nach allfaelliger Glaettung (smooth)
       //Aufstieg-Berechnung:
       if(zvorher!=0)
	{
	 if(z>zvorher) aufstieg += z-zvorher;//TODO
	}
       zvorher=z;
       //Distanz-Berechnung:
       lat0=lat1; lon0=lon1;  lat1=y; lon1=x;
       if(firstpoint) {firstpoint=false;}
       else
	{
	 double ds;
	 if(bflag) ds=sqrt(sq(lat1-lat0)+sq(lon1-lon0));
	 else      ds=distanz(lat0,lon0,lat1,lon1);
	 if(isnan(ds)) printf("Error: ds=%f lat0=%f lon0=%f lat1=%f lon1=%f bflag=%d\n",
			      ds,lat0,lon0,lat1,lon1,bflag);//test
	 else {dist += ds; linienzug.strecke += ds;}
	 if(dist>=dist1km)
	  {
	   glm::vec3 v(x,y,z);
	   gpx_wegpunkte.push_back(v); iwpkt++;
	   dist1km += 1000; //alle 1km ein Faehnchen
	  }
	}
      }
     else if(strncmp(str,"<time>",6)==0)
      {
       //TODO: Zeit auslesen fuer Geschwindigkeitsberechnung
      }
     else if(strncmp(str,"<extensions>",12)==0) {}
     else if(strncmp(str,"<ns3:",5)==0) {}
     else if(strncmp(str,"<gpxtpx:",8)==0) {}
     else if(str[1]!='/') printf("unbekannter Eintrag in <trkpt>: \"%s\"\n",str);
    }
   else if(strncmp(str,"<trkpt ",7)==0)
    {
     nr++; //Anzahl Punkte zaehlen
     trkpt_flag=true;
     int n=latlon_read(&str[7],&x,&y);
     if(n!=2) fprintf(stderr,"Fehler: %s\n",str);
    }
   else if(strncmp(str,"<trkseg>",8)==0)
    {
     //printf("<trkseg> erkannt\n");//test
     linienzug.init();
    }
   else if(strncmp(str,"</trkseg>",9)==0)
    {
     //printf("</trkseg> erkannt\n");//test
     int n=linienzuege.size();
     linienzug.setcolor(n); //unterschiedliche Farben fuer verschiedene Tracks
     linienzuege.push_back(linienzug);
    }
  }
 if(iwpkt>=25)
  {//zu viele Faehnchen entfernen:
   //printf("test: iwpkt=%d\n",iwpkt);//test
   int i=0,j=0;
   if(iwpkt>=1250)
    {
     for(i=99;i<iwpkt;i+=100,j++) {gpx_wegpunkte[j] = gpx_wegpunkte[i];}
    }
   else if(iwpkt>=250)
    {
     for(i=49;i<iwpkt;i+=50,j++) {gpx_wegpunkte[j] = gpx_wegpunkte[i];}
    }
   else if(iwpkt>=50)
    {
     for(i=9;i<iwpkt;i+=10,j++) {gpx_wegpunkte[j] = gpx_wegpunkte[i];}
    }
   else
    {
     for(i=4;i<iwpkt;i+=5,j++) {gpx_wegpunkte[j] = gpx_wegpunkte[i];}
    }
   iwpkt = (j>25)? 25 : j;
  }
 if(startpos_auto==1) startpos_auto=2;
 printf("GPX-Datei eingelesen: dist=%f  iwpkt=%d\n",dist,iwpkt);//test
 if(linienzuege.size()!=1) printf(" %d GPX-Tracks\n",(int)linienzuege.size());//test
 for(uint32 i=0;i<linienzuege.size();i++)//test
  {
   printf(" %d.Track: %.3f km\n",i+1,linienzuege[i].strecke*0.001);//test
  }
 printf(" Gesamtstrecke: %.3f km\n",dist*0.001);//test
 printf(" Gesamtaufstieg: ca. %d Meter\n",int(aufstieg));//TODO
 return nr;
}

int minmax(float x,float& xmin, float& xmax)
{
 int j=0;
 if(x<xmin) {xmin=x; j++;}
 if(x>xmax) {xmax=x; j+=2;}
 return j;
}

void rasterdatei_erstellen(const char *ziel,const char *quelle)
{
 FILE *fp=myfopen(quelle,"r"); if(!fp) return;
 int c, npkt=0;
 if(hatendung(quelle,".csv"))
  {
#ifdef _DEBUG
   printf("testpunkt: csv-Datei einlesen\n");//test
#endif
   do {c=getc(fp);} while(c==' '); //erstes Zeichen lesen
   ungetc(c,fp);
   if(!isdigit(c) && c!='-' && c!='.') //wenn erste Zeile nicht mit einer Zahl beginnt
    {
     char zeile[400], *s;
     getline(fp,zeile,400); //erste Zeile einlesen
#ifdef _DEBUG
     printf("erste Zeile von csv-Datei: \"%s\"\n",zeile);//test
#endif
     for(s=zeile; *s!=0;)
      {
       while(*s==' ') s++;
       if(strncmp(s,"xscal:",6)==0) sscanf(&s[6],"%f",&xscalfaktor);
       else if(strncmp(s,"yscal:",6)==0) sscanf(&s[6],"%f",&yscalfaktor);
       else if(strncmp(s,"hscal:",6)==0) sscanf(&s[6],"%f",&hscal);
       else if(strncmp(s,"hoffset:",8)==0) sscanf(&s[8],"%f",&zoffset);
       else if(strncmp(s,"logflags:",9)==0) sscanf(&s[9],"%d",&logflags);
       while(*s!=0 && *s!=',') {s++;}
       if(*s==',') s++;
      }
#ifdef _DEBUG
     printf("logflags=%d\n",logflags);//test
#endif
    }
   double x=0,y=0,z=0;
   double xmi=1e20,xma= -1e20, ymi=1e20,yma= -1e20, zmi=1e20,zma= -1e20;
   int i;
   for(i=0;;i++)
    {
     int n=fscanf(fp,"%lf%*c%lf%*c%lf",&x,&y,&z);
     if(n!=3) break;
#ifdef _DEBUG
     if(i<3) printf("%d: x=%g y=%g z=%g\n",i,x,y,z);//test
#endif
     if(x<xmi) {xmi=x;} if(x>xma) {xma=x;}
     if(y<ymi) {ymi=y;} if(y>yma) {yma=y;}
     if(z<zmi) {zmi=z;} if(z>zma) {zma=z;}
    }
   npkt=i;
#ifdef _DEBUG
   printf("letzter Punkt: i=%d x=%g y=%g z=%g\n",i-1,x,y,z);//test
   printf("%d Punkte gelesen\n",npkt);//test
   printf("Extremwerte: x: %g %g  y: %g %g  z:%g %g\n",xmi,xma,ymi,yma,zmi,zma);//test
#endif
  }
 else //kein .csv
  {
   while((c=getc(fp))!=EOF) //Anzahl Zeilen zaehlen == Anzahl xyz-Punkte
    {if(c=='\n') npkt++;}
  }
 fclose(fp);
 float minaquadrat = (xyztrackwidth==0) ? 16 : xyztrackwidth*xyztrackwidth; //TODO: Defaultwert?
 rasterdatei_erstellen(npkt,quelle,ziel,minaquadrat,false);
}

void rasterdatei_erstellen(int npunkte,const char *quelle,const char *ziel,float minaquadrat,bool geoflag)
{
//#ifdef _DEBUG
 printf("rasterdatei_erstellen(%d, %s, %s, %f)\n",npunkte,quelle,ziel,minaquadrat);//test
//#endif
 FILE *fp=myfopen(quelle,"r"); if(!fp) return;//Fehlermeldung schon in myfopen()
 float x=0,y=0,z=0;
 if(npunkte==1) npunkte++;//TODO: es muessen mindestens 2 Punkte sein
 Vector *punkte=new Vector[npunkte];
 float xmin=1e20,xmax= -1e20, ymin=1e20,ymax= -1e20, zmin=1e20,zmax= -1e20;
 char zeile[200];
 int i,j;
 Vector tiefsterpunkt,hoechsterpunkt;
 int izmin=0,izmax=0;
 //float rasterhoehen[MAXB+1][MAXB+1]; //anders machen wenn MAXB eine Variable ist
 BigMatrix2 rasterhoehen(MAXB+1,MAXB+1); //so ist Zugriff anders:
 // x=rasterhoehen[i][j] wird zu x=rasterhoehen.get[i,j]
 // rasterhoehen[i][j]=x wird zu rasterhoehen.set[i,j,x]
 //printf("BigMatrix2 rasterhoehen(MAXB+1,MAXB+1)  fertig initialisiert\n");//test

 bool rasterhoehen_gesetzt=false;
 if(hatendung(quelle,".csv"))
  {
   int c=' ';
   while(c==' ' || c=='\t') c=getc(fp);//Leerzeichen ueberlesen
   if(!isdigit(c) && c!='-' && c!='.')
    {while((c=getc(fp))!=EOF && c!='\n') {}}//Kommentarzeile ueberlesen
   else ungetc(c,fp);
#ifdef _DEBUG
   float xvorher=0, yvorher=0;//test
   float dxmin=1e20,dxmax= -1e20,dymin= 1e20,dymax= -1e20;//test
#endif
   for(i=0; c!=EOF;)
    {
     int n=fscanf(fp,"%f%*c%f%*c%f",&x,&y,&z);
     if(n==EOF) break;
     if(n!=3) {printf("Fehler beim Einlesen von \"%s\": n=%d\n",quelle,n);}
     if(logflags&1) x = log10f(x);
     if(logflags&2) y = log10f(y);
     if(logflags&4) z = log10f(z);
#ifdef _DEBUG
     if(i!=0 && x!=xvorher)//test
      {
       if(x-xvorher < dxmin) dxmin = x-xvorher;
       if(x-xvorher > dxmax) dxmax = x-xvorher;
      }
     if(i!=0 && y>yvorher)//test
      {
       if(y-yvorher < dymin) dymin = y-yvorher;
       if(y-yvorher > dymax) dymax = y-yvorher;
      }
     xvorher=x; yvorher=y;//test
#endif
     punkte[i].set(x,y,z);
     j=minmax(x,xmin,xmax);
     j=minmax(y,ymin,ymax);
     j=minmax(z,zmin,zmax);
     if(j&1) {izmin=i;}
     if(j&2) {izmax=i;}
     i++;
     c=getc(fp); //Komma oder Zeilentrenner ueberlesen
    }
#ifdef _DEBUG
   printf("%d Punkte von csv-Datei eingelesen\n",i);//test
   printf(" dx: %f ... %f\n",dxmin,dxmax);//test
   printf(" dy: %g ... %g\n",dymin,dymax);//test
#endif
  }
 else //if(hatendung(quelle,".gpx") || hatendung(quelle,".xyz")) //kein .csv
  {// bei .gpx wurde schon in .xyz umgerechnet
   for(i=0; getline(fp,zeile,200); i++)
    {
     const char *s;
     for(s=zeile; *s==' ' || *s=='\t' || (*s&0xFF)>=0x80; s++) {}
     if(*s==0) continue;
     int n=sscanf(s,"%f%*c%f%*c%f",&x,&y,&z);
     if(n!=3) {printf("Fehler beim Einlesen von \"%s\": zeile=\"%s\" n=%d i=%d\n",quelle,zeile,n,i);}
     punkte[i].set(x,y,z);
     j=minmax(x,xmin,xmax);
     j=minmax(y,ymin,ymax);
     j=minmax(z,zmin,zmax);
     if(j&1) {izmin=i;}
     if(j&2) {izmax=i;}
    }
   if(i==1) //TODO: sicherstellen dass es mindestens 2 Punkte sind
    {
     x *= 1.000005; //y *= 1.000001;
     punkte[i].set(x,y,z);
     /*
     j=minmax(x,xmin,xmax);
     j=minmax(y,ymin,ymax);
     j=minmax(z,zmin,zmax);
     if(j&1) {izmin=i;}
     if(j&2) {izmax=i;}
     i++;
     */
     //printf("Punkt zugefuegt %f %f %f\n",x,y,z);//test
    }
  }
 tiefsterpunkt=punkte[izmin];
 hoechsterpunkt=punkte[izmax];
 fclose(fp);
 //#ifdef _DEBUG
 if(vflag) {
 printf("Bereich x: %f bis %f\n",xmin,xmax);//test
 printf("Bereich y: %g bis %g\n",ymin,ymax);//test
 printf("Bereich z: %f bis %f\n",zmin,zmax);//test
 printf("Index tiefster Punkt: %d\n",izmin);//test
 printf("Index hoechster Punkt: %d\n",izmax);//test
 }
 //#endif
 
 float xmitte=(xmin+xmax)/2, ymitte=(ymin+ymax)/2;
 float geoscal = 1.0;
 //float mingpxdist=0.002; //in main.cpp definiert
 if(xmax-xmin < mingpxdist)
  {
   if(xmax<xmin) printf("Fehler: xmax<xmin\n");//test
   xmin = xmitte - mingpxdist/2;
   xmax = xmitte + mingpxdist/2;
  }
 if(ymax-ymin < mingpxdist)
  {
   if(ymax<ymin) printf("Fehler: ymax<ymin\n");//test
   ymin = ymitte - mingpxdist/2;
   ymax = ymitte + mingpxdist/2;
  }
 if(bflag) swissflag=true;
 else
  {
   swissflag = (xmitte>=5.895 && xmitte<=10.745 && ymitte>=45.807 && ymitte<=47.822);
   if(swissflag)
    {//alles mit wgs2swiss() umrechnen, und dann bflag setzen:
     wgs2swiss(xmin,ymin); wgs2swiss(xmax,ymax);
     xmitte=(xmin+xmax)/2; ymitte=(ymin+ymax)/2;
     wgs2swiss(tiefsterpunkt.x, tiefsterpunkt.y);
     wgs2swiss(hoechsterpunkt.x, hoechsterpunkt.y);
     for(i=0;i<npunkte;i++) wgs2swiss(punkte[i].x,punkte[i].y);
     for(uint32 i=0; i<linienzuege.size(); i++)
      for(uint32 j=0; j<linienzuege[i].numPunkte; j++)
       {wgs2swiss(linienzuege[i].punkte[j].x, linienzuege[i].punkte[j].y);}
     for(i=0;i<iwpkt;i++)
       {wgs2swiss(gpx_wegpunkte[i].x, gpx_wegpunkte[i].y);}
     bflag=1;
    }
   else if(geoflag) {geoscal = cos(ymitte*GRAD);}
  }
 if(vflag) printf("geoscal=%f\n",geoscal);//test
 if(hscal==0) {zoffset = 0.01*(zmax-zmin); hscal = (zmax==zmin) ? 1.0 : hmax/(zmax-zmin);}
 koord.init(xmitte,ymitte,zmin-zoffset,hscal);
 float dx=0,dy=0;
 if(geoflag)
      {
       dx=(xmax-xmin)/20; dy=(ymax-ymin)/20; //Rand fuer GPX-Track
      }
 else {dx=(xmax-xmin)/20; dy=(ymax-ymin)/20;}//Rand fuer 2D-IR
 if(xscalfaktor==0) xscalfaktor = (MAXXY-MINXY)/((xmax-xmin)+2*dx); //TODO
 if(yscalfaktor==0) yscalfaktor = (MAXXY-MINXY)/((ymax-ymin)+2*dy); //TODO
 if(geoflag)
  {
   if((xmax-xmin)*geoscal >= ymax-ymin)
        {yscalfaktor = xscalfaktor/geoscal;}
   else {xscalfaktor = yscalfaktor*geoscal;}
  }
 yscalfaktor = -yscalfaktor; //fuer Umrechnung von User-y nach Opengl-z Koordinate
 //printf("Fehlersuche: xscalfaktor=%f yscalfaktor=%f\n",xscalfaktor, yscalfaktor);//test
 //printf("     yscalfaktor*geoscal=%f\n",yscalfaktor*geoscal);//test
 tiefsterpunkt.scal();
 hoechsterpunkt.scal();
#ifdef _DEBUG
 printf("Skalierfaktoren: %f %f hscal=%f\n",xscalfaktor,yscalfaktor,hscal);//test
 printf("Skalierte Werte:\n");//test
 printf("Tiefster Punkt:  x=%f y=%f z=%f\n",tiefsterpunkt.x, tiefsterpunkt.y, tiefsterpunkt.z);//test
 printf("Hoechster Punkt: x=%f y=%f z=%f\n",hoechsterpunkt.x, hoechsterpunkt.y, hoechsterpunkt.z);//test
#endif

#ifdef _DEBUG
 float hy1=1e20,hy2= -1e20;//test
#endif
 for(i=0;i<npunkte;i++)
  {
   punkte[i].scal(); //Punkte skalieren und somit in OpenGL-Koordinatensystem umwandeln
#ifdef _DEBUG
   if(punkte[i].y<hy1) hy1=punkte[i].y;//test
   if(punkte[i].y>hy2) hy2=punkte[i].y;//test
#endif
  }
 //Linienzuege skalieren nach weiter hinten verschoben
#ifdef _DEBUG
 printf("alle %d Punkte skaliert\n",npunkte);//test
 printf("kleinster groesster skaliert: %f %f\n",hy1,hy2);//test
#endif
 if(geoflag)//test
  {
   const double r=6.38e6; //Erdradius in Meter
   //double sx=cos(47.35*GRAD)*1000.0*sin(0.001*GRAD)*r;
   double sx = (bflag) ? 1.0 : geoscal*1000.0*sin(0.001*GRAD)*r;
   double xscaltest=sx/xscalfaktor;
#ifdef _DEBUG
   double sz=1000.0*sin(0.001*GRAD)*r;
   double zscaltest=sz/yscalfaktor;
   double hscaltest= (hscal==0) ? 1.0 : 1.0/hscal;
   printf("Skalierfaktoren fuer GPX-Umrechnung in Meter: %f %f %f\n",xscaltest,zscaltest,hscaltest);//test
   double dist = sqrt(sq(punkte[0].x-punkte[npunkte-1].x)+sq(punkte[0].z-punkte[npunkte-1].z));
   double hdist = punkte[npunkte-1].y - punkte[0].y;
   float h1=koord.yopengl2userz(punkte[0].y);
   float h2=koord.yopengl2userz(punkte[npunkte-1].y);
   printf("Hoehe Start: %f  Ziel: %f\n",h1,h2);//test
   printf("Distanz Luftlinie Start->Ziel: %f Meter\n",dist*xscaltest);//test
   printf("Hoehendifferenz Start->Ziel: %f Meter\n",hdist*hscaltest);//test
#endif
   printf("Vorschlag fuer gpxhscal: %f\n",1.0/(xscaltest*hscal));//TODO: falsche Formel!?
  }

 //#define MIT_ZWISCHENPUNKTEN
 if(geoflag) //TODO: Zwischenpunkte berechnen bei GPX-Track
  {
   float dsmax=0;
   for(int i=1;i<npunkte;i++)//test
    {
     float ds=sqrt(sq(punkte[i].x-punkte[i-1].x)+sq(punkte[i].z-punkte[i-1].z));
     if(ds>dsmax) dsmax=ds;
    }
   printf("maximale Distanz zwischen 2 Punkten: %f\n",dsmax);//test
#ifdef MIT_ZWISCHENPUNKTEN
   int k=int(dsmax/0.5+0.5);
   if(k>1)
    {
     printf("es werden jeweils %d Zwischenwert eingefuegt. altes npunkte=%d\n",k-1,npunkte);//test
     Vector *puneu=new Vector[npunkte*k];
     int j=0;
     for(int i=0;i<npunkte-1;i++)
      {
       glm::vec3 delta = (punkte[i+1]-punkte[i]); delta *= 1.0f/k;
       puneu[j++] = punkte[i];
       for(int ik=1;ik<k;ik++)
	{
	 puneu[j] = punkte[i];
	 puneu[j++] += delta*((float)(ik));
	}
      }
     puneu[j++] = punkte[npunkte-1];
     delete[] punkte;
     punkte=puneu;
     npunkte=j;
     printf("Zwischenpunkte berechnet: npunkte=%d\n",npunkte);//test
    }
#endif
  }
 
 if(argflag['D']) //Swissalti3d
  {
   if(!swissflag) printf("Hoehendaten (Swissalti3d) nur fuer Schweizerkarte vorhanden\n");
   else
    {
     //printf("Testpunkt Swissalti3d: hscal=%f znullpunkt=%f\n",koord.hscal,koord.znullpunkt);//test
     int ix,iy; float xo,yo; double x,y;
     float minh=1e4, maxh=0;
     float dxyo=1.0*100/N100;
     for(iy=0,yo=MINXY; iy<=MAXB; iy++,yo+=dxyo)
      {
       y=koord.zopengl2usery(yo);
       for(ix=0,xo=MINXY; ix<=MAXB; ix++,xo+=dxyo)
	{
	 x=koord.xopengl2userx(xo);
	 float hoehe_ueber_meer = swissalti3d(x,y);
	 if(hoehe_ueber_meer < minh) minh=hoehe_ueber_meer;
	 if(hoehe_ueber_meer > maxh) maxh=hoehe_ueber_meer;
	 //rasterhoehen[ix][iy] = koord.zuser2opengly(hoehe_ueber_meer);
	 rasterhoehen.set(ix,iy, hoehe_ueber_meer);
	}
      }
     //Hoehenskalierung anpassen:
     koord.znullpunkt = minh;
     koord.hscal = hscal = hmax/(maxh-minh);
     for(iy=0; iy<=MAXB; iy++)
      for(ix=0; ix<=MAXB; ix++)
	{
	 //rasterhoehen[ix][iy] = koord.zuser2opengly(rasterhoehen[ix][iy]);
	 rasterhoehen.set(ix,iy, koord.zuser2opengly(rasterhoehen.get(ix,iy)));
	}
     rasterhoehen_gesetzt=true;
    }
  }
 
 //Startposition setzen: mit Umrechnung von xyz-Koordinatensystem in OpenGL2-Koordinatensystem:
 if(startpos_auto==2)
  {startposx=punkte[0].x; startposz=punkte[0].z; startposy=punkte[0].y+10;}
 else if(startpos_auto==3)
  {startposx=punkte[npunkte-1].x; startposz=punkte[npunkte-1].z; startposy=punkte[npunkte-1].y+10;}
 else if(startpos_auto==4)
  {startposx=hoechsterpunkt.x; startposz=hoechsterpunkt.z; startposy=hoechsterpunkt.y+10;}
 //if(argflag['D']) startposy = landschaft1.yf(startposx,startposz)+0.1;//test
 //erst in main.cpp machen, landschaft1 ist hier noch nicht initialisiert
 //printf("startpos: %f %f %f\n",startposx,startposy,startposz);//test
 int n1=npunkte-1;
 gpx_startpos = punkte[0];
 gpx_zielpos = punkte[n1];
 
 if(geoflag)
  {
   for(uint32 i=0; i<linienzuege.size(); i++)
   {if(linienzuege[i].typ==0) //wenn noch in GPX-Koordinaten (sollte bei allen der Fall sein)
    {
     //printf("linienzuege[%d] wird von GPX in OpenGL-Koordinaten umgewandelt\n",i);//test
     for(uint32 j=0; j<linienzuege[i].numPunkte; j++)
      {
       if(argflag['D'])
	{//Hoehe auf Landschaftshoehe anpassen
	 double x=linienzuege[i].punkte[j].x;
	 double y=linienzuege[i].punkte[j].y;
	 linienzuege[i].punkte[j].z = swissalti3d(x,y);
	}
       vec3scal(linienzuege[i].punkte[j]);
       if(argflag['D']) linienzuege[i].punkte[j].y += 0.1;
       else linienzuege[i].punkte[j].y += 0.5;//alles etwas hoeher damit nicht in Boden verschwindend
      }
     linienzuege[i].typ=1; //jetzt in OpenGL-Koordinaten
    }
    else printf("Fehler: linienzuege[%d] schon skaliert!?\n",i);//test
   }
   printf("alle Linienzuege skaliert.\n");//test
  }
 
 if(!rasterhoehen_gesetzt)
  {
   printf("Testpunkt rasterhoehen von GPX-Track setzen\n");//test
   for(int iy=0; iy<=MAXB; iy++)
    for(int ix=0; ix<=MAXB; ix++)
     {
      //rasterhoehen[ix][iy] = naechster_punkt_suchen(ix,iy,punkte,npunkte,minaquadrat);
      rasterhoehen.set(ix,iy, naechster_punkt_suchen(ix,iy,punkte,npunkte,minaquadrat));
     }
   //TODO: Kanten glaetten
  }
 float Gxp=0,Gyp=0;
 float Gxpj[11],Gypj[11];//TODO
 bool kartenzusammensetzen = (qflag!=0 && qflagproz==101);//TODO
 if(qflag)
  {
   int i=(npunkte-1)*qflagproz/100; if(i>=npunkte) i=npunkte-1; else if(i<0) i=0;
   //printf("qflag=%d qflagproz=%d i=%d\n",qflag,qflagproz,i);//test
   Gxp = koord.xopengl2userx(punkte[i].x);
   Gyp = koord.zopengl2usery(punkte[i].z);
   //printf(" Gxp=%f Gyp=%f\n",Gxp,Gyp);//test
   if(kartenzusammensetzen)
    {for(int q=0,j=0; q<=100; q+=10,j++)
      {
       i=(npunkte-1)*q/100;
       Gxpj[j] = koord.xopengl2userx(punkte[i].x);
       Gypj[j] = koord.zopengl2usery(punkte[i].z);
      }
    }
  }
 delete[] punkte;
 
 FILE *fp2=myfopen(ziel,"w");
 if(fp2==NULL) {printf("Fehler2\n"); return;}//test
 fprintf(fp2,"xminxmax: %d %d\n",MINXY,MAXXY);
 fprintf(fp2,"yminymax: %d %d\n",MINXY,MAXXY);
 if(N100!=100)
  fprintf(fp2,"dxdy: %f %f\n",100.0/N100,100.0/N100);
 fprintf(fp2,"maxhoehe: %f\n",zmax);
#ifdef _DEBUG
 int np=0;//test
#endif
 for(int iy=0; iy<=MAXB; iy++)
  for(int ix=0; ix<=MAXB; ix++)
   {
    fprintf(fp2,"%f\n",rasterhoehen.get(ix,iy));
#ifdef _DEBUG
    np++;//test
#endif
   }
 fclose(fp2);
 gpx_gipfelpos = hoechsterpunkt;
 gpx_tiefpunkt = tiefsterpunkt; //wird in objekte.h noch an Landschaft1 angepasst
 //#ifdef _DEBUG
 if(vflag) {
 printf("gpx_gipfelpos gesetzt: x=%f y=%f z=%f\n",gpx_gipfelpos.x,gpx_gipfelpos.y,gpx_gipfelpos.z);//test
 printf("gpx_tiefpunkt gesetzt: x=%f y=%f z=%f\n",gpx_tiefpunkt.x,gpx_tiefpunkt.y,gpx_tiefpunkt.z);//test
 printf("Testpunkt: gpx_wegpunkte[%d] werden skaliert\n",iwpkt);//test
 }
 //#endif
 for(int i=0;i<iwpkt;i++)
  {
   vec3scal(gpx_wegpunkte[i]);
  }
#ifdef _DEBUG
 printf("%d Punkte geschrieben\n",np);//test
 printf("rasterdatei erstellt.\n");//test
#endif
 if(geoflag)
  {//Werte der Eckpunkte fuer Kartenueberlagerung:
   Gxmin = koord.xopengl2userx(-100.0f); Gxmax = koord.xopengl2userx(100.0f);
   Gymin = koord.zopengl2usery(100.0f);  Gymax = koord.zopengl2usery(-100.0f);
   if(bflag)
    {//fuer die Kartengrenze auf GPX-Koordinaten umrechnen
     swiss2wgs(Gxmin,Gymin);
     swiss2wgs(Gxmax,Gymax);
     if(qflag) swiss2wgs(Gxp,Gyp);
     if(kartenzusammensetzen)
      {
       for(int j=0; j<=10; j++)
	swiss2wgs(Gxpj[j],Gypj[j]);
      }
    }
   denkmal_scalefaktor_setzen(Gxmax-Gxmin);//TODO
   if(kflag==1)
    {//wenn Kartenname direkt angegeben wurde
     karte_laden("init",1);
     karte_laden("init",2); //beide Kartentabellen laden
    }
   else if(kflag>=11)
    {//wenn Kartennummer direkt angegeben wurde
     sprintf(kartenname,"karte%d.jpg",kflag);
     int nr=karte_download(0,kflag,0,0,kartenname,400);
     if(nr==0) {printf("Fehler: %s nicht gefunden\n",kartenname);}
    }
   else if(kflag==2)
    {
     bool gefunden=false;
     for(int knr=1; gefunden==false && knr<=2; knr++)
      {//zuerst in kartentabelle suchen:
       karte_laden("init",knr);
       int j;
       for(j=0;j<numKarten;j++)
	{
	 if(Gxmin >= kartentabelle[j].lonmin && Gxmin <= kartentabelle[j].lonmax &&
	    Gxmax >= kartentabelle[j].lonmin && Gxmax <= kartentabelle[j].lonmax &&
	    Gymin >= kartentabelle[j].latmin && Gymin <= kartentabelle[j].latmax &&
	    Gymax >= kartentabelle[j].latmin && Gymax <= kartentabelle[j].latmax)
	  {strcpy(kartenname, kartentabelle[j].name); break;}
	}
       if(j==numKarten) //wenn nicht in kartentabelle gefunden:
	{
	 //automatisch richtige Swisstopo-Karte auswaehlen und ev. downloaden:
	 //int nr=karte_download(Gxmin,Gxmax,Gymin,Gymax,kartenname,400);
	 int nr;
	 if(qflag)
	  {
	   if(kartenzusammensetzen) //TODO
	    {
	     for(j=0;j<=10;j++)
	      {
	       nr=karte_download(Gxpj[j],Gxpj[j],Gypj[j],Gypj[j],kartenname,400);
	       nr=karten_zusammensetzen(j,nr,kartenname);//TODO
	       //printf("Nummer von zusammgesetzter Karte: %d\n",nr);//test
	      }
	    }
	   else
	    nr=karte_download(Gxp,Gxp,Gyp,Gyp,kartenname,400);
	  }
	 else nr=karte_download(Gxmin,Gxmax,Gymin,Gymax,kartenname,400);
	 if(nr==0) {gefunden=false;}
	 //else {gefunden=true;}
	 else {gefunden=true; printf("passende Karte: %s\n",kartenname);}//test
	}
       //else {gefunden=true;}
       else {gefunden=true; printf("passende Karte aus Tabelle: %s\n",kartenname);}//test
      }
     if(!gefunden) printf("keine passende Karte gefunden\n");//test
    }
  }
}

//float sq(float x) {return x*x;} //schon in landschaften.h definiert, als double

float naechster_punkt_suchen(int ix,int iy,Vector punkte[],int npunkte,float minaquadrat)
{ //naechst gelegener Punkt zu ix,iy suchen und entsprechender z-Wert zurueckgeben
 float aquadrat, best[4]={FMAX,FMAX,FMAX,FMAX}, zbest[4]={0,0,0,0};
 for(int i=0;i<npunkte;i++)
  {
   //aquadrat = sq(punkte[i].x-(ix+MINXY))+sq(punkte[i].y-(iy+MINXY)); //alte Version
   aquadrat = sq(punkte[i].x-(ix+MINXY))+sq(punkte[i].z-(iy+MINXY));
   //da punkte[i] jetzt schon skaliert ist, ist y jetzt die Hoehe und z was vorher y-Koordinate war
   
   //float abstand = sqrt(aquadrat);
   if(aquadrat<best[0])
    {
     for(int j=4;--j>=1;) {best[j]=best[j-1]; zbest[j]=zbest[j-1];} //nachruecken
     //best[0]=aquadrat; zbest[0]=punkte[i].z; //neuer Bestwert //alte Version
     best[0]=aquadrat; zbest[0]=punkte[i].y; //neuer Bestwert (Hoehenwert)
    }
  }
 if(best[0]<=minaquadrat) return zbest[0]; //test
 else return 0.0;//test
 double gewichtungsumme=0,eins=1.0; //TODO: besserer Algorithmus suchen
 for(int j=0;j<4;j++)
  {
   if(best[j]<0.001) return zbest[j];
   if(best[j]<minaquadrat) gewichtungsumme += eins/best[j];
  }
 double z=0;
 for(int j=0;j<4;j++)
  {
   if(best[j]<minaquadrat)  z += zbest[j]/(gewichtungsumme*best[j]);
  }
 return (float)z;
}

void get_topo3dlink(char *ziel,const char *csvdatei,int x0,int y0)
{
 FILE *fp=myfopen(csvdatei,"r"); if(fp==NULL) {printf("csvdatei nicht gefunden\n"); *ziel=0; return;}
 char zeile[400], x0y0string[40];
 sprintf(x0y0string,"%04d-%04d",x0,y0);
 while(getline(fp,zeile,400))
  {
   if(strstr(zeile,x0y0string)!=NULL)
    {mystrncpy(ziel,zeile,400); fclose(fp); return;}
  }
 fclose(fp);
 printf("get_topo3dlink() Fehler: Link mit \"%s\" nicht gefunden\n",x0y0string);//test
 *ziel=0;//nicht gefunden
}

static float tmphdat[500][500];

float zeile_swissalti3d_hoehe_lesen(FILE *fp)
{
 char zeile[80];
 float x=0,y=0,hoehe=0;
 if(!getline(fp,zeile,80)) return 0.0;
 if(!isdigit(*zeile)) getline(fp,zeile,80);
 sscanf(zeile,"%f %f %f",&x,&y,&hoehe);
 if(hoehe== -9999.0) hoehe=193-5; //wenn Hoehendaten fehlen 5m unter tiefster Punkt der Schweiz setzen
 if(hoehe<100.0 || hoehe>6000)
  {printf("Fehler: unrealistische Hoehe: %f\n",hoehe);//test
   printf(" in Zeile: \"%s\"\n",zeile);//test
   hoehe=193-5;
  }
 return hoehe;
}

float swissalti3d(double x,double y)
{
 static int geladen_x0= -1, geladen_y0= -1;
 unsigned char c1,c2;
 if(x<2e6) x+=2.6e6;
 if(y<1e6) y+=1.2e6;
 //static int testflag=10;//test
 //bool testflag2 = (x>2.6e6+87600 && x<2.6e6+87612 && y>1.2e6+47148 && y<1.2e6+47169); //test
 int x0=x/1000, y0=y/1000;
 if(geladen_x0!=x0 || geladen_y0!=y0)
  {
   char toponame[400], tmp[400];
   sprintf(toponame,"karten/topo/topo%04d-%04d.hdat",x0,y0);
   if(vorhanden(toponame))
    {
     FILE *fp=myfopen(toponame,"rb");
     for(int i=0;i<500;i++)
      for(int j=0;j<500;j++)
       {
	c1=getc(fp); c2=getc(fp);
	uint16_t hoehe_in_dezimeter = (c1<<8)+c2;
	tmphdat[i][j] = hoehe_in_dezimeter*0.1;
       }
     fclose(fp);
     geladen_x0=x0; geladen_y0=y0;
    }
   else
    {
     get_topo3dlink(tmp,"ch.swisstopo.swissalti3d.csv",x0,y0);
     if(*tmp==0)
      {
       float hoehe = 193-5; //5m tiefer als tiefster Punkt der Schweiz: 193m Langensee
       for(int j=500;--j>=0;)
	for(int i=0;i<500;i++)
	  tmphdat[i][j] = hoehe;
       printf("keine Hoehendaten fuer diesen Bereich vorhanden - auf %.1f gesetzt\n",hoehe);//test
      }
     else
      {
       https_get(tmp,"tmp.zip");
       if(!vorhanden("tmp.zip")) {printf("konnte tmp.zip nicht erhalten\n"); exit(1);}//test
#ifdef _WIN32
       sprintf(tmp, "7za.exe x tmp.zip");
#else
       sprintf(tmp, "unzip tmp.zip");
#endif
       int err=system(tmp);
       if(err!=0) {printf("Fehler mit \"unzip tmp.zip\"\n"); exit(1);}//test
       unlink("tmp.zip");
       sprintf(tmp, "SWISSALTI3D_2_XYZ_CHLV95_LN02_%04d_%04d.xyz", x0,y0);
       FILE *fp=myfopen(tmp,"r");
       if(fp==NULL) {printf("Feher: Datei von unzip falsch?\n"); exit(1);}//test
       for(int j=500;--j>=0;)
	for(int i=0;i<500;i++)
	 {
	  float hoehe=zeile_swissalti3d_hoehe_lesen(fp);
	  if(hoehe==0) {printf("Fehler in Datei \"%s\"\n",tmp); exit(1);}//test
	  tmphdat[i][j] = hoehe;
	 }
       fclose(fp);
       unlink(tmp); //"SWISSALTI...xyz" loeschen
      }
     geladen_x0=x0; geladen_y0=y0;
     FILE *fp=myfopen(toponame,"wb"); if(fp==NULL) return 0;
     for(int i=0;i<500;i++)
      for(int j=0;j<500;j++)
       {
	uint16_t hoehe_in_dezimeter = int(tmphdat[i][j]*10+0.5);
	c1 = (hoehe_in_dezimeter >> 8); c2 = (hoehe_in_dezimeter & 0xFF);
	putc(c1,fp); putc(c2,fp);
       }
     fclose(fp);
    }
  }
 x -= x0*1000.0;
 y -= y0*1000.0;
 // x und y sind jetzt im Bereich zwischen 0 und 1000
 // die Daten in tmphdat[i][j] entsprechen den Hoehen an den Punkten 1, 3, 5 ... 999 in beiden Richtungen
 // Um die richtige Hoehe zu finden wird die Umkehrung der folgenden Formel benoetigt:
 //  tmphdat[i][j] = Hoehe(x=1+2*i, y=1+2*j)
 double x10=(x-1.0)/2;
 double y10=(y-1.0)/2;
 int i0,i1,j0,j1;
 if(x10<0) {i0=i1=0;} else {i0=int(x10); i1=i0+1; if(i1==500) i1--;}
 if(y10<0) {j0=j1=0;} else {j0=int(y10); j1=j0+1; if(j1==500) j1--;}
 x10 = 1+2*i0;
 y10 = 1+2*j0;
 float h0,h1;
 /* */
 h0=tmphdat[i0][j0]; h1=tmphdat[i1][j0];
 float hj0 = h0+(x-x10)*(h1-h0)/2.0; //Interpolation in x-Richtung an Stelle j0
 h0=tmphdat[i0][j1]; h1=tmphdat[i1][j1];
 float hj1 = h0+(x-x10)*(h1-h0)/2.0; //Interpolation in x-Richtung an Stelle j1
 float hoehe = hj0+(y-y10)*(hj1-hj0)/2.0; //Interpolation in y-Richtung
 /* if(testflag && testflag2)
  {
   printf("Testpunkt Interpolation: %f %f %f %f\n",tmphdat[i0][j0],tmphdat[i1][j0],tmphdat[i0][j1],tmphdat[i1][j1]);
   printf("  i0=%d j0=%d i1=%d j1=%d\n",i0,j0,i1,j1);//test
   printf("  aktueller Punkt: x=%f y=%f\n",x,y);
   printf("  abgerundete Pkt: x10=%f y10=%f\n",x10,y10);
   printf("  Interpoliert x-Richt.: hj0=%f hj1=%f\n",hj0,hj1);
   printf("  Interpoliert y-Richt.: hoehe=%f\n",hoehe);
   testflag--;
   } */
 /* */
 /* test: zuerst in y-Richtung interpolieren * /
 h0=tmphdat[i0][j0]; h1=tmphdat[i0][j1];
 float hj0 = h0+(y-y10)*(h1-h0)/2.0; //Interpolation in y-Richtung an Stelle i0
 h0=tmphdat[i1][j0]; h1=tmphdat[i1][j1];
 float hj1 = h0+(y-y10)*(h1-h0)/2.0; //Interpolation in y-Richtung an Stelle i1
 float hoehe = hj0+(x-x10)*(hj1-hj0)/2.0; //Interpolation in x-Richtung
 / * test: zuerst in y-Richtung interpolieren */
 
 return hoehe;
}

bool vorhanden_in_liste(int n,int nrliste[],int imax)
{
 for(int i=0;i<imax;i++)
  if(n==nrliste[i]) return true;
 return false;
}

int karten_zusammensetzen(int j,int nr,char *kartenname) //j muss von 0 bis 10 laufen, Rueckgabewwert=nr
{
 static int nrliste[11], imax=0;
 static char *nameliste[11];
 int i;
 //printf("karten_zusammensetzen(j=%d, nr=%d, kartenname=\"%s\")\n",j,nr,kartenname);//test
 if(j==0) imax=0;
 for(i=0; i<imax; i++) {if(nrliste[i]==nr) break;}
 if(i==imax) //wenn noch nicht in Liste
  {
   nrliste[imax] = nr; //in Liste aufnehmen
   nameliste[imax] = new char[strlen(kartenname)+1];
   strcpy(nameliste[imax],kartenname);
   imax++;
  }
 if(j==10) //wenn alle Kartenteile fertig gesucht
  {
   if(imax<2) return nr;//kein Zusammensetzen noetig
   printf("Karten zusammensetzen:\n");//test
   for(i=0;i<imax;i++) printf(" nr=%d name=\"%s\"\n",nrliste[i],nameliste[i]);//test
   if(qflag==25) //Karten 1:25000
    {
     char command[1024];
     int inr1=0, inr2= -1;
     if(imax>9) {printf("bisher auf 9 Karten beschraenkt\n"); imax=9;}//test
     for(i=1; i<imax; i++) //kleinste Nummer suchen
      {if(nrliste[i] < nrliste[inr1]) inr1=i;}
     for(i=0; i<imax; i++) //groesste Nummer suchen
      {if(i!=inr1 && (inr2<0 || nrliste[i] > nrliste[inr2])) inr2=i;}
     int nr1=nrliste[inr1], nr2=nrliste[inr2];//kleinste und groesste Nummer 
     if(imax==2) //nur 2 Karten zusammensetzen
      {
       printf("genau 2 Karten zusammensetzen\n");//test
       sprintf(kartenname,"karten/karte%d%d.jpg",nr1,nr2);
       nr = nr1*10000+nr2;
       if(nr2-nr1 == 1) //horizontale Anordnung
	{
	 sprintf(command,"%s %s %s -geometry +0+0 -resize 50%% %s", MONTAGE_COMMAND,
		 nameliste[inr1], nameliste[inr2], kartenname);
	}
       else if(nr2-nr1 == 20) //vertikale Anordnung
	{
	 sprintf(command,"%s %s %s -tile x2 -geometry +0+0 -resize 50%% %s", MONTAGE_COMMAND,
		 nameliste[inr1], nameliste[inr2], kartenname);//TOTO: eventuell -tile x4?
	}
       else {printf("Fehler: kann Karten %d und %d nicht zusammensetzen\n",nr1,nr2); return nr;}//test
      }
     else //mehr als 2 Karten zusammensetzen: 4, 6, oder 9
      {
       printf("mehr als 2 Karten zusammensetzen: imax=%d\n",imax);//test
       int inr3=0,inr4= -1;
       for(i=1; i<imax; i++) //kleinste Nummer in hinterer Stelle suchen
	{if(nrliste[i]%20 < nrliste[inr3]%20) inr3=i;}
       for(i=0; i<imax; i++) //groesste Nummer in hinterer Stelle suchen
	{if(i!=inr3 && (inr4<0 || nrliste[i]%20 > nrliste[inr4]%20)) inr4=i;}
       int nr3=nrliste[inr3], nr4=nrliste[inr4];
       int jxmax=nr4%20-nr3%20;
       int jymax=nr2/20-nr1/20;
       int nr5 = nr1/20*20 + nr3%20;
       for(int jx=0;jx<=jxmax;jx++)
	for(int jy=0;jy<=jymax;jy++)
	 {
	  int n=nr5+jx+jy*20;
	  if(!vorhanden_in_liste(n,nrliste,imax))
	   {
	    if(imax>=9) printf("Fehler: zu viele Teilkarten\n");
	    else
	     {
	      nrliste[imax]=n; nameliste[imax]=new char[400];
	      karte_download(0,n,0,0,nameliste[imax],400);
	      imax++;
	     }
	   }
	 }
       //neu sortieren:
       i=0;
       for(int jy=0;jy<=jymax;jy++)
	for(int jx=0;jx<=jxmax;jx++)
	 {
	  int n=nr5+jx+jy*20;
	  nrliste[i] = n;
	  sprintf(nameliste[i],"karten/karte%d.jpg",n);
	  printf("neu sortiert i=%d: nr=%d name=\"%s\"\n",i,n,nameliste[i]);//test
	  i++;
	 }
       sprintf(kartenname,"karten/karte%d%d.jpg",nrliste[0],nrliste[imax-1]);
       if(imax==4)
	sprintf(command,"%s %s %s %s %s -tile x2 -geometry +0+0 -resize 50%% %s", MONTAGE_COMMAND,
		nameliste[0], nameliste[1], nameliste[2], nameliste[3], kartenname);//TODO
       else if(imax==9)
	sprintf(command,"%s %s %s %s %s %s %s %s %s %s -tile x3 -geometry +0+0 -resize 33%% %s", MONTAGE_COMMAND,
		nameliste[0], nameliste[1], nameliste[2], nameliste[3], nameliste[4], nameliste[5],
		nameliste[6], nameliste[7], nameliste[8], kartenname);//TODO
       else if(imax==6)
	{
	 int xn = (jxmax>jymax) ? 2 : 3;
	 sprintf(command,"%s %s %s %s %s %s %s -tile x%d -geometry +0+0 -resize 33%% %s", MONTAGE_COMMAND,
		 nameliste[0], nameliste[1], nameliste[2],
		 nameliste[3], nameliste[4], nameliste[5], xn, kartenname);//TODO
	}
      }
     if(!vorhanden(kartenname))
      {
       printf("system(%s)\n",command);//test
       printf("Karten werden zusammengesetzt ...\n");//test
       system(command);
      }
     else printf("zusammengesetzte Karte \"%s\" schon vorhanden\n",kartenname);//test
    }
   else
    {
     printf("Karten zusammensetzen geht bisher nur fuer 1:25'000 Karten\n");//test
     //TODO
     //TODO: nr neu setzen, kartenname neu setzen
    }
   //aufraeumen:
   for(i=0;i<imax;i++)
    {delete[] nameliste[i];}
   imax=0;
  }
 return nr;
}
