program triangle
c	find triangulations given #vertices and Euler characteristic.
	implicit none
	integer Edges, simp, Nedges, NF, NV, NE, from, vert, edge,
     .		Vertices, Nvertices, V,E,F, Nb_edges, orient,
     .		first_edge, set_connect  
	logical dupl_edges, want_orientable
	common / / simp (100, 3), Edges (100, 175, 2), Nedges(100), 
     .		from(100), NF, NV, NE, vert, edge, dupl_edges,
     .		Vertices (100,100), Nvertices (100), V,E,F, Nb_edges,
     .		want_orientable, orient (100), first_edge (100,2),
     .		set_connect 

	character*1 choice
	integer EC, I(100), J(100), K, L
	integer cnt

c	function definitions
	logical checkface, trim_tree, orientable, sym_check

	vert = 1		!tags that keep track of your path in the tree.
	edge = 2

	write (6,*) 'Number of vertices?'
	read (5,*) V
	write (6,*) 'Euler characteristic?'
	read (5,*) EC
	write (6,*) 'Orientable?  Y or N'
	read (5,'(a)') choice
	if (choice.eq.'Y'.or.choice.eq.'y') then
	    want_orientable = .true.
	    write (6,*) 'Will look for an orientable manifold'
	    write (6,*)
	else
	    want_orientable = .false.
	    write (6,*) 'Will look for a non-orientable manifold'
	    write (6,*)
	end if
	write (6,*) 'Fixed # vertices connected to a given vertex?'
	write (6,*) '    0 means no check for this.'
	read (5,*) set_connect

	write (6,*) 
	
	if (V.le.EC) stop 'No go'
	if (V.lt.4) stop 'No go'

	F = (V-EC)*2
	E = (V-EC)*3

	
	NF = 1		!also is level of recursion
	NV = 3
	simp (1,1) = 1
	simp (1,2) = 2
	simp (1,3) = 3
	
 10	continue
c	write (7,*) 'Simplex ', cnt+1
c	if (cnt+1.eq.1999) then
c	    write (7,*) ((simp(K,L), L = 1,3), K = 1, NF-1)
c	    pause
c	end if
	call boundary 				!sets Edges
c						!array & finds #edges for
c						!Euler char.

	cnt = cnt + 1
	if (mod(cnt, 1 000 000).eq.0) write (6,*) 'cnt is ', cnt

c	write (7,*) 'nf,ne,nv,nedges(nf), nvertices(nf), nb_edges'
c	write (7,*) nf,ne,nv,nedges(nf), nvertices(nf), nb_edges
c	write (7,'(a, 2i3)') 'first_edge ', (first_edge (NF, K), K = 1, 2)
c	write (7,'(a,20(i3))') '   Simplex', (simp(nf,K), K=1,3)
c	write (7,'(a,20(2i3,2x))') 
c     .	 '   Edges', ((edges(nf, K,L), L=1,2), K=1, nedges(nf))
c	write (7,'(a,20(i3))') 
c     .	'   Vertices', (Vertices(nf, K), K=1, nvertices(nf))
c	write (7,*)

c   Does the structure violate symmetry conditions?
	if (.not.sym_check()) then
	    NF = NF - 1
	    if (from(NF+1).eq.Vert) NV = NV - 1
	    if (NF.eq.0) then
		write (6,*) cnt, ' nodes of the tree were looked at'
		stop 'No triangulation'
	    end if
	    goto 20			!back in recursion
	end if

	if ((.not.dupl_edges).and.
     .	    NF.eq.F.and.NE.eq.E.and.NV.eq.V.and.Nedges(NF).eq.0) then
	    if (.not.want_orientable) then
		orient(1) = 1
		do K = 2, NF
		    orient (K) = 0
		end do
		if (orientable()) then
		    NF = NF - 1
		    if (from(NF+1).eq.Vert) NV = NV - 1
		    goto 20		!back in recursion
		end if
	    end if
	    write (6,*) 'Triangulation:'
	    do K = 1, NF
		write (6,'(3i3)') (simp (K, L), L= 1, 3)
	    end do
	    write ( 6,*) cnt, ' nodes of the tree were looked at'
	    stop
	end if
c
	if (trim_tree()) then
	    NF = NF - 1
	    if (from(NF+1).eq.Vert) NV = NV - 1
	    if (NF.eq.0) then
		write (6,*) cnt, ' nodes of the tree were looked at'
		stop 'No triangulation'
	    end if
	    goto 20	!back one recursion level
	end if
c
c	Loop through edges list and vertices list, make new simplexes.
	    J(NF) = 0			!loop index on vertices to connect to
 21	    continue			!top of loop
	    J(NF) = J(NF) + 1
	    if (J(NF).gt.Nvertices(NF)) goto 29	!out of loop
	    if (checkface (first_Edge(NF, 1),
     .			   first_EDGE (NF, 2), 
     .			   Vertices(NF, J(NF)))) goto 10	
c						!down one level of recursion
c	    checkface makes new simplex if everything is OK.
 20	    continue		!returning from recursion
	    goto 21		!bottom of J(NF) loop
 29	    continue		!exit from J(NF) loop
	    
	if (NF.gt.1) then		!go back to node of tree you came from
	    NF = NF - 1
	    if (from(NF+1).eq.Vert) NV = NV - 1
	    goto 20	!back one recursion level
	end if
c
	write (6,*) 'No triangulation'
	write (6,*) cnt, ' nodes of the tree were looked at'
	stop
	end
c
	subroutine boundary
	implicit none
	integer Edges, simp, Nedges, NF, NV, NE, from, vert, edge,
     .		Vertices, Nvertices, V,E,F, Nb_edges, orient,
     .		first_edge, set_connect  
	logical dupl_edges, want_orientable
	common / / simp (100, 3), Edges (100, 175, 2), Nedges(100), 
     .		from(100), NF, NV, NE, vert, edge, dupl_edges,
     .		Vertices (100,100), Nvertices (100), V,E,F, Nb_edges,
     .		want_orientable, orient (100), first_edge (100,2),
     .		set_connect 

	integer t1, t2, t3, I, J, code_ij

	integer code_edge (0:175*175)

	dupl_edges = .false.
	NE = 0
	Nb_edges = 0

	do I = 1, NV*NV
	    code_edge(I) = 0
	end do

	do I = 1, NF
	    t1 = simp(I,1) -1 + NV*(simp(I,2)-1)		!store edges in 
								!base NV code.
	    code_edge ( t1) = code_edge (t1) + 1	!simplexes have vertices
c								!in ascending order

	    t2 = simp(I,1) -1 + NV*(simp(I,3)-1)
	    code_edge (t2) = code_edge (t2) + 1

	    t3 = simp(I,2) -1 + NV*(simp(I,3)-1)
 	    code_edge (t3) = code_edge (t3) + 1
	end do

	Nedges(NF) = 0
	do I = 1, NV*NV
	    if (code_edge(I).gt.2) then
		dupl_edges = .true.
		return
	    end if
	    if (code_edge(I).gt.0) NE = NE + 1
	    if (mod (code_edge (I), 2).ne.0) then
	    	Nedges(NF) = Nedges(NF) + 1
		Edges (NF, Nedges(NF), 1) = mod (I, NV) + 1
		Edges (NF, Nedges(NF), 2) = I - Edges (NF, Nedges(NF), 1) + 1
		Edges (NF, Nedges(NF), 2) = Edges (NF, Nedges(NF), 2)/ NV + 1
	    end if
	end do

	Nvertices (NF) = 0	!list of vertices on the edge
	do I = 1, 100
	    Vertices(NF, I) = 0
	end do

	do I = 1, Nedges (NF)
	    do J = 1,2
	        Vertices(NF, Edges(NF, I,J)) = 
     .		Vertices (NF, Edges(NF, I,J)) + 1
	    end do
	end do

	do I = 1, NV
	    if (Vertices (NF, I).eq.1) 
     .		write (6,*) 'Sheesh!  Vertex is at only one edge, nf is', nf
	    if (Vertices (NF,I).gt.0) then
		Nvertices (NF) = Nvertices (NF) + 1
		Vertices (NF,Nvertices(NF)) = I
	    end if
	end do

	first_edge (NF, 1) = Vertices (NF, 1)
	first_edge (NF, 2) = V
	do I = 1, Nedges (NF) 
	    if (Edges (NF, I, 1).eq.Vertices (NF, 1)) 
     .		first_edge (NF, 2) = min (first_edge (NF, 2), Edges (NF, I, 2))						!find the "smallest" edge
	end do

c	Now check to see how many edges are between vertices on the
c	boundary.for trimming the tree later.
	do I = 1, Nvertices (NF)
	    do J = I+1, Nvertices (NF)
		code_ij = Vertices(NF,I)-1 + 
     .		    NV * (Vertices(NF,J)-1)
		if (Code_edge(code_ij).ne.0) nb_edges = nb_edges + 1
	    end do
	end do

	if (NV.lt.V) then
	    Nvertices (NF) = Nvertices (NF) + 1
	    Vertices (NF,Nvertices(NF)) = NV + 1
	end if

	return
	end

	logical function Checkface (v1,v2, v3)
	implicit none
	integer v1, v2, v3
	integer Edges, simp, Nedges, NF, NV, NE, from, vert, edge,
     .		Vertices, Nvertices, V,E,F, Nb_edges, orient,
     .		first_edge, set_connect  
	logical dupl_edges, want_orientable
	common / / simp (100, 3), Edges (100, 175, 2), Nedges(100), 
     .		from(100), NF, NV, NE, vert, edge, dupl_edges,
     .		Vertices (100,100), Nvertices (100), V,E,F, Nb_edges,
     .		want_orientable, orient (100), first_edge (100,2),
     .		set_connect 

	integer s1, s2, s3, I, J, t1, t2
	
c	a function.
	logical pinch

	checkface = .false.
	if (v3.eq.v1.or.v3.eq.v2) return
c
	if (v3.lt.v1) then	!sort vertices
	    s1 = v3
	    s2 = v1
	    s3 = v2
	elseif (v3.lt.v2) then
	    s1 = v1
	    s2 = v3
	    s3 = v2
	else
	    s1 = v1
	    s2 = v2
	    s3 = v3
	end if

c	Is there already a face for these three vertices?
	do I = 1, NF
	    if (simp(I,1).eq.s1.and.simp(I,2).eq.s2.and.simp(I,3).eq.s3) 
     .	    return
	end do
c
c	check to see if we'd be making a pinchpoint
	if (pinch(v1,v2,v3)) then
c	    write (7,*) 'Pinchpoint'
c	    write (7,'(a,100(2i3,2x))') 
c     .		'Edges are', ((Edges(NF, I, J), J = 1, 2), I = 1, Nedges(NF))
c	    write (7,*) 'Pinchpoint at', v2
c	    write (7,*) 'Edge & vertex that''re being tossed out', v1,v2,v3
c	    write (7,*)
	    return
	end if
	if (pinch(v2,v1,v3)) then
c	    write (7,*) 'Pinchpoint'
c	    write (7,'(a,100(2i3,2x))') 
c     .		'Edges are', ((Edges(NF, I, J), J = 1, 2), I = 1, Nedges(NF))
c	    write (7,*) 'Pinchpoint at', v1
c	    write (7,*) 'Edge & vertex that''re being tossed out', v1,v2,v3
c	    write (7,*)
	    return
	end if
	if (pinch(v1,v3,v2)) then
c	    write (7,*) 'Pinchpoint'
c	    write (7,'(a,100(2i3,2x))') 
c     .		'Edges are', ((Edges(NF, I, J), J = 1, 2), I = 1, Nedges(NF))
c	    write (7,*) 'Pinchpoint at', v3
c	    write (7,*) 'Edge & vertex that''re being tossed out', v1,v2,v3
c	    write (7,*)
	    return
	end if
c	make a new face.
	checkface = .true.
	NF = NF + 1
	simp (NF, 1) = s1
	simp (NF, 2) = s2
	simp (NF, 3) = s3
	if (v3.eq.NV+1) then
	    NV = NV+1
	    from (NF) = vert
	else
	    from (NF) = edge
	end if

	return
	end

	logical function trim_tree ()
	implicit none
	integer Edges, simp, Nedges, NF, NV, NE, from, vert, edge,
     .		Vertices, Nvertices, V,E,F, Nb_edges, orient,
     .		first_edge, set_connect  
	logical dupl_edges, want_orientable
	common / / simp (100, 3), Edges (100, 175, 2), Nedges(100), 
     .		from(100), NF, NV, NE, vert, edge, dupl_edges,
     .		Vertices (100,100), Nvertices (100), V,E,F, Nb_edges,
     .		want_orientable, orient (100), first_edge (100,2),
     .		set_connect 

	integer S1, S2, S3, S4, DV, DE, DF, Dedges

c	function definition.
	logical orientable

	trim_tree = .true.
	if (NF.eq.F) return
	if (dupl_edges) return
	if (nedges(nf).eq.0) return

c  Now.  There are 4 basic triangling operatings:
c	S1:  add vertex:  NF -> NF+1, NE -> NE+2, NV -> NV+1, Nedges -> 
c			  Nedges+1
c	S2:  add edge:    NF -> NF+1, NE -> NE+2, Nedges -> Nedges+1
c	S3:  fill in boundary:  NF -> NF+1, NE -> NE+1, Nedges -> Nedges-1
c	S4:  fill in a triangle: NF -> NF+1, Nedges -> Nedges-3
c
c	This means that we can decide whether the current list of simplexes 
c	can lead to the desired triangulation.

	DV = V-NV
	DE = E-NE
	DF = F-NF
	if (abs(DV-DE+DF).gt.DF) return   !Euler characteristic too much off
	Dedges = - Nedges(NF)

	S1 = DV		!# of S1 operations
	do S2 = 0, max (FLOAT(DF-S1), (DE-2*S1)/2.0)+1	!try values for S2
	    S3 = DE - 2*S1 - 2*S2
	    if (S3.lt.0) goto 9
	    S4 = DF - S1 - S2 - S3
	    if (S4.lt.1) goto 9		!must have at least one of these.
	    goto 10
 9	    continue
	end do
	return			!trim the tree
 10	continue
	
c   Now, are there enough exposed vertices left over to give you the
c   necessary extra edges & faces?

	if (Vertices(NF, Nvertices(NF)).eq.NV+1) then
	    DV = DV + Nvertices(NF) - 1		!# of vertices on the boundary
	else					!plus unused ones.
	    DV = DV + Nvertices(NF)
	end if

	if ((DV*(DV-1)/2) + NE - Nb_edges.lt.E) return	!can only make
								!new edges from
								!these vertices
	if (want_orientable.and..not.orientable()) return	

	trim_tree = .false.
	return
	end

	logical function Pinch (v1,v2, v3)
	implicit none
	integer v1, v2, v3
	integer Edges, simp, Nedges, NF, NV, NE, from, vert, edge,
     .		Vertices, Nvertices, V,E,F, Nb_edges, orient,
     .		first_edge, set_connect  
	logical dupl_edges, want_orientable
	common / / simp (100, 3), Edges (100, 175, 2), Nedges(100), 
     .		from(100), NF, NV, NE, vert, edge, dupl_edges,
     .		Vertices (100,100), Nvertices (100), V,E,F, Nb_edges,
     .		want_orientable, orient (100), first_edge (100,2),
     .		set_connect 

	integer s1,s2,s3, I, match, edge_cnt
	integer edge_list(175,2)
	logical edge_tag(175)

c	This function tells if v2 would be a pinch point with this new 
c	triangulation.

	pinch = .false.
c	is (v2,v3) an edge on the boundary?
	s2 = min (v2,v3)
	s3 = max (v2,v3)

	do I = 1, Nedges (NF)
	    if (Edges(NF,I,1).eq.s2.and.Edges(NF,I,2).eq.s3) goto 10
	end do
	return		! (v2,v3) is not on the boundary
 10	continue
c
	do I = 1, Nedges (NF)
	    if (Edges (NF,I,1).eq.v2.and.
     .		(Edges(NF,I,2).ne.v1.and.Edges(NF,I,2).ne.v3)) goto 15  !keep on
 										     !looking.
	    if (Edges (NF,I,2).eq.v2.and.
     .		(Edges(NF,I,1).ne.v1.and.Edges(NF,I,1).ne.v3)) goto 15
	end do
 14	return 		!no pinch point
 15	continue		!keep on looking

c	is (v1,v2) an edge on the boundary?
	s1 = min (v1,v2)
	s2 = max (v1,v2)

	do I = 1, Nedges (NF)
	    if (Edges(NF,I,1).eq.s1.and.Edges(NF,I,2).eq.s2) goto 20
	end do
	return		! (v1,v2) is not on the boundary
 20	continue
c
c	Now are we forming a cycle of simplexes around v2?
	edge_cnt = 0
	do I = 1, NF
	    if (simp(I,1).eq.v2) then
		edge_cnt = edge_cnt + 1
	  	edge_list (edge_cnt,1) = simp(I,2)
	  	edge_list (edge_cnt,2) = simp(I,3)
	    elseif (simp(I,2).eq.v2) then
		edge_cnt = edge_cnt + 1
	  	edge_list (edge_cnt,1) = simp(I,1)
	  	edge_list (edge_cnt,2) = simp(I,3)
	    elseif (simp(I,3).eq.v2) then
		edge_cnt = edge_cnt + 1
	  	edge_list (edge_cnt,1) = simp(I,1)
	  	edge_list (edge_cnt,2) = simp(I,2)
	    end if
	end do

	do I = 1, edge_cnt
	    edge_tag(I) = .true.
	end do

c	Now see if there's a cycle around v2.  Starting with v1, can you get to
c	v3?
	match = v1
 30	continue
	do I = 1, edge_cnt
	    if (edge_tag(I)) then
	    if (edge_list(I,1).eq.match) then
	 	match = edge_list (I,2)
		edge_tag(I) = .false.
		if (match.eq.v3) then		!has to be a pinch point, we 
							!know some single edges with
							!v2 remain. 
		    pinch = .true.
		    return
		end if
		goto 30		
	    elseif (edge_list(I,2).eq.match) then
	 	match = edge_list (I,1)
		edge_tag(I) = .false.
		if (match.eq.v3) then
		    pinch = .true.
		    return
		end if
		goto 30
	    end if
	    end if
	end do		
	
	return	!no pinch point is being made
	end

	logical function orientable ()
	implicit none
	integer Edges, simp, Nedges, NF, NV, NE, from, vert, edge,
     .		Vertices, Nvertices, V,E,F, Nb_edges, orient,
     .		first_edge, set_connect  
	logical dupl_edges, want_orientable
	common / / simp (100, 3), Edges (100, 175, 2), Nedges(100), 
     .		from(100), NF, NV, NE, vert, edge, dupl_edges,
     .		Vertices (100,100), Nvertices (100), V,E,F, Nb_edges,
     .		want_orientable, orient (100), first_edge (100,2),
     .		set_connect 

	integer bdry (100*100)
	integer  t1, t2, t3, I, J

	orientable = .true.
	if (NF.eq.1) then
	    orient (NF) = 1
	    return
	end if

	orient (NF) = 0		!unset this to reset in this function.

c	first sum up the simplexes for which we already have orientation.
	do I = 1, NV*NV
	    bdry (I) = 0
	end do

	do I = 1, NF
	    if (orient(I).eq.0) goto 9			!end of simplexes for
c								!which orientation was
c								!found

	    t1 = simp(I,1) -1 + NV*(simp(I,2)-1)	!store edges in 
								!base NV code.
	    bdry ( t1) = bdry (t1) + orient(I)		!simplexes have vertices
c								!in ascending order

	    t2 = simp(I,1) -1 + NV*(simp(I,3)-1)
	    bdry (t2) = bdry (t2) - orient(I)

	    t3 = simp(I,2) -1 + NV*(simp(I,3)-1)
 	    bdry (t3) = bdry (t3) + orient(I)
	end do
 9	continue

	do J = I, NF
	    if (orient(J).eq.0) then 
c	Now add on the J'th simplex
		t1 = simp(J,1) -1 + NV*(simp(J,2)-1)
		t2 = simp(J,1) -1 + NV*(simp(J,3)-1)
		t3 = simp(J,2) -1 + NV*(simp(J,3)-1)
		if (bdry(t1).eq.1) then
	    	    orient (J) = -1
		elseif(bdry(t1).eq.-1) then
	    	    orient (J) = 1
		elseif (bdry (t2).eq.1) then
	    	    orient (J) = 1
		elseif (bdry (t2).eq.-1) then
	    	    orient (J) = -1
		elseif (bdry (t3).eq.1) then
	    	    orient (J) = -1
		elseif (bdry (t3).eq.-1) then
	    	    orient (J) = 1
		end if
	
		if (orient(J).ne.0) then
		    bdry (t1) = bdry (t1) + orient (J)
		    bdry (t2) = bdry (t2) - orient (J)
		    bdry (t3) = bdry (t3) + orient (J)
		    if (bdry (t1).eq.2.or.bdry(t1).eq.-2) then
			orientable = .false.
		 	return
		    end if
		    if (bdry (t2).eq.2.or.bdry(t2).eq.-2) then
			orientable = .false.
		 	return
		    end if
		    if (bdry (t3).eq.2.or.bdry(t3).eq.-2) then
			orientable = .false.
		 	return
		    end if
		else
		    write (6,*) 'Orientation wasn''t found ...'
		end if
	    end if
	end do
	return
	end

	logical function sym_check ()
	implicit none
	integer Edges, simp, Nedges, NF, NV, NE, from, vert, edge,
     .		Vertices, Nvertices, V,E,F, Nb_edges, orient,
     .		first_edge, set_connect  
	logical dupl_edges, want_orientable
	common / / simp (100, 3), Edges (100, 175, 2), Nedges(100), 
     .		from(100), NF, NV, NE, vert, edge, dupl_edges,
     .		Vertices (100,100), Nvertices (100), V,E,F, Nb_edges,
     .		want_orientable, orient (100), first_edge (100,2),
     .		set_connect 

	integer connect (100), I, J
	logical on_edge (100)

c    Do symmetry checking on the triangulation.

	sym_check = .true.

     	if (set_connect.ne.0) then 	!check to see each vertex isn't
					!connected to too many others.

	    do I = 1, 100
	 	connect (I) = 0
	    end do

	    do I = 1, NF	
		do J = 1,3
		    connect (Simp (I, J)) = connect (Simp (I, J)) + 1
		end do
	    end do

c	Set on_edge tag true if vertex is on the boundary
	    do I = 1, NV
		on_edge (I) = .false.
	    end do

	    do I = 1, NVertices (NF)
		on_edge (Vertices (NF, I)) = .true.
	    end do		
		
c	If vertex is on the boundary, don't want set_connect faces already
c	attached to it.
	    do I = 1, NV
		if (on_edge(I)) then
		    if (connect(I).ge.set_connect) then
		    	sym_check = .false.
		 	return
		    end if
		else 
c	if vertex is not on the boundary, want set_connect faces attached.
		    if (connect(I).ne.set_connect) then
			sym_check = .false.
			return
		    end if
		end if
	    end do
	end if

	return
	end