function [phi,lambda,h]=ch2wgs(X,Y,Z)
% CH2WGS   Umrechnung Schweizer Landeskoordinaten in das WGS84-System.
%
%   [PHI,LAMBDA,H]=CH2WGS(X,Y,Z) berechnet die Breite PHI, 
%   Laenge LAMBDA und Hoehe H im WGS84-System aus den Schweizer
%   Landeskoordinaten.
%
%   Die Umrechnung basiert auf folgendem Material des Bundesamtes fuer
%   Landestopographie der Schweiz, Wabern:
%     1. Schweizerisches Projektionssystem. Formeln fuer die 
%        Umrechnung von Landeskoordinaten in geographische 
%	 Koordinaten und umgekehrt, 1984
%     2. Transformation von Landeskoordinaten CH-1903 in WGS-84 
%        Koordinaten, 1990
%     3. Ergaenzung zur Formelzusammenstellung fuer die Umrechnung
%        von WGS84-Koordinaten in Schweizerische Projektionskoordinaten
%
%   Siehe auch WGS2CH.
   
% Copyright (c) 2001 by AMRON
% Norbert.Marwan AT gmx.net, Potsdam
% Last Revision: 2001-05-03


%% PROJEKTION ZYLINDER -> KUGEL

	B0=0.81947406867611;
	L0=0.1298452241431;
	b0=0.81869435858167;
	e=sqrt(0.006674372230614);
	a=6377397.155;
	K=0.0030667323772751;
	alpha=1.00072913843038;
	R=6378815.90365;
	
	phi=2*(atan(exp((X-200000)/R))-pi/4);
	lambda=(Y-600000)/R;

%% PSEUDOMERIDIAN SYSTEM
	lCH=atan(sin(lambda)/(-sin(b0)*tan(phi)+cos(b0)*cos(lambda)));
	phiCH=asin(cos(b0)*sin(phi)+sin(b0)*cos(phi)*cos(lambda));


%% KUGEL -> ELLIPSE	

	phi=B0;
	er=1;
	while er>eps
	er=phi;
	S=(1/alpha)*(log(tan(pi/4+phiCH/2))-K)+((e/2)*log( (1+(e*sin(phi))) / (1-(e*sin(phi))) ));
	phi=2*(atan(exp(S))-pi/4);
	er=abs(er-phi);
	end

	l=lCH/alpha+L0;
	h=Z;

%% ELLIPSE CH -> GEOZENTR. CH
	e=0.006674372230614;
	a=6377397.155;

	Rn=a/sqrt(1-e*sin(phi)^2);

	xCH=(Rn+h)*cos(phi)*cos(l);
	yCH=(Rn+h)*cos(phi)*sin(l);
	zCH=(Rn*(1-e)+h)*sin(phi);
%% GEOZENTR CH -> GEOZENTR WGS
	dX=660.075;
	dY=13.551;
	dZ=369.344;
	M=1.00000566;
	alpha=(2.485/10000)*pi/200;
	beta=(1.783/10000)*pi/200;
	gamma=(2.939/10000)*pi/200;

	xWGS=dX+M* (...
	   cos(beta)*cos(gamma)*xCH + ...
	   (cos(alpha)*sin(gamma)+sin(alpha)*sin(beta)*cos(gamma))*yCH + ...
	   (sin(alpha)*sin(gamma)-cos(alpha)*sin(beta)*cos(gamma))*zCH);

	yWGS=dY+M* (...
	   -cos(beta)*sin(gamma)*xCH + ...
	   (cos(alpha)*cos(gamma)-sin(alpha)*sin(beta)*sin(gamma))*yCH + ...
	   (sin(alpha)*cos(gamma)-cos(alpha)*sin(beta)*sin(gamma))*zCH);

	zWGS=dZ+M* (...
	   sin(beta)*xCH + ...
	   -sin(alpha)*cos(beta)*yCH + ...
	   cos(alpha)*cos(beta)*zCH);
%% GEOZENTR. WGS -> ELLIPS. WGS
	a=6378137;
	e=0.00669438;

	lambda=atan(yWGS/xWGS);
	
	phi=46.952405555*pi/180;
	er=1;
	while er>eps
	er=phi;
	Rn=a/sqrt(1-e*sin(phi)^2);
	h=sqrt(xWGS^2+yWGS^2)/cos(phi)-Rn;
	phi=atan((zWGS/sqrt(xWGS^2+yWGS^2))/(1-Rn*e/(Rn+h)));
	er=abs(er-phi);
	end
	
	phi=phi*180/pi;
	lambda=lambda*180/pi;
	h=round(h*10)/10;

