program cl4_check

c 11/24/99
c to remake
c rlogin parep2.lbl.gov -l merrill
c cd $MDOCS/cl4_check
c module load lang
c f77 -Bstatic -g -o cl4_check.exe cl4_check.f
c dbx cl4_check.exe
c func causecalc
c display ii
c stop at 90
c cd $MYMAPS/tr940115/erd_sf90
c run <orig_pop3069.cl4

c parep2.lbl.gov:$MDOCS/cl4_check/cl4_check.f 11/24/99
c from sy$seedis:[seedis.map_routines]cl4_check.for 9/6/96
c check validity of a file in cl4 format
c for each polygon:
c checks for negative area (OK in first polygon)
c checks for causeways (OK)
c checks for crossing segments (not OK)
c calculates area of each polygon

c reads from unit 5 (standard input)
c writes to unit 6 (standard output)

c	warning - maximum of 4000 polygons, 900 points per polygon
c	warning - not implemented for more than two causeways in a polygon
c
c	number of points in current polygon
	integer npts

c	max number of polygons
	data maxpoly /4000/
c	max number of points in any one polygon
      	data maxnpts /900/
c	point arrays for one polygon
        integer jx(900),jy(900)
c	point array for all polygons
	integer jxarray(4000,900), jyarray(4000,900)
c	twice area of each polygon
	integer iarea
c	net area of map (should be zero if boundary polygon present)
	integer iareatot
c	number of self-intersections for each segment of a polygon
	integer mcross(900)
c	number of causeways for each segment of a polygon
	integer mcause0(900)
c	number of self-intersections in a polygon
	integer ncross
c	number of causeways in a polygon
	integer ncause
c	number of self-intersections in whole map
	integer ncrosstot
c	number of causeways in whole map
	integer ncausetot

c	5 geocodes for each polygon - maximum of 4000 polygons
	integer jarray(5,4000)
c	number of points in each polygon
	integer nptsarray(4000)

c	sequence numbers of points in orginal polygon
	integer iseq0(900)
c	no of points in each subpolygon of a polygon
	integer mpts(11)
c	causeway ids in each subpolygon
	integer mcause(900,11)
c	sequence numbers of points in each subpolygon
	integer iseq(900,11)
c
c Input.
c  for005:  CL4 format from RLInt (input)
c  for006:  output file
c
c
c INCLUDED files:      FORMAT.INC

c	input CL4 file
cVMS      open(unit=5,readonly,type='OLD')
      open(unit=5,type='OLD')

c	output file
cVMS      open(unit=6,carriagecontrol='LIST',status='NEW')
      open(unit=6,status='NEW')

c	implemented only for 4-geocode levels eg clsc,nmcdtr80,etc
	level=41

c	logical unit of CL4 input file, already open
	iunitcl4=5
c	logical unit of output file
	iotmsg=6

	call read_cl4(iunitcl4,level,maxpoly,maxnpts,x0,y0,scalar,
     x	ixmin,ixmax,iymin,iymax,fconv,maxd,ngeoplus2,jarray,nptsarray,
     x	jxarray,jyarray,nbpt)

c nickel header as stored in CL4 file
c	write(iotmsg,710)x0,y0,scalar,ixmin,ixmax,iymin,iymax

c	total new area of map
	iareatot = 0
c	total causeways in whole map
	ncausetot = 0
c	total intersections in whole map
	ncrosstot = 0

c	max and min x,y of map
	xmin=99999
	ymin=99999
	xmax=-99999
	ymax=-99999

	do 1402 i=1,maxd
c	begin loop over polygons
	j1 = jarray(1,i)
	j2 = jarray(2,i)
	j3 = jarray(3,i)
	j4 = jarray(4,i)
	j5 = jarray(5,i)

c	number of points in polygon i
	npts = nptsarray(i)
c	write (6,*) ' cl4_check/i,npts=',i,npts
c	number of causeways in this polygon
	ncause = 0
c	number of intersections in this polygon
	ncross = 0

	do 1403 ii=1,npts
c	number of causeways for segment ii
	mcause0(ii)=0
c	number of crossings in segment ii
	mcross(ii)=0
1403	continue

c	CL4 format: positive areas counterclockwise
c	NICKEL format: positive areas clockwise
c	jx,jy array: positive areas clockwise

c	calculate 2*area of this polygon
	call areacalc(i,npts,jxarray,jyarray,iarea)

c	store ncause = no of causeways in this polygon

c	THE FOLLOWING DOES NOT WORK IN UNIX!  IS THIS A SYSTEM BUG?
c	call causecalc(i,npts,jxarray,jyarray,ncause,mcause0)
	call causecalc(i,nptsarray(i),jxarray,jyarray,ncause,mcause0)

	do 1400 ii=1,nptsarray(i)
c	pass over points in same order
		iseq0(ii)=ii
1400	continue

c	Now we have
c		i=polygon sequence number
c		npts=pts in polygon
c		ncause=no of causeways in polygon
c		for ii=1,npts:
c			jxarray(i,ii),jyarray(i,ii) = coords
c			mcause0(ii)=causeway id of pt ii
c			iseq0(ii)=sequence number (in jxarray) of ii

c	simple case: zero causeways in polygon i
	if (ncause.eq.0) then
c		store ncross = number of crossing segments in polygon
c		for ii=1,npts:
c			mcross(ii) = intersections with segment ii
c		subpolygon sequence no within polygon
 		call crosscalc(i,nptsarray(i),iseq0,
     x		jxarray,jyarray,ncross,mcross)
	else if (ncause.gt.0) then
		call splitpoly(nptsarray(i),ncause,mcause0,iseq0,
     x		mpts,mcause,iseq,iotmsg)

c		write(iotmsg,1300) ncause
c		write(iotmsg,1300) 
c		write(iotmsg,1300) npts
c		write(iotmsg,1300) (iseq0(jj),jj=1,npts)
c		write(iotmsg,1300) (mcause0(jj),jj=1,npts)
c		write(iotmsg,1300)
c		write(iotmsg,1300) mpts(1)
c		write(iotmsg,1300) (iseq(jj,1),jj=1,mpts(1))
c		write(iotmsg,1300) (mcause(jj,1),jj=1,mpts(1))
c		write(iotmsg,1300)
c		write(iotmsg,1300) mpts(2)
c		write(iotmsg,1300) (iseq(jj,2),jj=1,mpts(2))
c		write(iotmsg,1300) (mcause(jj,2),jj=1,mpts(2))
c		write(iotmsg,1300)
		if (ncause.eq.2) then
c		write(iotmsg,1300) mpts(3)
c		write(iotmsg,1300) (iseq(jj,3),jj=1,mpts(3))
c		write(iotmsg,1300) (mcause(jj,3),jj=1,mpts(3))
c		write(iotmsg,1300)
c		write(iotmsg,1300) mpts(4)
c		write(iotmsg,1300) (iseq(jj,4),jj=1,mpts(4))
c		write(iotmsg,1300) (mcause(jj,4),jj=1,mpts(4))
c		write(iotmsg,1300)
		endif
1300		format (1x,'c',25i3)

c		check for intersections within each of the subpolygons,
c		then for intersections of the two subpolygons
		if (ncause.eq.1) then
c			one causeways, check subpoly 1 and 2
 			call crosscalc(i,mpts(1),iseq(1,1),
     x			jxarray,jyarray,ncross,mcross)
 			call crosscalc(i,mpts(2),iseq(1,2),
     x			jxarray,jyarray,ncross,mcross)
 			call crosscalc2(i,1,2,mpts(1),iseq(1,1),
     x			jxarray,jyarray,ncross,mcross)
		endif
		if (ncause.eq.2) then
c			two causeways, check subpoly 1,3,4
c			then intersections of subpolygon pairs
 			call crosscalc(i,mpts(1),iseq(1,1),
     x			jxarray,jyarray,ncross,mcross)
 			call crosscalc(i,mpts(3),iseq(1,3),
     x			jxarray,jyarray,ncross,mcross)
 			call crosscalc(i,mpts(4),iseq(1,4),
     x			jxarray,jyarray,ncross,mcross)
 			call crosscalc2(i,1,3,mpts(1),iseq(1,1),
     x			jxarray,jyarray,ncross,mcross)
 			call crosscalc2(i,1,4,mpts(1),iseq(1,1),
     x			jxarray,jyarray,ncross,mcross)
 			call crosscalc2(i,3,4,mpts(1),iseq(1,1),
     x			jxarray,jyarray,ncross,mcross)
		endif
	else
c		(ncause.gt.2) not implemented
	endif

c	maximum and minimum x,y of map
	do 1404 ii=1,npts
	jx(ii)=jxarray(i,ii)
	jy(ii)=jyarray(i,ii)
	ixmin=min(ixmin,jx(ii))
	ixmax=max(ixmax,jx(ii))
	iymin=min(iymin,jy(ii))
	iymax=max(iymax,jy(ii))
1404	continue

c	total area of map
	iareatot = iareatot + iarea
c	total number of causeways in map
	ncausetot = ncausetot+ncause
c	total number of illegal crossings in map
	ncrosstot = ncrosstot+ncross

c	write out negative-area polygons,
c	polygons with causeways, and polygons
c	with invalid crossing segments
	if (iarea.le.0.or.ncause.gt.0.or.ncross.gt.0) then
		write (iotmsg,206) 
     x		i,j1,j2,j3,j4,j5,npts,iarea,ncause,ncross
	endif
206	format(1x,7i6,i12,2i4)

	if (ncross.gt.0.or.ncause.gt.0) then
c		loop over points in polygon i
		do 207 ii=npts,1,-1
c		ii minus one
		iim1 = ii-1
		if (ii.eq.1) then
			iim1=npts
		endif
		if (mcause0(ii).eq.0.and.mcross(ii).eq.0) then
			write(iotmsg,2060) ii,jx(ii),jy(ii),
     x			jx(iim1),jy(iim1),mcause0(ii),mcross(ii)
		else if (mcause0(ii).eq.0.and.mcross(ii).gt.0) then
			write(iotmsg,2061) ii,jx(ii),jy(ii),
     x			jx(iim1),jy(iim1),mcause0(ii),mcross(ii)
		else if (mcause0(ii).gt.0.and.mcross(ii).eq.0) then
			write(iotmsg,2062) ii,jx(ii),jy(ii),
     x			jx(iim1),jy(iim1),mcause0(ii),mcross(ii)
		else
			write(iotmsg,2063) ii,jx(ii),jy(ii),
     x			jx(iim1),jy(iim1),mcause0(ii),mcross(ii)
		endif
2060		format (i4,1x,4i6,'      ',i2,'      ',i2)
2061		format (i4,1x,4i6,'      ',i2,' cross',i2)
2062		format (i4,1x,4i6,' cause',i2,'      ',i2)
2063		format (i4,1x,4i6,' cause',i2,' cross',i2)
207		continue
	endif

c	end of polygon loop i
1402	continue

c nickel header with ixmin etc. recalculated
	write(iotmsg,710)x0,y0,scalar,ixmin,ixmax,iymin,iymax

	write(iotmsg,701) iareatot,ncausetot,ncrosstot
701	format(i12,' = twice net map area',
     x	/,i12,' = number of causeways',
     x	/,i12,' = number of intersections')


cVMS	call exit
	call exit(0)

	INCLUDE 'format.inc'

	end

	include 'fread_cl4.sub'
	include 'areacalc.sub'
	include 'causecalc.sub'
	include 'splitpoly.sub'
	include 'splitpoly0.sub'
	include 'crosscalc.sub'
	include 'crosscalc2.sub'
	include 'xsect.sub'