#!/usr/bin/ruby

require 'shapelib'

sfp_main = ShapeLib::ShapeFile::new("airspace.shp", :Polygon, 
                               [['name', :String, 64],
                                ['type', :String, 20],
                                ['class', :String, 5]])

sfp_atz = ShapeLib::ShapeFile::new("atz.shp", :Polygon, 
                               [['name', :String, 64],
                                ['type', :String, 20],
                                ['class', :String, 5]])

sfp_points = ShapeLib::ShapeFile::new("points.shp", :Point, 
                               [['name', :String, 64],
                                ['type', :String, 20]])


class String
    def StartsWith(s)
        return false if !index(s)
        return self[(s.size)..-1]
    end
end

def parseLat(s)
    result = s[1..2].to_f + s[3..4].to_f/60.0 + s[5..6].to_f/3600.0
    result = -result if(s[0] == ?S)
    return result
end

def parseLon(s)
    result = s[1..3].to_f + s[4..5].to_f/60.0 + s[6..7].to_f/3600.0
    result = -result if(s[0] == ?W)
    return result
end

def pointDistance(lat1, lon1, lat2, lon2)
    lat1 = lat1 * 2 * Math::PI / 360
    lon1 = lon1 * 2 * Math::PI / 360
    lat2 = lat2 * 2 * Math::PI / 360
    lon2 = lon2 * 2 * Math::PI / 360

    distance=2*Math::asin(Math::sqrt(Math::sin((lat1-lat2)/2) ** 2 +
                 Math::cos(lat1)*Math::cos(lat2)*(Math::sin((lon1-lon2)/2) ** 2)));
    return distance*3438.461;
end

def bearing(lat1, lon1, lat2, lon2)
    lat1 = lat1 * 2 * Math::PI / 360
    lon1 = lon1 * 2 * Math::PI / 360
    lat2 = lat2 * 2 * Math::PI / 360
    lon2 = lon2 * 2 * Math::PI / 360

    angle = Math::atan2(Math::sin(lon2-lon1)*Math::cos(lat2),Math::cos(lat1)*Math::sin(lat2)-Math::sin(lat1)*Math::cos(lat2)*Math::cos(lon2-lon1));
    angle += 2*Math::PI while(angle < 0) 
    angle -= 2*Math::PI while(angle > 2*Math::PI) 
    return angle;
end

def newPoint(centrelat, centrelon, angle, r)
    centrelat = centrelat * 2 * Math::PI / 360
    centrelon = centrelon * 2 * Math::PI / 360
    newlat = Math::asin(Math::sin(centrelat) * Math::cos(r) + \
                  Math::cos(centrelat) * Math::sin(r) * Math::cos(angle));
    newlon=-centrelon-Math::asin(Math::sin(angle)*Math::sin(r)/Math::cos(newlat))+Math::PI;

    newlon -= Math::PI;
    newlon -= 2*Math::PI while(newlon > Math::PI) 
    newlon += 2*Math::PI while(newlon < -Math::PI) 
    newlon = -newlon;
    newlat = newlat * 360 / (2*Math::PI);
    newlon = newlon * 360 /(2*Math::PI);
    return [newlon, newlat] #Note order - this is used in shapefile
end

$SL_DEGREES = (5.0 * (2*Math::PI)/360.0)

def curveToSegments(fromlat, fromlon, centrelat, centrelon,tolat,tolon,reverse)
    points = Array.new
    startangle = bearing(centrelat, centrelon, fromlat, fromlon)
    endangle = bearing(centrelat, centrelon, tolat, tolon)
    endangle += 2*Math::PI if(startangle >= endangle && !reverse)
    startangle += 2*Math::PI if(startangle <= endangle && reverse)

    r = pointDistance(centrelat, centrelon, tolat, tolon)/3438.461
    anglejump = (endangle-startangle);
    nsegs = (anglejump/$SL_DEGREES).to_i.abs
    anglejump = anglejump>0 ? +$SL_DEGREES : -$SL_DEGREES
    nsegs-=1; startangle+=anglejump;
    loop {
        break if nsegs<=0
        points << newPoint(centrelat, centrelon, startangle, r);
        startangle +=  anglejump; nsegs-=1
    }
    return points;
end

def writeShapeFile(sfp, title, type, airspace_class, point, base, tops)
    return if(type=="GSEC")

    ibase=-1;
    ibase = 0 if base.StartsWith("SFC")
    ibase = base.to_i if base.index("ALT") || base.index("AGL")
    ibase = -(base[2..-1].to_i) if base.StartsWith("FL")

    puts "Can't interpret #{base}" if(ibase==-1)
    return if ibase<-55 || ibase>5000 # Base >5000ft or FL55

    #Check for things in the index which we want to represent with icons
    if(title.index(' RTM') || title.index('GVS') || title.index('LASER'))
        return # FIXME Write to sfp_points
    end
    type = "HIRTA" if(title.index('HIRTA'))
    type = "CTR" if type=="CTA/CTR" && airspace_class!="G" && ibase==0

    points = Array.new
    point_array = point.split(' ')
    while(point_array.size > 0)
        p = point_array.shift
        if(p=="CIRCLE")
            radius = point_array.shift.StartsWith("RADIUS=").to_f
            centrelat = parseLat(point_array.shift.StartsWith("CENTRE="))
            centrelon = parseLon(point_array.shift)
            #toplon, toplat = newPoint(centrelat, centrelon, 0, radius/3438.461)
            toplon = centrelon; toplat = centrelat + radius/60.0
            newpoints=curveToSegments(toplat, toplon, centrelat, centrelon,
                                      toplat, toplon, false)
            points << [toplon, toplat] # Complete the circle
            points.concat(newpoints)
            points << [toplon, toplat] # Complete the circle
            p=false
        elsif(p=="CLOCKWISE" || p=="ANTI-CLOCKWISE")
            radius = point_array.shift.StartsWith("RADIUS=").to_f
            centrelat = parseLat(point_array.shift.StartsWith("CENTRE="))
            centrelon = parseLon(point_array.shift)
            nextlats = point_array.shift.StartsWith("TO=")
            nextlons = point_array.shift
            nextlat = parseLat(nextlats)
            nextlon = parseLon(nextlons)

            newpoints = curveToSegments(points.last[1], points.last[0],
                                        centrelat, centrelon,
                                        nextlat,nextlon,
                                        p=="ANTI-CLOCKWISE")
            points.concat(newpoints)

            point_array.unshift(nextlons) #Put the TO= point back in for the
                                          #Final line
            p = nextlats
        end
        if p # We might have finished off above (eg. for a circle)
            #We have a straightforward point.
            lats = p
            lons = point_array.shift
            points << [parseLon(lons), parseLat(lats)] #Note order for shapefile
        end
    end
    #Check for non airspace stuff we want to deal with separately
    if(title.index('GVS') || title.index('LASER'))
        return # FIXME Write to sfp_points
    end

poly = ShapeLib::Polygon::new(*points)
    poly['name'] = title
    poly['type'] = type
    poly['class'] = airspace_class
    sfp.write poly
end

fp = File.open('/home/james/.MicroNotam/uk.air')

type=""; point=""; airspace_class=""; tops=""; base=""
title = nil
fp.each_line do |l|
    i = l.index('#')
    l = l[0..i] if i
    l.chop! # Removes either \n or #
    next if l.size==0

    a = l.StartsWith("TYPE="); type = a if a
    a = l.StartsWith("CLASS="); airspace_class = a if a
    a = l.StartsWith("TOPS="); tops = a if a
    a = l.StartsWith("BASE="); base = a if a
    a = l.StartsWith("POINT="); point += a + " " if a
    a = l.StartsWith("CIRCLE"); point += l + " " if a
    a = l.StartsWith("CLOCKWISE"); point += l + " " if a
    a = l.StartsWith("ANTI_CLOCKWISE"); point += l + " " if a
    a = l.StartsWith("TITLE=")
    if (a && title!=nil)
        sfp = sfp_main
        sfp = sfp_atz if type=="CTA/CTR" && airspace_class=="G"
        writeShapeFile(sfp, title, type, airspace_class, point, base, tops)
        point=""
    end
    title = a if a
end
sfp_main.close
sfp_atz.close