## makefile -- for M-750 air route web page
## Copyright       : http://www.fsf.org/copyleft/gpl.html
## Author          : Dan Jacobson -- http://jidanni.org/comm/air/m750/
## Created On      : March 2003
## Last Modified On: Thu Mar 23 09:10:36 2006
## Update Count    : 459

SANAS_Lon=121.692222
SANAS_Lat= 24.896944
PILOX_Lon=121.024722
PILOX_Lat= 24.415000
ANLOT_Lon=120.522500
ANLOT_Lat= 23.940833
TOROX_Lon=119.908056
TOROX_Lat= 23.353889
Dan_Lon=  120.866040
Dan_Lat=   24.181705
#Note we are careful to use all WGS84 datum data

#Related to north_pass.html
np: north_pass
	@./$<
	@echo; echo Fields:
	@sed -n '/print/,$${;s/.* #//p;}' $<
star:
	@echo Declination, right ascension, constellation, name, magnitude
	@awk 'BEGIN{FS=",";OFS="\t"};$$4~/^[567].:/&&\
	$$5<3{print $$4,$$3,$$1,$$5}' ~/.xephem/YBS-Altnames.edb|\
	awk -F: '$$1>=58&&$$1<=75'|sort -t: -n|column -t

# 2005: to do: remove usage of proprietary xephem program.
#This produces a straight long across xephem's whole sky view
T=/tmp/tmp.hzn.aid
short_horizon:basic
	seq 230 $$((360+32)) > $T;\
	set -- $$(echo -e \
	$(ANLOT_Lon)\\t$(ANLOT_Lat)\\n$(PILOX_Lon)\\t$(PILOX_Lat)|\
	project -T$(Dan_Lon)/$(Dan_Lat) -C0/90 -Fpq);\
	project -C$$1/$$2 -E$$3/$$4 -G0.02 -Lw|\
	awk 'BEGIN{OFMT="%.6f"};{print 360-$$1,$$2}'|sample1d -N$T > $T.1
	for height in 41 38 37 34 33 30 29; do let $$((height*=1000));\
	awk 'BEGIN{pif=45/atan2(1,1);\
	print "#M750 air route `horizon` at elev '$$height' ft above sea level"\
	;print "#From '$$PWD'"};\
	{print $$1%360,pif*atan2('$$height'*0.3048-777,(90-$$2)*60*1852)}' \
	$T.1 > ~/.xephem/m750-$${height}ft.hzn;done;rm $T $T.1
	#not a full circle, yes.
horizon:basic
	seq 220 $$((360+60)) > $T
	awk -F'[=#]' '/^[A-Z]*_Lon=/{a=$$2};/^[A-Z]*_Lat=/{print a,$$2};\
	/TOROX_Lat/{exit}' makefile|\
	sample1d -I.01|\
	project -T$(Dan_Lon)/$(Dan_Lat) -C0/90 -Fxypq|\
	awk 'BEGIN{pif=45/atan2(1,1)};{x=$$1;y=$$2;az=360-$$3;\
	dist=(90-$$4)*60*1852;h2="";for(h=20000;h<=45000;h+=1000){ \
	h2=h2 " " pif*atan2(h*.3048-dist*dist*.0000000675-777,dist)}};\
	{print az,dist,h2}'|sample1d -N$T|awk '{$$1=int($$1%360)};1'|\
	awk '{f=2;for(h=20000;h<=45000;h+=1000){\
	f++;fn="'$$HOME'/.xephem/m750-"h"ft.hzn";\
	print $$1,$$f > fn}};{printf "%d %.1f\n",$$1,$$2/1000}' >\
	azimuth_vs_distance
#offical M750 just needs 41 38 37 34 33 30 29 *1000
#675: earth's curvature

#Related to trail.html
how_far_north:basic
	echo $(Dan_Lon) >/tmp/$(Dan_Lon)
	project -C$(PILOX_Lon)/$(PILOX_Lat) -E$(ANLOT_Lon)/$(ANLOT_Lat) \
	-G0.05 -Lw|awk '{print $$1,$$2}'|sample1d -Fl -N/tmp/$(Dan_Lon)|\
	awk '{print ($$2-$(Dan_Lat))*60*1852}' #answer then printed in meters
	rm /tmp/$(Dan_Lon)
point_closest_to_dan:
	@echo 'Position(lon, lat), Distance(km):'
	@set -- $$(echo $(Dan_Lon) $(Dan_Lat)|\
	project -C$(PILOX_Lon)/$(PILOX_Lat) \
	-E$(ANLOT_Lon)/$(ANLOT_Lat) -Q -Frsq);\
	echo $$@;echo 'Bearing (forward, back), \
	another distance version, in meters:';\
	echo -e $(Dan_Lon) $(Dan_Lat)\\n$$1 $$2|\
	cs2cs -s +proj=latlong|grep -v pj_latlong|\
	awk '{a=a $$1" "$$2" "};END{print a}'|geod -I -p +ellps=WGS84
	#(Bearing is perpendicular to route at that point)
points_90_degrees_to_dan:basic
	echo $(Dan_Lon) >/tmp/$(Dan_Lon)
	project -C$(PILOX_Lon)/$(PILOX_Lat) \
	-E$(ANLOT_Lon)/$(ANLOT_Lat) -G0.05 -Lw|\
	awk '{print $$1,$$2}'|sample1d -Fl -N/tmp/$(Dan_Lon)|tee $@
	echo $(Dan_Lat) >/tmp/$(Dan_Lat)
	project -C$(PILOX_Lon)/$(PILOX_Lat) \
	-E$(ANLOT_Lon)/$(ANLOT_Lat) -G0.05 -Lw|\
	awk '{print $$1,$$2}'|sample1d -Fl -T1 -N/tmp/$(Dan_Lat)|tee -a $@
length_of_points_90_degrees_to_dan:points_90_degrees_to_dan
	set -- `cat $^`; project -C$$1/$$2 -E$$3/$$4 -G99 -Q
speed:#using 60 sec over 13 km as example
	python -c $$'k=13;sec=60;s=k/(sec/3600.0)\nfor \
	i in ("km/h",1),("knots",1.852),("mi/h",1.609):\
	print "%d %s;"%(s/i[1], i[0]),'
Sound_Speed=668.68#knots=1238.4km/h=0.344km/s at 20 deg C [should use colder]
#velocity = 331 m/s + (0.6 m/s/C)*Temperature OK, try 0 C someday
Plane_Speeds=range(150, 500, 50)#knots
see_hear_angle:#at moment when either is at nearest point to me on trail
	python -c "exec '''import math\nprint 'km/h, knots, angle(deg)'\
	\nfor i in $(Plane_Speeds):\
	print int(round(i*1.852)),i,\
	round(360/2/math.pi*math.asin(i/$(Sound_Speed)),1)'''"
basic:
	gmtset D_FORMAT = %.8f	#used haphazardly by me currently
arc_length:
	horiz=$$(echo $(Dan_Lon) $(Dan_Lat)|\
	project -C$(PILOX_Lon)/$(PILOX_Lat) \
	-E$(ANLOT_Lon)/$(ANLOT_Lat) -Q -Fq);\
	python -c "exec '''ft2km=3280.84;h2=$$horiz**2\nimport\
	math\nprint 'height(ft, km), length(arc minutes)'\nfor \
	i in range(25000,41000,5000):\
	print '%d %4.1f %.1f'\
	% (i,i/ft2km,360*60*0.070/(2*math.pi*math.sqrt((i/ft2km)**2+h2)))'''"

## 2005.12.19: to plot on paper maps see
## http://jidanni.org/comm/air/routes stuff instead!

#and duh, to84c isn't defined, etc.
#plot air route M-750
p=+proj=tmerc +lon_0=121 +x_0=250000 +k=0.9999 -f %.3f
# m750pa: /tmp/notpa
# 	project -C$(PILOX_Lon)/$(PILOX_Lat) \
# 	-E$(ANLOT_Lon)/$(ANLOT_Lat) -G0.05 -Lw|\
# 	awk '{print $$1,$$2}'|\
# 	cs2cs +proj=latlong +datum=WGS84 +to $p \
# 	+ellps=aust_SA +towgs84=$(to84c)|\
# 	sample1d -Fl -N$?|perl -walne 'printf "%.0f %.0f\n", $$F[0],$$F[1]'
# /tmp/notpa:
# 	seq 202000 2000 250000 > $@
# m750at: /tmp/notat
# 	project -C$(ANLOT_Lon)/$(ANLOT_Lat) \
# 	-E$(TOROX_Lon)/$(TOROX_Lat) -G0.05 -Lw|\
# 	awk '{print $$1,$$2}'|\
# 	cs2cs +proj=latlong +datum=WGS84 +to $p \
# 	+ellps=aust_SA +towgs84=$(to84c)|\
# 	sample1d -Fl -N$?|perl -walne 'printf "%.0f %.0f\n", $$F[0],$$F[1]'
# /tmp/notat:
# 	seq 160000 2000 201000 > $@
#Will other routes come near to me?
#http://www.vatroc.org/chinese/info/aipnav.htm
#might be TWD67's, using instead pictures on
#http://www.vatroc.org/chinese/teach/rcssrckh.html
TONPU_Lon=120.795#47.7
TONPU_Lat=24.221667#13.3
XEROX_Lon=121.103333# 06.2
XEROX_Lat=24.628333#37.7
#jeppesen:
DALIN_Lon=120.28#22.8
DALIN_Lat=23.593333# 35.6
XEROX-DALIN: /tmp/not
	gmtset D_FORMAT = %.8f
	project -C$(XEROX_Lon)/$(XEROX_Lat) \
	-E$(DALIN_Lon)/$(DALIN_Lat) -G0.05 -Lw|\
	awk '{print $$1,$$2}'|\
	cs2cs +proj=latlong +datum=WGS84 +to $p \
	+ellps=aust_SA +towgs84=$(to84c)|\
	sample1d -Fl -N$?|perl -walne 'printf "%.0f %.0f\n", $$F[0],$$F[1]'
TONPU-DALIN: /tmp/not
	gmtset D_FORMAT = %.8f
	project -C$(TONPU_Lon)/$(TONPU_Lat) \
	-E$(DALIN_Lon)/$(DALIN_Lat) -G0.05 -Lw|\
	awk '{print $$1,$$2}'|\
	cs2cs +proj=latlong +datum=WGS84 +to $p \
	+ellps=aust_SA +towgs84=$(to84c)|\
	sample1d -Fl -N$?|perl -walne 'printf "%.0f %.0f\n", $$F[0],$$F[1]'
#both go past xinshe...
tonpu:
	echo $(TONPU_Lon)E $(TONPU_Lat)N 0|\
	cs2cs +proj=latlong +datum=WGS84 +to $p +ellps=aust_SA +towgs84=$(to84c)
#it seems my tonpu isn't very accurate.