#!/usr/bin/perl

#############################################################################
## >>> 2014 <<<<
## Martin Hauke (mardnh@gmx.de)
##
## Older GSM BTS's from O2 Germany (formerly MNC 07 - now MNC 03) send out their BTS
## coordinates via smscb (channel 221).
## With Alex Badea's patches for libosmocore/osmocombb applied you can log those smscb
## messages via osmocom-bb's cell_log application.

## This script parses osmocom-bb cell_log output and writes the relevant info to a
## nice kml-file.
##
## usage example:
## cat cell.log | o2-osmocom-bb-smscb2kml.pl > nice_map.kml
##
##
## The used Converter for Gauss - Krueger coordinates is licensed under the GPL and was
## written by
## - Andreas Achtzehn (achtzehn@linux-society.de)
## - Nobert Hüttisch (nobbi@nobbi.com)


#use strict;

my %seen;

my ($line, $time, $arfcn, $mcc, $mnc, $lac, $cid, $cb);
my ($viag_gk, $sr, $sx);
my $regex = qw(^CellBroadcast: TIME=(?<time>\d+)\sARFCN=(?<arfcn>\d+)\sMCC=(?<mcc>\d+)\sMNC=(?<mnc>\d+)\sLAC=(?<lac>.*)\sCID=(?<cid>.*)\sCB=\"(?<cb>.*)\"$);

&printKMLheader;

foreach $line (<STDIN>) {
	chomp($line);
	if ( $line =~ /$regex/ ) {
		#print "$line\n";
		$time  = "$+{time}";
		$arfcn = "$+{arfcn}";
		$mcc   = "$+{mcc}";
		$mnc   = "$+{mnc}";
		$lac   = "$+{lac}";
		$cid   = "$+{cid}";
		$cb    = "$+{cb}";
		$cid_dec = hex($cid);
		$lac_dec = hex($lac);
		$time_x  = scalar localtime($time);

		# check if O2-Germany (O2 MNC was until mid 2016 - since then 03)
		if ($mnc == "03" && $mcc == "262") {
			# did we seen that entry point already?
			next if $seen{"$arfcn,$mcc,$mnc,$lac,$cid"}++;
			$cb =~ /(?<gk>\d+).*/;
			$viag_gk = "$+{gk}";
			$xlat=substr($viag_gk,0,6)*10;
			$xlon=substr($viag_gk,6,6)*10;
			getCoordinate($xlat, $xlon);
			#print "Google Maps: https://maps.google.com/maps?z=12&h=m&q=loc:$wgs84lat+$wgs84lon&z=15\n";
			#print "$time $arfcn $mcc $mnc $lac $cid $viag_gk $wgs84lat $wgs84lon\n";
			
			# print placemarks
			print qq{\t\t<Placemark>\n};
			print qq{\t\t\t<name>CID:$cid_dec</name>\n};
			print qq{\t\t\t<description><![CDATA[CID: <b>$cid ($cid_dec)</b><br> LAC: <b>$lac ($lac_dec)</b><br>ARFCN: <b>$arfcn</b><br> MCC: <b>$mcc</b><br>MNC: <b>$mnc</b><br>TIME: <b>$time_x</b> ]]></description>\n};
			#print qq{\t\t\t</description>\n};
			#print "<styleUrl>#rxl69</styleUrl>"
			print qq{\t\t\t<Point>\n};
			print qq{\t\t\t\t<coordinates>$wgs84lon,$wgs84lat</coordinates>\n};
			print qq{\t\t\t</Point>\n};
			print qq{\t\t</Placemark>\n};
		};
	};
};

&printKMLfooter;

sub getCoordinate {
	$sr = shift;
	$sx = shift;
	&potsdam;
	&modolensky;
}

sub potsdam {
	# Errechnen der Position nach Potsdam-Koordinatensystem
	$bm  = int($sr/1000000);
	$y   = $sr-($bm*1000000+500000);
	$si  = $sx/111120.6196;
	$px  = $si+0.143885358*sin(2*$si*0.017453292)+0.00021079*sin(4*$si*0.017453292)+0.000000423*sin(6*$si*0.017453292);
	$t   = (sin($px*0.017453292))/(cos($px*0.017453292));
	$v   = sqrt(1+0.006719219*cos($px*0.017453292)*cos($px*0.017453292));
	$ys  = ($y*$v)/6398786.85;
	$lat = $px-$ys*$ys*57.29577*$t*$v*$v*(0.5-$ys*$ys*(4.97-3*$t*$t)/24);
	$dl  = $ys*57.29577/cos($px*0.017453292) * (1-$ys*$ys/6*($v*$v+2*$t*$t-$ys*$ys*(0.6+1.1*$t*$t)*(0.6+1.1*$t*$t)));
	$lon = $bm*3+$dl;

	$ph =int($lat);
	$pz =($lat-$ph)*60;
	$pm =int($pz);
	$ps =($pz-$pm)*60;

	$lh =int($lon);
	$lz =($lon-$lh)*60;
	$lm =int($lz);
	$ls =($lz-$lm)*60;
}

sub modolensky {
	$potsd_a  = 6377397.155;
	$wgs84_a  = 6378137.0;
	$potsd_f  = 1/299.152812838;
	$wgs84_f  = 1/298.257223563;

	$potsd_es = 2*$potsd_f - $potsd_f*$potsd_f;

	$potsd_dx = 606.0;
	$potsd_dy = 23.0;
	$potsd_dz = 413.0;
	$pi = 3.141592654;
	$latr = $lat/180*$pi;
	$lonr = $lon/180*$pi;

	$sa = sin($latr);
	$ca = cos($latr);
	$so = sin($lonr);
	$co = cos($lonr);

	$bda  = 1-$potsd_f;

	$delta_a = $wgs84_a - $potsd_a;
	$delta_f = $wgs84_f - $potsd_f;

	$rn = $potsd_a / sqrt(1-$potsd_es*sin($latr)*sin($latr));
	$rm = $potsd_a * ((1-$potsd_es)/sqrt(1-$potsd_es*sin($latr)*sin($latr)*1-$potsd_es*sin($latr)*sin($latr)*1-$potsd_es*sin($latr)*sin($latr)));

	$ta = (-$potsd_dx*$sa*$co - $potsd_dy*$sa*$so)+$potsd_dz*$ca;
	$tb = $delta_a*(($rn*$potsd_es*$sa*$ca)/$potsd_a);
	$tc = $delta_f*($rm/$bda+$rn*$bda)*$sa*$ca;
	$dlat = ($ta+$tb+$tc)/$rm;

	$dlon = (-$potsd_dx*$so + $potsd_dy*$co)/($rn*$ca);

	$wgs84lat = ($latr + $dlat)*180/$pi;
	$wgs84lon = ($lonr + $dlon)*180/$pi;
	$ph=int($wgs84lat);
	$pz=($wgs84lat-$ph)*60;
	$pm=int($pz);
	$ps=int(($pz-$pm)*60);
	
	$lh=int($wgs84lon);
	$lz=($wgs84lon-$lh)*60;
	$lm=int($lz);
	$ls=int(($lz-$lm)*60);

	return ($wgs84lat, $wgs84lon )
}

sub printKMLheader {
	print qq{<?xml version="1.0" encoding="UTF-8"?>\n};
	print qq{<kml xmlns="http://www.opengis.net/kml/2.2">\n};
	print qq{\t<Document>\n};
};

sub printKMLfooter {
	print qq{\t</Document>\n};
	print qq{</kml>\n};
};

# EOF
