function [X,Y,Z]=wgs2ch(phi,lambda,h)
% WGS2CH   Umrechnung von WGS84-Koordinaten in Schweizer Landeskoordinaten.
%
%   [X,Y,Z]=WGS2CH(PHI,LAMBDA,H) berechnet Schweizer Landeskoordinaten
%   X, Y und Z aus Breite PHI, Laenge LAMBDA und Hoehe H im WGS84-System.
%
%   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 CH2WGS.
   
% Copyright (c) 2001 by AMRON
% Norbert.Marwan AT gmx.net, Potsdam
% Last Revision: 2001-05-03


	phi=phi*pi/180;
	lambda=lambda*pi/180;

%%%%%%%%%%%

%% ELLIPS WGS -> GEOZENTR. WGS
	a=6378137.000;
	e=0.00669438000;

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

	xWGS=(Rn+h)*cos(phi)*cos(lambda);
	yWGS=(Rn+h)*cos(phi)*sin(lambda);
	zWGS=(Rn*(1-e)+h)*sin(phi);

%% GEOZENTR. WGS -> GEOZENTR. CH
	dX=-660.075;
	dY=-13.551;
	dZ=-369.340;
	M=0.99999436;
	alpha=(-2.485/10000)*pi/200;
	beta=(-1.783/10000)*pi/200;
	gamma=(-2.939/10000)*pi/200;

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

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

	zCH=dZ+M* (...
	   sin(beta)*xWGS + ...
	   -sin(alpha)*cos(beta)*yWGS + ...
	   cos(alpha)*cos(beta)*zWGS);

%% GEOZENTR. CH -> ELLIPS. CH
	a=6377397.155;
	e=0.006674372231;

	lCH=atan(yCH/xCH);
	
	phiCH=46.952405555*pi/180;
	er=1;
	while er>eps
	er=phiCH;
	Rn=a/sqrt(1-e*sin(phiCH)^2);
	hCH=sqrt(xCH^2+yCH^2)/cos(phiCH)-Rn;
	phiCH=atan((zCH/sqrt(xCH^2+yCH^2))/(1-Rn*e/(Rn+hCH)));
	er=abs(er-phiCH);
	end


%% ELLIPSE -> 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;
	
	S=(alpha*log(tan((pi/4)+(phiCH/2)))) - ((alpha*e/2)*log( (1+(e*sin(phiCH))) / (1-(e*sin(phiCH))) ))+K;

	phiCH=2*(atan(exp(S))-pi/4);
	lCH=alpha*(lCH-L0);

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

%% PROJEKTION AUF ZYLINDER
	Y=round((600000+R*lambda)*10)/10;
	X=round((200000+(R/2)*log((1+sin(phi))/(1-sin(phi))))*10)/10;
	Z=round(hCH*10)/10;

	
