c
c Plotting of arrows from a vector field data file
c
c Fortran version written by H. Nakanishi, need to be compiled with "-lpepl".
c
	program flow
c
c
c
	character*12 if
	dimension x(10000),y(10000),z1(10000),z2(10000),r(10000)
c
c       Use subroutines to do the work
c
	call init(if,rmin,afac,aoff,xmin,xmax,ymin,ymax,x,y,z1,z2,r,n)
	call disp(if,rmin,afac,aoff,xmin,xmax,ymin,ymax,x,y,z1,z2,r,n)
	stop
	end
c
	subroutine init(if,rmin,afac,aoff,xmin,xmax,ymin,ymax
     c		,x,y,z1,z2,r,n)
c
c       Initialize variables
c
	character*12 if
        dimension x(1),y(1),z1(1),z2(1),r(1)
	print *,'Flow: input file name ?'
	read(5,15) if
15	format(a12)
	open(1,file=if)
	rewind(1)
	rmin=0.
	rmax=0.
	n=0
	do 10 i=1,10000
	read(1,*,end=20) x(i),y(i),z1(i),z2(i)
	if(i.eq.1) then
		xmin=x(i)
		xmax=x(i)
		ymin=y(i)
		ymax=y(i)
	endif
	n=n+1
	r(i)=sqrt(z1(i)**2+z2(i)**2)
	if(r(i).lt.rmin) rmin=r(i)
	if(r(i).gt.rmax) rmax=r(i)
	if(x(i).lt.xmin) xmin=x(i)
	if(x(i).gt.xmax) xmax=x(i)
	if(y(i).lt.ymin) ymin=y(i)
	if(y(i).gt.ymax) ymax=y(i)
10	continue
20	close(1)
	print *,'max vector length = ',rmax,', min = ',rmin
	if(rmax.ne.rmin) then
		print *,'map the max and min lengths to: '
		read(5,*) amax,amin
		afac=(amax-amin)/(rmax-rmin)
		aoff=amin
	else
		print *,'map the arrow length to: '
		read(5,*) amin
		afac=0.
		aoff=amin
	endif
	return
	end
c
	subroutine disp(if,rmin,afac,aoff,xmin,xmax,ymin,ymax
     c		,x,y,z1,z2,r,n)
c
c
c
	character*12 if
	character*18 if2
	dimension x(1),y(1),z1(1),z2(1),r(1)
	call usrmon(.true.,.false.,-0.5,9.,-0.5,12.)
	call pltlun(19,.true.,.false.)
	call pltlfn('graph.out')
	call plots
	call plot(0.,3.5,-3)
	call text(0.5,6.3,0.2,'Vector Field Diagram',0.,20,0)
	if2="Data: "//if
	call text(1.5,5.8,0.2,if2,0.,12,0)
	call plot(0.,0.,3)
	call plot(5.,0.,2)
	call plot(5.,5.,2)
	call plot(0.,5.,2)
	call plot(0.,0.,2)
	if(xmin.ne.xmax) then
		xfac=5./(xmax-xmin)
		xoff=0.
	else
		xfac=1.
		xoff=2.5
	endif
	if(ymin.ne.ymax) then
		yfac=5./(ymax-ymin)
		yoff=0.
	else
		yfac=0.
		yoff=2.5
	endif
	do 10 i=1,n
		x0=xfac*(x(i)-xmin)+xoff
		y0=yfac*(y(i)-ymin)+yoff
		r0=sqrt(z1(i)**2+z2(i)**2)
		a0=afac*(r0-rmin)+aoff
		if(z1(i).ne.0) then
			th=atan(z2(i)/z1(i))*180./3.14159
		else
			th=90.
			if(z2(i).lt.0) th=270.
		endif
		if(z1(i).lt.0) th=th+180.
		if(th.lt.0) th=th+360.
		call marker(x0,y0,a0,20,th,-1)
10	continue
	call number(0.,-0.5,0.1,xmin,0.,'1p,e12.5')
	call number(4.,-0.5,0.1,xmax,0.,'1p,e12.5')
	call number(-0.3,0.,0.1,ymin,90.,'1p,e12.5')
	call number(-0.3,4.,0.1,ymax,90.,'1p,e12.5')
	call plot(0.,0.,999)
	return
	end
