Option Declare

Const RADIAN = Pi / 180  
Const MEAN_RADIUS = 6371.299 ' mittlerer Erdradius in km     
Const EQUATOR_RADIUS = 6378.14 ' Äquatoradius in km
Const FLATTENING_FACTOR = 1 / 298.257 ' Abplattung der Erde

Sub Click(Source As Button)
	
	'Mettmann    06E59      51N15
	'Düsseldorf   06E47      51N14	
	Msgbox  dblDistanceBetween2PointsOnEarth ( "51N15" , "06E59" , "51N15" , "06E59" )
	
End Sub

Function dblDistanceBetween2PointsOnEarth(  c1N As String, c1E As String, c2N As String, c2E As String ) As Double 
	'Die Funktion liefert die Entfernung zweier Punkte auf der Erdoberfläche in km. 
	'Das Verfahren berücksichtigt die Erdabplattung und ist auch für sehr kleine Entfernungen genau.
	'Keine Berücksichtigung der Höhe über N.N.! 
	'Quelle: Jean Meeus - Astronomische Algorithmen (Auflage 1992, Verlag Johann Ambrosius Barth), Seite 93f, 118ff 
	
	'Variables 
	Dim phi1_rad    		As Double
	Dim lambda1_rad 		As Double
	Dim phi2_rad    		As Double
	Dim lambda2_rad 		As Double
	Dim F           		As Double
	Dim sinF2       		As Double
	Dim cosF2       		As Double
	Dim G           		As Double
	Dim sinG2       		As Double
	Dim cosG2       		As Double
	Dim dl          		As Double
	Dim s           		As Double
	Dim C          			As Double
	Dim om          		As Double
	Dim R           		As Double
	Dim cosD        		As Double
	Dim D           		As Double
	Dim H1          		As Double
	Dim H2          		As Double 
	Dim DEG 			As String
	Dim MMSS 			As String
	Dim DIRECTION  			As String
	Dim RES				As Double
	
	'Convert to radians     
	DIRECTION = Mid$(c1N,3,1)
	MMSS= atWord (c1N,DIRECTION,2)
	DEG = atWord (c1N,DIRECTION,1)
	Res = Radians (DEG, MMSS)
	
	Select Case DIRECTION
	Case "N"
		phi1_rad = RES
	Case "S"
		phi1_rad = -RES	
	End Select
	
	DIRECTION = Mid$(c1E,3,1)
	MMSS= atWord (c1E,DIRECTION,2)
	DEG = atWord (c1E,DIRECTION,1)
	Res = Radians (DEG, MMSS)
	
	Select Case DIRECTION	
	Case "E"
		lambda1_rad = RES
	Case "W"
		lambda1_rad = -RES
	End Select
	
	DIRECTION = Mid$(c2N,3,1)
	MMSS= atWord (c2N,DIRECTION,2)
	DEG = atWord (c2N,DIRECTION,1)
	Res = Radians (DEG, MMSS)
	
	Select Case DIRECTION
	Case "N"
		phi2_rad = RES
	Case "S"
		phi2_rad = - RES
	End Select
	
	DIRECTION = Mid$(c2E,3,1)
	MMSS= atWord (c2E,DIRECTION,2)
	DEG = atWord (c2E,DIRECTION,1)
	Res = Radians (DEG, MMSS)
	
	Select Case DIRECTION	
	Case "E"
		lambda2_rad = RES
	Case "W"		
		lambda1_rad = - RES
	End Select	
	
	'Check on precise calculation     
	cosD = Sin(phi1_rad) * Sin(phi2_rad) + Cos(phi1_rad) * Cos(phi2_rad) * Cos(lambda1_rad - lambda2_rad) 
	If cosD < 0.999995 Then         
	'Distance is long 
		F = (phi1_rad + phi2_rad) / 2         
		sinF2 = (Sin(F)) ^ 2         
		cosF2 = (Cos(F)) ^ 2
		G = (phi1_rad - phi2_rad) / 2
		sinG2 = (Sin(G)) ^ 2        
		cosG2 = (Cos(G)) ^ 2         
		dl = (lambda1_rad - lambda2_rad) / 2 
		s = sinG2 * (Cos(dl)) ^ 2 + cosF2 * (Sin(dl)) ^ 2       
		C = cosG2 * (Cos(dl)) ^ 2 + sinF2 * (Sin(dl)) ^ 2         
		om = Atn(Sqr(s / C)) 
		R = Sqr(s * C) / om         
		H1 = (3 * R - 1) / (2 * C)         
		H2 = (3 * R + 1) / (2 * s) 
		dblDistanceBetween2PointsOnEarth = 2 * om * EQUATOR_RADIUS *  (1+ FLATTENING_FACTOR * H1 * sinF2 * cosG2  - FLATTENING_FACTOR * H2 * cosF2 * sinG2 ) 
	Else        
	 'Distance is short 
		dblDistanceBetween2PointsOnEarth = Sqr(((lambda2_rad - lambda1_rad) * Cos((phi1_rad + phi2_rad) / 2)) ^ 2 + (phi2_rad - phi1_rad) ^ 2) * MEAN_RADIUS 
	End If
End Function

Function atWord ( sourceString As String, separator As String, number As Integer ) As String
	Dim searchString As String, substring As String
	Dim i As Integer, pos As Integer
	searchString=SourceString & separator
	For i = 1 To number
		pos=Instr(searchString, separator)
		If pos=0 Then Exit For
		substring=Left(searchString,pos-1)
		searchString=Mid(searchString, pos+1)
	Next
	If pos > 0 Then
		atWord=substring
	Else
		atWord=""
	End If
End Function


Function Radians ( DEG As String, MMSS As String ) As Double
	If Len (MMSS) = 4 Then
		Radians = (DEG + (Left$(MMSS,2)  / 60) + (Right$(MMSS,2)  / 3600)) * RADIAN
	Else
		Radians = (DEG + (MMSS  / 60)) *RADIAN
	End If
End Function