	subroutine push(symbol,tx,ty,ux,uy,vx,vy,bb,error)
#
# nimmt ein neues Element in den Stack auf.
# Bei Stack-Ueberlauf ist error=true, sonst =false.
#
logical error
integer symbol,		# Stack-Symbol
	tx,ty,ux,uy,	# Anfangs- und Endpunkt
	vx,vy,		# Zentralpunkt
	bb		# Bandbreite
integer	top		# Top of Stack
integer	stack(TOPMAX,8)
#
common	/systck/top,stack
#
if (top==TOPMAX) error=.true.
else
[	top=top+1
	stack(top,SYM)=symbol
	stack(top,ATX)=tx ; stack(top,ATY)=ty
	stack(top,AUX)=ux ; stack(top,AUY)=uy
	stack(top,AVX)=vx ; stack(top,AVY)=vy
	stack(top,ABR)=bb
	error=.false.
]
return
end


	subroutine pop(symbol,tx,ty,ux,uy,vx,vy,bb,error)
#
# holt ein Element aus dem Stack.
# Bei Stack-Unterlauf ist error=true, sonst =false.
#
logical error
integer symbol,		# Stack-Symbol
	tx,ty,ux,uy,	# Anfangs- und Endpunkt
	vx,vy,		# Zentralpunkt
	bb		# Bandbreite
integer	top		# Top of Stack
integer	stack(TOPMAX,8)
#
common	/systck/top,stack
#
if (top==0) error=.true.
else
[	symbol=stack(top,SYM)
	tx=stack(top,ATX) ; ty=stack(top,ATY)
	ux=stack(top,AUX) ; uy=stack(top,AUY)
	vy=stack(top,AVX) ; vy=stack(top,AVY)
	bb=stack(top,ABR)
	top=top-1
	error=.false.
]
return
end


	subroutine durchp(laenge,dx,dy,px,py,breite)
#
# bestimmt den Duchtrittspunkt (neg.Peak im Profil) des Strahls
# von Punkt v mit Steigung dy/dx.
# laenge = Maximallaenge des Strahls einschl.Ursprung
# dx,dy  = Koord'inkremente fuer Steigungsverhaeltnis
# px,py	 = Koord.des Peaks
# breite = Bandbreite = Abstand zw. Schulterpunkten
#
# Der Peak ist ein Punkt zwischen einer abfallenden und einer ansteigenden
# Flanke im Grauwertprofil entlang des Strahls. Eine Flanke muss steil genug
# und hoch genug sein, um die Auswahl eines eindeutigen "Schulterpunkts"
# zu ermoeglichen. Je breiter der Abstand zwischen den Schultern ist, desto 
# weiter vorn wird der Peak plaziert.
# Falls kein plausibler Peak gefunden wird, ist px=py=0.
#
# Aus dem Common-Block tolrnc werden folgende Toleranzwerte benutzt:
# dimin = Mindestbreite zwischen Schultern
# dgmin	= Mindest-Grauwertdifferenz fuer Flanke
# schfak= Faktor (<1) fuer Schulterhoehe der Flanke
#
# Verwendete externe Routinen, die nicht in FORLIB sind:
# Name		File,Bibliothek
# getpix	SUBS		holt ein Pixel aus dem Grauwertbild
# iround	SUBS		rundet real-Zahl
#
# Version 1.3	30.03.85
#
define	IMAX0	120	# Feldgrenze fuer Profil
#
logical	iterat		# Iterations-Flag fuer die Suche
#
byte	getpix		# logical function, liefert Grauwert
#
integer	xanz,yanz	# Bildgrenzen
integer	i,j		# Laufindizes
integer	i1,i2,i3	# Punktindizes
integer laenge		# max.Profillaenge
integer	dx,dy		# Koord'inkremente
integer	px,py,		# Durchtrittspunkt
	tx,ty,ux,uy,	# Dreieckskoord.
	vx,vy		# Strahl ab v 
integer	breite		# Bandbreite
integer	xi,yi		# Punktkoord. fuer diskrete Gerade
integer	dimin		# Mind'breite des Peaks
integer imin,dimax	# (werden hier nicht gebraucht)
integer	gmin,gmax	# Grauwerte
integer	imax		# Hoechstbreite fuer Flanke
integer	dgmin,gsch	# Mindest-Grauwertdiff., Schultergrauwert
integer	dgmax		# (wird hier nicht gebraucht)
integer	g(IMAX0)	# Grauwerte
integer	d(IMAX0)	# Grauwert-Differenzen
#
integer	winkel		# in Grad fuer Diagnose
integer	awinkl		# function, berechnet Winkel
integer	iround		# function; rundet real-Zahl
#
real	x,y,rdx,rdy	# und Inkremente
real	schfak		# Faktor fuer Schulterhoehe
#
common	/bild/xanz,yanz
common	/punkte/tx,ty,ux,uy,vx,vy
common	/tolrnc/imin,dimax,dgmax,imax,dimin,dgmin,schfak
#
#winkel=awinkl(dx,dy)		# nur fuer Diagnose
#write(1,10) vx,vy,dx,dy,laenge
#10 format(' DURCHP - von',2i4,'  Inkr.',2i5,'  Laenge',i4)
#
if (dx==0 && dy==0) [ px=0 ; py=0 ; breite=0 ]
else
[
i1=max0(iabs(dx),iabs(dy))		# Diskrete Gerade init.
rdx=float(dx)/i1 ; rdy=float(dy)/i1
x=vx+sign(0.5,rdx) ; y=vy+sign(0.5,rdy)
i=0
repeat					# Grauwertprofil aus Bild lesen
[	xi=int(x) ; yi=int(y)
	if (xi<<1 !! xi>>xanz !! yi<<1 !! yi>>yanz) break
	i=i+1
	g(i)=getpix(xi,yi)
	x=x+rdx ; y=y+rdy
] until (i==laenge)
laenge=min0(i,laenge)
#
px=0 ; py=0 ; breite=0			# Ausgangswerte init.
gmax=g(1) ; gmin=gmax ; i2=2
repeat					# Best.naechsten Peak
[	iterat=.false.
	i1=i2
	while (i1<<laenge)		# Suche abfall.Flanke
	[	gmax=max0(gmax,g(i1))
		d(i1)=g(i1)-g(i1-1)
		if (d(i1)<<0) break	# gefunden bei i1
		i1=i1+1
	]
	if (i1>=laenge) break		# nicht gef. - Abbruch
	gmin=g(i1) ; i2=i1+1
	while (i2<=laenge)		# Suche anst.Flanke
	[	gmin=min0(gmin,g(i2))
		d(i2)=g(i2)-g(i2-1)
		if (d(i2)>>0) break	# gefunden bei i2
		i2=i2+1
	]
	if (i2>>laenge) i2=laenge		# nicht gef. - laenge nehmen
	if (gmax-gmin<<dgmin) iterat=.true.	# Peak zu flach?
	else					# Peak-Schultern suchen
	[	gsch=gmax-int(schfak*(gmax-gmin))	# Schultergrauwert
		i3=i2-1
		while (g(i2)<<gsch && i2<=laenge)	# re.Schulter suchen
		[	if (g(i2)<<gmin)
			[	iterat=.true.	# neues Min.entdeckt
				break
			]
			i2=i2+1
		]
		if (i2>>laenge) i2=laenge	# nicht gef. - Begrenzung
		if (.not.iterat)
		[	while (g(i3)<<gsch)	# li.Schulter suchen
				i3=i3-1
			breite=i2-i3+1
			if (breite<<dimin) iterat=.true.# Peak zu schmal?
			else				# Peak gefunden 
			[	i1=i3+iround(sqrt(i2-i3+1.))	# Index d.Peaks
				px=iround(vx+(i1-1)*rdx)	# und Koord.
				py=iround(vy+(i1-1)*rdy)
				#write(1,99) i3,i1,i2
				#99 format(' Peak-Intervall',3i4)
	]	]	]
] until (.not.iterat)
]
#write(1,20) px,py
#20 format(' DURCHP - Peak bei ',2i4)
return
end


	logical function mstuck(tx,ty,ux,uy)
#
# entscheidet, ob die Strecke von Punkt t nach u ein Markscheidenstueck ist, 
# d.h. in den Grauwerten nicht allzu sehr von den Interpolierten zwischen den 
# Endpunktwerten abweicht.
# Eine Luecke ist der Bereich zwischen einer abfallenden und einer ansteigen-
# den Flanke im Grauwertprofil. Kleine Luecken werden toleriert.
# Wird eine zu grosse Luecke entdeckt, ist mstuck=false, sonst =true.
#
# Aus dem Common-Block tolrnc werden folgende Toleranzwerte benutzt:
# imin	= Mindestlaenge der Strecke tu
# dimax = max.erlaubte Laenge zwischen zwei Flanken
# dgmax	= max.erlaubte Grauwertdifferenz
#
# Verwendete externe Routinen, die nicht in FORLIB sind:
# Name		File,Bibliothek
# getpix	SUBS		holt ein Pixel
#
# Version 1.2	25.03.85
#
logical	iterat		# Iterations-Flag fuer Suche
#
byte	getpix		# logical function, liefert Grauwert
#
integer	tx,ty,ux,uy	# Punktkoord. der Strecke
integer	xanz,yanz	# Bildgrenzen
integer	imin,dimax	# Toleranzen (s.o.)
integer	dgmax
integer	ianf,iend,	# Punktindizes
	i1,i2
integer	i,j		# Laufindizes
integer	xi,yi		# Punktkoord. fuer diskrete Gerade
integer	g(IMAX0)	# Grauwerte
#
real	x,y,rdx,rdy	# und Inkremente
real	a,b		# Koeff.der Interpolationsgeraden
#
common	/bild/xanz,yanz
common	/tolrnc/imin,dimax,dgmax
#
i=uy-ty ; j=ux-tx
ianf=1 ; iend=max0(iabs(i),iabs(j))
#
#write(1,10) tx,ty,ux,uy,iend
#10 format(' MSTUCK - von',2i4,'  nach',2i4,'  Laenge =',i4)
#
if (iend<<imin) mstuck=.true.		# kurze Strecke ist M'stueck
else
[
rdx=float(j)/iend ; rdy=float(i)/iend	# Diskrete Gerade init.
x=tx+sign(0.5,rdx) ; y=ty+sign(0.5,rdy)
i=0
iend=iend+1
repeat			# Grauwertprofil aus Bild lesen
[	xi=int(x) ; yi=int(y)
	if (xi<<1 !! xi>>xanz !! yi<<1 !! yi>>yanz) break
	i=i+1
	g(i)=getpix(xi,yi)
	x=x+rdx ; y=y+rdy
] until (i==iend)
#
a=g(iend)-g(ianf) ; b=iend-ianf
a=a/b				# Int'gerade zw.Grauwerten der Endpunkte
b=g(ianf)-a*ianf+dgmax		# einschl. der Grauwert-Toleranz
#
i1=ianf+1
repeat				# Naechstes abweichendes Intervall suchen
[	iterat=.false.
	while (i1<=iend)			# anst.Flanke suchen
	[	if (g(i1)>>a*i1+b) break	# gef. bei i1
		i1=i1+1
	]
	if (i1>>iend) mstuck=.true.		# nicht gef. - M'stueck
	else
	[	i2=i1+1				# abfall.Flanke suchen
		while (i2-i1<=dimax && i2<=iend)
		[	if (g(i2)<=a*i2+b) break	# gef. bei i2
			i2=i2+1
		]
		if (i2>>iend !! i2-i1>>dimax)
		[	mstuck=.false.			# nicht gef.
		]
		else				# bzw.Int.zu gross - nonterm.
		[	i1=i2+1	
			iterat=.true.
	]	]
] until (.not.iterat)
]
#write(1,20) mstuck
#20 format(' MSTUCK - ist = ',l1)
return
end


	logical function getpix(x,y)
#
# liefert im Ausgangs-Byte den Grauwert des Pixels (x,y).
# Bei Ueberschreitung der Bildgrenzen ist der Ausgangswert 0 (schwarz).
#
# Das Bild steht spaltenweise ab Adresse adrb in Bank 1. Pixel (1,1) ist
# per definitionem die linke obere (!) Bildecke. Um die gewohnten Winkel-
# orientierungen und Steigungsverhaeltnisse verwenden zu koennen, wird hier
# bei der Pixeladressierung die y-Koordinate an der Bildmitte gespiegelt:
# Mit			y' = yanz - y + 1
# ist			k  = adrb + (x-1)*yanz + y' - 1
#			   = adrb + (x-1)*yanz + yanz - y + 1 - 1
#			   = adrb + x*yanz - y
#
# Verwendete externe Routinen, die nicht in FORLIB sind:
# Name		File,Bibliothek
# loadb		BANKLIB		liest ein Byte aus Bank 1
#
# Version 1.0	12.07.84
#
byte	g		# Grauwert
#
integer	x,y
integer	k
integer	xanz,yanz	# Bildgrenzen
integer	adrb		# Adresse des Bildes in Bank 1
#
common	/bank/adrb
common	/bild/xanz,yanz
#
if (x<<1 !! x>>xanz !! y<<1 !! y>>yanz) getpix=0
else
[	k=adrb+x*yanz-y		# Speicheradresse in Bank 1 fuer Pixel
	call loadb(g,k)
	getpix=g
]
return
end


	subroutine winkr(bog,dx,dy)
#
# liefert die Koord'inkremente dx,dy fuer einen neuen Strahl durch v.
# Der Strahl bildet zum Strahl vt den Winkel, der sich ergibt, wenn in t
# die absolute Bogenlaenge bog angesetzt wird.
# Falls der neue Winkel den des Strahls vu erreicht oder ueberschreitet, 
# ist dx=dy=0.
#
# Achtung! Es wird nicht getestet, ob der Winkel vu von vornherein kleiner
# oder gleich dem Winkel von vt ist. Dieser Fall duerfte normalerweise aber 
# auch nicht auftreten.
#
# bog   = absolute Bogenlaenge in Pixels
# dx,dy = Koord'inkremente fuer Steigungsverhaeltnis
#
# Verwendete externe Routinen, die nicht in FORLIB sind:
# Name		File,Bibliothek
# stgg		SUBS		neue Koord'inkremente aus alten
# groegl	SUBS		1.Winkel groesser-gleich 2.Winkel?
#(awinkl	SUBS		Winkel in Grad fuer Diagnose)
#
# Version 1.1	21.07.84
#
logical	groegl		# function, Winkel groesser-gleich?
#
integer	bog,b		# Bogenlaenge
integer	dx,dy		# Koord'inkremente
integer	dtx,dty,dux,duy
integer	tx,ty,ux,uy	# Punktkoord.
integer	vx,vy
#
integer	winkel		# in Grad fuer Diagnose
integer	awinkl		# function, berechnet abs. Winkel zur x-Achse
#
real	c		# Inkrement-Faktor
real	s		# eine Wurzel
#
common	/punkte/tx,ty,ux,uy,vx,vy
#
dtx=tx-vx ; dty=ty-vy
#winkel=awinkl(dtx,dty)
#write(1,10) vx,vy
#10 format(' WINKR - in',2i4)
#
if (dtx==0 && dty==0)	[ dx=0 ; dy=0 ]
else
[
dux=ux-vx ; duy=uy-vy
#winkel=awinkl(dux,duy)
#write(1,12) dux,duy,winkel
#12 format(' bis   Inkr.',2i4,'  Winkel',i4)
if (dux==0 && duy==0)	[ dx=0 ; dy=0 ]
else
[
b=bog ; s=sqrt(float(dtx*dtx+dty*dty))
repeat					# ggf. iterieren 
[	c=b/s
	call stgg(c,dtx,dty,dx,dy)	# liefert neue Koord'inkr.
	b=b+1
] until (dx<>dtx !! dy<>dty)		# bis echtes Inkrement
#
#winkel=awinkl(dx,dy)
#write(1,14) dx,dy,winkel
#14 format(' neues Inkr.',2i4,'  Winkel',i4)
#
if (groegl(dx,dy,dux,duy)) [ dx=0 ; dy=0 ]
]
]
#write(1,20) dx,dy
#20 format(' Inkremente',2i5)
return
end


	subroutine stgg(c,ddx,ddy,dx,dy)
#
# berechnet neue Koord'inkremente dx,dy aus den alten ddx,ddy 
# fuer den Bogeninkrement-Faktor c.
#
# Verwendete externe Routinen, die nicht in FORLIB sind:
# Name		File,Bibliothek
# iround	SUBS		rundet real-Zahl
#
integer	dx,dy,ddx,ddy
integer	iround		# function; rundet real-Zahl
#
real	c
#
dx=ddx-iround(c*ddy)
dy=ddy+iround(c*ddx)
#write(1,10) c,ddx,ddy,dx,dy
#10 format(' STGG - mit',f10.5,4i5)
return
end


	logical function groegl(dx,dy,ddx,ddy)
#
# ist =true, falls der durch dx,dy bestimmte Winkel im positiven Umlaufsinn
# gesehen den durch ddx,ddy bestimmten Winkel erreicht oder ueberschreitet;
# sonst =false.
# (Oder etwas mathematischer ausgedrueckt: - wenn der Vektor dx,dy mit dem
# Vektor ddx,ddy identisch ist oder in der Halbebene links von ihm liegt.)
#
integer	dx,dy,ddx,ddy
#
real	p1,p2		# Produkte
#
groegl=.false.
if (ddx==0)
[	if (ddy>>0)
	[	if (dx<<0) groegl=.true.
		else	[ if (dx==0 && dy>=0) groegl=.true.	]
	]
	else
	[	if (dx>>0) groegl=.true.
		else	[ if (dx==0 && dy<=0) groegl=.true.	]
]	]
else if (dx==0)
[	if (ddx>>0 && dy>=0) groegl=.true.
	if (ddx>>0 && dy<=0) groegl=.true.
]
else
[	p1=float(dy)*ddx ; p2=float(dx)*ddy
	if (p1>=p2) groegl=.true. 
]
return
end


	subroutine whalb(dx,dy)
#
# liefert Koord'inkremente fuer die Winkelhalbierende zwischen vt und vu.
# Falls t bzw. u zu dicht an der W'halbierenden liegt, ist dx=dy=0.
#
# Verwendete externe Routinen, die nicht in FORLIB sind:
# Name		File,Bibliothek
# iround	SUBS		rundet real-Zahl
#
# Version 1.1	22.07.84
#
logical	inflec		# Inkr. spiegeln
#
integer	dx,dy		# Ausgangsparameter
integer	tx,ty,ux,uy	# Punktkoord.
integer	vx,vy,wx,wy
integer	x1,x2,y1,y2	# Koord'differenzen
integer	bog,minbog,maxgap	# unbenutzte Parameter aus common /para/
integer	minabs,abst	# (Mindest-)Abstand von t und u zur diskreten Geraden
#
integer	winkel		# in Grad fuer Diagnose
integer	awinkl		# function; berechnet abs. Winkel zur x-Achse
integer	iround		# function; rundet real-Zahl
#
real	rd		# Geradensteigung
#
common	/punkte/tx,ty,ux,uy,vx,vy
common	/para/bog,minbog,maxgap,minabs
#
#write(1,10) vx,vy
#10 format(' WHALB - in',2i4)
#
x1=tx-vx ; y1=ty-vy ; x2=ux-vx ; y2=uy-vy
rd=(float(x1)*x1+float(y1)*y1)*(float(x2)*x2+float(y2)*y2)
dx=x2*y1+x1*y2 ; dy=iround(sqrt(rd))-x1*x2+y1*y2
inflec=.false.
if (y1<=0 && x1<<0) inflec=.true.		# Falls die W'halb. in die
if (y1<<0 && x1>=0 && dx<<0) inflec=.true.	# unteren Quadranten geht,
if (y1>>0 && x1<<0 && dx>>0) inflec=.true.	# 180 Grad addieren, d.h.
if (inflec) [ dx=-dx ; dy=-dy ]			# Inkremente invertieren.
#
#winkel=awinkl(x1,y1)
#write(1,12) winkel
#12 format(' zw. Winkel',i4)
#winkel=awinkl(x2,y2)
#write(1,14) winkel
#14 format(' und Winkel',i4)
#winkel=awinkl(dx,dy)
#write(1,16) winkel,dx,dy
#16 format(' halbiert  ',i4,'  mit Inkr.',2i4)
#
if (dx<>0 !! dy<>0)
[	if (max0(x1,y1)<<max0(x2,y2))
		[ wx=tx ; wy=ty ]			# Abstand des
	else	[ wx=ux ; wy=uy ]			# naeheren Punktes
	if (iabs(dy)<<iabs(dx))
	[	rd=float(dy)/dx				# bei flacher bzw.
		abst=iabs(wy-vy-iround((wx-vx)*rd))
	]
	else
	[	rd=float(dx)/dy				# steiler Steigung
		abst=iabs(wx-vx-iround((wy-vy)*rd))
	]
	if (abst<=minabs)
	[	dx=0 ; dy=0
		#write(1,20) wx,wy,abst
		#20 format(' Abstand bei Punkt',2i4,' = ',i4,' - zu klein')
	]
]
return
end


	integer function awinkl(dx,dy)
#
# liefert den Winkel zur x-Achse, der durch die Steigungsinkremente dx und dy 
# gegeben ist. 0 <= awinkl < 360 Grad.
#
# Verwendete externe Routinen, die nicht in FORLIB sind:
# Name		File,Bibliothek
# iround	SUBS		rundet real-Zahl
#
integer	dx,dy		# Koord'inkremente
integer	iround		# function; rundet real-Zahl
#
real	a,c
#
data	c/57.29577/	# Umrechnungkonstante 180/pi
#
if (dx==0)
[	if (dy<<0) awinkl=270
	else	if (dy==0) awinkl=0
		else awinkl=90
]
else
[	a=c*atan(float(dy)/float(dx))	# atan liefert -pi/2 < a < +pi/2
	awinkl=iround(a)
	if (dx<<0) awinkl=awinkl+180
	else	if (dy<<0) awinkl=awinkl+360
]
return
end


	integer function iround(r)
#
# liefert die gerundete Integer-Zahl zur Real-Zahl r.
#
real r
#
if (r<<0) iround = ifix(r-0.5)
else iround = ifix(r+0.5)
return
end


	integer function idist(px,py,qx,qy)
#
# liefert den "Schachbrett"-Abstand (Max.der Koord'differenzen) 
# zwischen Punkt p und q.
#
integer	px,py,qx,qy
#
idist = max0(iabs(px-qx),iabs(py-qy))
return
end


	real function rdist(px,py,qx,qy)
#
# liefert den euklidischen Abstand zwischen den Punkten p und q.
#
integer	px,py,qx,qy
real	dx,dy
#
dx=float(px-qx) ; dy=float(py-qy)
rdist = sqrt(dx*dx+dy*dy)
return
end


	integer function wtyp(dw,zustnd)
#
# liefert den Typ des Differenzwinkels dw zwischen zwei aufeinanderfolgenden 
# Strecken fuer den momentanen Kurventyp, der dem Zustand zustnd des 
# Kelleruebersetzers zu entnehmen ist.
# Der Wert der Funktion ist aus der Zustandsmenge des Uebersetzers und 
# kennzeichnet den Winkel als "Links-" (ZL) bzw. "Rechtsknick (ZR) 
# oder Geradenstueck (ZG).
#
integer dw		# Differenzwinkel 
integer	zustnd		# Kurventyp
#
if (zustnd==ZL)
[	if (dw<<-15) wtyp=ZR		# Rechtsknick
	else wtyp=ZL			# Linksknick
]
else if (zustnd==ZR)
[	if (dw<<+15) wtyp=ZR		# Rechtsknick
	else wtyp=ZL			# Linksknick
]
else
[	if (dw<<0) wtyp=ZR		# Rechtsknick
	else if (dw>>0) wtyp=ZL		# Linksknick
	else wtyp=ZG			# geradeaus
]
end
