# duplicates the first part of xsections.aml
# This version reads a profile and builds xsection lines.
# It builds a .gen file, and also puts the info into draw.aml for
# quicker examination.

import sys, string, math
print "xsections.py, draft version .9, hgreen@u.washington.edu"
interval =  10000	# meters between sections
baseline =  15000	# meters of half the baseline for establishing the angle
halfwidth = 50000	# meters of xsection width
# equatorial radius 6378137, polar = 6356752.3, Interpolate for 30 deg
radian = 180./math.pi
degmeters = 6371009 * math.pi / 180.0  # degree on earth ~= 30 deg N
f=open('1_prf.xy','r')
g=open('sections.gen','w')
draw=open('draw.aml','w')
lines = f.readlines()
idxy=[]
# calculate distance in meters, while leaving x&y geographic
for i in range(len(lines)-1,-1,-1):
  rec=string.split(lines[i])
  x=float(rec[4])
  y=float(rec[5])
  if i==len(lines)-1:	# bottom point
    dist=0.
  else:
    dist+=degmeters*((y-oldy)**2 + ((x-oldx)*math.cos(y/radian))**2)**0.5
  idxy.append([dist,x,y])
  oldx=x
  oldy=y

# As halfbase may be larger or smaller than interval, make 3 passes through the list
first=[]	# first baseline point
d=0.
i=0
sought=0.
while i < len(lines):
  while d < sought:
    d=idxy[i][0]
    i+=1
    if i==len(lines):
      break
  if i!=len(lines):
    first.append([d,idxy[i][1],idxy[i][2]])
    sought+=interval

second=[]       # second baseline point
d=0.
i=1
sought=baseline
while i < len(lines):
  while d < sought:
    d=idxy[i][0]
    i+=1
    if i==len(lines):
      break
  if i!=len(lines):
    second.append([d,idxy[i][1],idxy[i][2]])
    sought+=interval

cross=[]	# center of Xsection
d=0.
i=1
sought=baseline/2.
while i < len(lines):
  while d < sought:
    d=idxy[i][0]
    i+=1
    if i==len(lines):
      break
  if i!=len(lines):
    cross.append([d,idxy[i][1],idxy[i][2]])
    sought+=interval

for i in range(len(second)):
  dx=(second[i][1]-first[i][1]) * math.cos(cross[i][2]/radian)
  dy=second[i][2]-first[i][2]
  angle=math.atan2(dx,dy)+math.pi/2.
  dx2=math.sin(angle)*halfwidth/(degmeters*math.cos(cross[i][2]/radian))
  dy2=math.cos(angle)*halfwidth/degmeters
  g.write("%i\n%.5f %.5f\n%.5f %.5f\nend\n" % (i,cross[i][1]+dx2,cross[i][2]+dy2,cross[i][1]-dx2,cross[i][2]-dy2))
  draw.write("line %.5f %.5f %.5f %.5f\n" % (cross[i][1]+dx2,cross[i][2]+dy2,cross[i][1]-dx2,cross[i][2]-dy2))
g.write("end\n")
g.close()
draw.close()