C--------------- Assistant Subroutines
C----- ASub.1. Draw thr arrow of the Magnetic Vector  -----------------------
	SUBROUTINE arrow(rdat,adat,nx1,nx2,ny1,ny2,tmax)
	REAL rdat(170,128),adat(170,128)
	real xxt(2),yyt(2),rx(3),ry(3)
	do i=ny1,ny2,3
        do j=nx1,nx2,3
	temp=2.0*rdat(j,i)/tmax
	if(temp.gt.0.035) then
          fa=adat(j,i)
          aff=10.0*.5
          ar1=aff*temp*cosd(fa)
          ar2=aff*temp*sind(fa)
          yyt(1)=i+0.5*ar2
          xxt(1)=j+0.5*ar1
          yyt(2)=i-0.5*ar2
          xxt(2)=j-0.5*ar1
          rx(1)=j+0.3*aff*temp*cosd(fa+20)
          ry(1)=i+0.3*aff*temp*sind(fa+20)
          rx(2)=xxt(1)
          ry(2)=yyt(1)
          rx(3)=j+0.3*aff*temp*cosd(fa-20)
          ry(3)=i+0.3*aff*temp*sind(fa-20)
          call pgsci(10)
          call pgline(2,xxt,yyt)
          call pgline(3,rx,ry)
        end if
	end do
	end do
C------ the following for the Transverse Field Cailbration
	do k=1,2
           tline=10.0*(2000./k)/tmax 
           xxt(1)=40*k+0.5*tline
	   xxt(2)=40*k-0.5*tline
 	   yyt(1)=125.
	   yyt(2)=125.
	   call pgsci(k)
	   call pgslw(4)
	   call pgline(2,xxt,yyt)
	   call pgslw(1)
	end do
	return
	END

C------- ASub.2. Draw the Observed Transverse Field with no Arrow
        SUBROUTINE obs_arr(rdat,adat,nx1,nx2,ny1,ny2,tmax)
        REAL rdat(170,128),adat(170,128)
        real xxt(2),yyt(2)
        do i=ny1,ny2,3
        do j=nx1,nx2,3
        temp=2.0*rdat(j,i)/tmax
        if(temp.gt.0.035) then
          fa=adat(j,i)
          aff=10.0*.5
          ar1=aff*temp*cosd(fa)
          ar2=aff*temp*sind(fa)
          yyt(1)=i+0.5*ar2
          xxt(1)=j+0.5*ar1
          yyt(2)=i-0.5*ar2
          xxt(2)=j-0.5*ar1
          call pgline(2,xxt,yyt)
        end if
        end do
        end do
C------ the following for the Transverse Field Cailbration
	do k=1,2
           tline=10.0*(2000./k)/tmax 
           xxt(1)=40*k+0.5*tline
	   xxt(2)=40*k-0.5*tline
 	   yyt(1)=125.
	   yyt(2)=125.
	   call pgsci(k)
	   call pgslw(4)
	   call pgline(2,xxt,yyt)
	   call pgslw(1)
	end do
        return
        END

C----- ASub.3. Arrow for a specific box 
	SUBROUTINE boxarr(bb,nbx,nby,ll)
	real bb(150),xxt(2),yyt(2),rx(3),ry(3)
	real nbx(150),nby(150)
	
 	do n=1,ll
	 fa=bb(n) 
	 aff=10.0
          ar1=aff*cosd(fa)
          ar2=aff*sind(fa)
          yyt(1)=nby(n)+0.5*ar2
          xxt(1)=nbx(n)+0.5*ar1
          yyt(2)=nby(n)-0.5*ar2
          xxt(2)=nbx(n)-0.5*ar1
	  rfa=fa+20
	  pfa=fa-20
	  if(rfa.gt.360.) rfa=rfa-360.
	  if(pfa.gt.360.) pfa=pfa-360.
          rx(1)=nbx(n)+0.3*aff*cosd(rfa)
          ry(1)=nby(n)+0.3*aff*sind(rfa)
          rx(2)=xxt(1)
          ry(2)=yyt(1)
          rx(3)=nbx(n)+0.3*aff*cosd(pfa)
          ry(3)=nby(n)+0.3*aff*sind(pfa)
          call pgsci(1)
          call pgline(2,xxt,yyt)
          call pgline(3,rx,ry)
	 end do
	call pgpoint(ll,nbx,nby,8)
	return
	END

C---ASub.4  Draw thr arrow of the Magnetic Vector for FSub.[labcn]-------------
	SUBROUTINE narrow(rdat,adat,nx1,nx2,ny1,ny2,tmax,box,nbox)
	REAL rdat(170,128),adat(170,128),box(4,nbox)
	real xxt(2),yyt(2),rx(3),ry(3)
	integer visible

        write(*,*)box,nbox
	do i=ny1,ny2,3
        do j=nx1,nx2,3
	temp=2.0*rdat(j,i)/tmax
	if(temp.gt.0.035) then
          fa=adat(j,i)
          aff=10.0*.5
          ar1=aff*temp*cosd(fa)
          ar2=aff*temp*sind(fa)
          yyt(1)=i+0.5*ar2
          xxt(1)=j+0.5*ar1
          yyt(2)=i-0.5*ar2
          xxt(2)=j-0.5*ar1
          rx(1)=j+0.3*aff*temp*cosd(fa+20)
          ry(1)=i+0.3*aff*temp*sind(fa+20)
          rx(2)=xxt(1)
          ry(2)=yyt(1)
          rx(3)=j+0.3*aff*temp*cosd(fa-20)
          ry(3)=i+0.3*aff*temp*sind(fa-20)
          call pgsci(10)
	  visible=1
          do k=1,nbox
             if((j.ge.box(1,k)).and.(j.le.box(2,k)).and.
     1            (i.ge.box(3,k)).and.(i.le.box(4,k))) then
                  visible=0
 	      end if
	  end do
          if(visible.eq.1) then
                  call pgline(2,xxt,yyt)
                  call pgline(3,rx,ry)
	  end if
        end if
	end do
	end do
C------ the following for the Transverse Field Cailbration
	do k=1,2
           tline=10.0*(2000./k)/tmax 
           xxt(1)=40*k+0.5*tline
	   xxt(2)=40*k-0.5*tline
 	   yyt(1)=125.
	   yyt(2)=125.
	   call pgsci(k)
	   call pgslw(4)
	   call pgline(2,xxt,yyt)
	   call pgslw(1)
	end do
	return
	END



C----- ASub.5. Draw thr arrow of shifting the Magnetic Vector ----------------
	SUBROUTINE mvarrow(rdat,adat,nx1,nx2,ny1,ny2,tmax,dxmv,dymv)
	REAL rdat(170,128),adat(170,128)
	real xxt(2),yyt(2),rx(3),ry(3)
	do i=ny1,ny2,3
        do j=nx1,nx2,3
	temp=2.0*rdat(j,i)/tmax
	if(temp.gt.0.035) then
          fa=adat(j,i)
          aff=10.0*.5
          ar1=aff*temp*cosd(fa)
          ar2=aff*temp*sind(fa)
	  xj=j+dxmv
	  yi=i+dymv
          yyt(1)=yi+0.5*ar2
          xxt(1)=xj+0.5*ar1
          yyt(2)=yi-0.5*ar2
          xxt(2)=xj-0.5*ar1
          rx(1)=xj+0.2*aff*temp*cosd(fa+40)
          ry(1)=yi+0.2*aff*temp*sind(fa+40)
          rx(2)=xxt(1)
          ry(2)=yyt(1)
          rx(3)=xj+0.2*aff*temp*cosd(fa-40)
          ry(3)=yi+0.2*aff*temp*sind(fa-40)
          call pgline(2,xxt,yyt)
          call pgline(3,rx,ry)
        end if
	end do
	end do
C------ the following for the Transverse Field Cailbration
	do k=1,2
           tline=10.0*(2000./k)/tmax 
           xxt(1)=40*k+0.5*tline
	   xxt(2)=40*k-0.5*tline
 	   yyt(1)=125.
	   yyt(2)=125.
	   call pgslw(4)
	   call pgline(2,xxt,yyt)
	   call pgslw(1)
	end do
	return
	END

C------ASub.6  Smooth 3 X 4 Let (512 X 512)===>(170 X 128)
	subroutine smooth(dat,flat)
        real dat(512,512),flat(170,128),max,min
	max=0.
	min=0.
        do i=1,512,4
           ii=(i+3)/4
           do j=1,510,3
              jj=(j+2)/3
              c=dat(j,i)+dat(j+1,i)+dat(j+2,i)
              c=c+dat(j,i+1)+dat(j+1,i+1)+dat(j+2,i+1)
              c=c+dat(j,i+2)+dat(j+1,i+2)+dat(j+2,i+2)
              c=c+dat(j,i+3)+dat(j+1,i+3)+dat(j+2,i+3)
           flat(jj,128-ii+1)=c/12.
	  if(max.lt.c/12.) max=c/12.
	  if(min.gt.c/12.) min=c/12.
           end do
        end do
	write(*,*)'max=  min= '
	write(*,*)max,min
        END


C----------ASub.7. Interpolating the Curve of Longitude Field Cailbration -----
        REAL FUNCTION INTERP(XX)
        REAL X(17),Y(17)
        DATA X/20.,40.,80.,160.,280.,424.,556.,677.,788.,890.,
     1   984.,1071.,1152.,1320.,1500.,1700.,2000./
        DATA Y/20.,40.,80.,160.,320.,640.,960.,1280.,1600.,1920.,
     1   2240.,2560.,2880.,3335.,3665.,3900.,4000./
 
        I=3
        N=17
111     IF(XX.LE.X(I)) THEN
           X1=XX-X(I-1)
           X2=X(I)-XX
           IF(X1.LE.X2) I=I-1
        ELSE IF(I.NE.(N-1)) THEN
           I=I+1
           GOTO 111
        END IF
        SUM=0.0
        DO K=I-1,I+1
           SS=1.0
           DO J=I-1,I+1
              IF(J.NE.K) THEN
                SS=SS*(XX-X(J))/(X(K)-X(J))
              END IF
           END DO
           SUM=SUM+SS*Y(K)
        END DO
        INTERP=SUM
        END

C------ASub.8  2D-linea Interpolation for FSub. [sub_prj]---------------------
	subroutine insert4(bb,xp,yp,pint,M,N)
	real bb(M,N),xp,yp,pint
	real p(4),pt1,pt2,xx,yy
	integer ixp,jyp	

	ixp=int(xp)
	jyp=int(yp)
	xx=xp-ixp
	yy=yp-jyp	
	p(1)=bb(ixp,jyp)
	p(2)=bb(ixp+1,jyp)
	p(3)=bb(ixp,jyp+1)
	p(4)=bb(ixp+1,jyp+1)
	pt1=p(1)+(p(2)-p(1))*xx
	pt2=p(3)+(p(4)-p(3))*xx
	pint=pt1+(pt2-pt1)*yy
	return
	END
 
C------ ASub.9. Select an Irregular Area with Cursor
	SUBROUTINE boundary(jx1,jx2,iy1,kl)
        INTEGER jx1(200),iy1(200),jx2(200),iy2(200)
        INTEGER jxt(200),iyt(200),iia(20)
        REAl xpt(50),ypt(50)
C       ***** jx1()--left side,jx2()--right side,iy1()--height,kl--points num
C       in height
 
	write(*,*)'---Select the Rigion You interest! '
	write(*,*)'--- Pgplot cursor Commond: A--select, X--exit'
c       now epecify the source !!!
c       first to find a window to correct !!!
c
          kstart=1
          npt=0
          call pgsci(5)
          call pglcur(50,npt,xpt,ypt)
          xpt(npt)=xpt(1)
          ypt(npt)=ypt(1)
c
c       make a new version to find both side of points at one time
c       find left boundary
c
        do k=1,npt
          kp1=k+1
          i1=ypt(k)
          i2=ypt(kp1)
          if(i2.gt.i1) then
           lpt=k
           goto 222 ! lpt the lowest point in this area
          endif
         enddo
c         write(*,*)'lowest point lpt=',lpt
222      do n=1,lpt-1
          np1=n+1
          i1=ypt(n)
          i2=ypt(np1)
          j1=xpt(n)
          j2=xpt(np1)
          call fside(i1,i2,j1,j2,jxt,iyt,iia(n))
c        write(*,*)'i1,i2,j1,j2',i1,i2,j1,j2
c        write(*,*)'left jxt,iytn=',n,jxt,iyt
          k1=kstart+iia(n)-2
           do k=kstart,k1
            kt=k-kstart+1
            jx1(k)=jxt(kt)
            iy1(k)=iyt(kt)
           enddo
           kstart=k1+1
         enddo
         kl=k1+1
         jx1(kl)=nint(xpt(lpt))
         iy1(kl)=nint(ypt(lpt))
         xup=xpt(1)
         yup=ypt(1)
         xlow=xpt(lpt)
         ylow=ypt(lpt)
c
c       find the right side points
c
         kstart=1
c
         do n=npt,lpt+1,-1
          np1=n-1
          i1=ypt(n)
          i2=ypt(np1)
          j1=xpt(n)
          j2=xpt(np1)
          call fside(i1,i2,j1,j2,jxt,iyt,iia(n))
c        write(*,*)'i1,i2,j1,j2',i1,i2,j1,j2
c        write(*,*)'right jxt,iytn=',n,jxt,iyt
          k1=kstart+iia(n)-2
           do k=kstart,k1
            kt=k-kstart+1
            jx2(k)=jxt(kt)
            iy2(k)=iyt(kt)
           enddo
           kstart=k1+1
	 enddo
 
         kr=k1+1
         jx2(kr)=nint(xpt(lpt))
         iy2(kr)=nint(ypt(lpt))
         nsid=kr     !point numbers in the left or right broad
        write(*,*)'nsid=',nsid
	RETURN
	END 

C------ ASub.10  Find the Boundary Points for Irregular Area
        SUBROUTINE fside(i1,i2,j1,j2,jxt,iyt,num)
        INTEGER jxt(200),iyt(200)
        num=abs(i1-i2)+1
        if(j1.eq.j2) then
          do i=i1,i2,-1
            k=abs(i-i1)+1
            jxt(k)=j1
            iyt(k)=i
          end do
        else
          qk=float(i2-i1)/float(j2-j1)
          do i=i1,i2,-1
             k=abs(i-i1)+1
             dx=abs(i-i1)/qk
             jxt(k)=nint(j1-dx)
             iyt(k)=i
          end do
        end if
        return
        END


C------ASub.11------  A SUBROUTINE FOR FAST FOURIER TRANSFORM
	subroutine d2ft(m,n,eta,a,b)
	dimension a(n,n),b(n,n)
	real a,b
	do kv=1,n
	n1=n/2
	nm1=n-1
	j=1
	do 103 i=1,nm1
	if (j.le.i) go to 100
	xt1=a(j,kv)
	xt2=b(j,kv)
	a(j,kv)=a(i,kv)
	b(j,kv)=b(i,kv)
	a(i,kv)=xt1
	b(i,kv)=xt2
100	k=n1
101     if (j.le.k) go to 102
	j=j-k
	k=k/2
	go to 101
102	j=j+k
103 	continue

	pi=3.1415926536
c	type *, pi
	do 200 l=1,m
	e=2.**l
	e1=e/2
	le=2**l
	le1=le/2
	u1=1.
	u2=0.0
	w1=cos(pi/e1)
	w2=sin(pi/e1)
	do 200 j=1,le1
	do 104 i=j,nm1,le
	ii=i+le1
	t1=a(ii,kv)*u1-eta*b(ii,kv)*u2
	t2=b(ii,kv)*u1+eta*a(ii,kv)*u2
	a(ii,kv)=a(i,kv)-t1
	b(ii,kv)=b(i,kv)-t2
	a(i,kv)=a(i,kv)+t1
	b(i,kv)=b(i,kv)+t2
104     continue
	r1=u1*w1-u2*w2
	u2=u1*w2+u2*w1
	u1=r1
200     continue
	end do
	an=float(n)
	if (eta.eq.-1.0) then
	 do kv=1,n
	  do i=1,n
	   a(i,kv)=a(i,kv)/an
	   b(i,kv)=b(i,kv)/an
	  end do
	end do
	end if	

c	vertical fourier transform
	do kh=1,n
	n1=n/2
	nm1=n-1
	j=1
	do 203 i=1,nm1
	if (j.le.i) go to 205
	xt1=a(kh,j)
	xt2=b(kh,j)
	a(kh,j)=a(kh,i)
	b(kh,j)=b(kh,i)
	a(kh,i)=xt1
	b(kh,i)=xt2
205	k=n1
201     if (j.le.k) go to 202
	j=j-k
	k=k/2
	go to 201
202	j=j+k
203 	continue
c
	do 300 l=1,m
	e=2.**l
	e1=e/2
	le=2**l
	le1=le/2
	u1=1.
	u2=0.0
	w1=cos(pi/e1)
	w2=sin(pi/e1)
	do 300 j=1,le1
	do 204 i=j,nm1,le
	ii=i+le1
	t1=a(kh,ii)*u1-eta*b(kh,ii)*u2
	t2=b(kh,ii)*u1+eta*a(kh,ii)*u2
	a(kh,ii)=a(kh,i)-t1
	b(kh,ii)=b(kh,i)-t2
	a(kh,i)=a(kh,i)+t1
	b(kh,i)=b(kh,i)+t2
204     continue
	r1=u1*w1-u2*w2
	u2=u1*w2+u2*w1
	u1=r1
300     continue
	end do
	if (eta.eq.-1.0) then
	 do j=1,n
	  do kh=1,n
	   a(kh,j)=a(kh,j)/an
	   b(kh,j)=b(kh,j)/an
	  end do
	end do
	end if	
	return     
	end

C--------ASub.12 Transfer data from Float Type to 4-digit Character Type -----	
	subroutine string(ff,str)
	real ff
	character str*5
	
	write(*,*)'value== ',ff
	nff=nint(ff)
	nch1=nff/1000
	nres=mod(nff,1000)
	str(1:1)=char(48+nch1)
	nch2=nres/100
	nres=mod(nres,100)
	str(2:2)=char(48+nch2)
	nch3=nres/10
	nres=mod(nres,10)
	str(3:3)=char(48+nch3)
	nch4=nres
	str(4:4)=char(48+nch4)
	write(*,*)'string*4== ',str
	return
	end


C-----ASub.13---  A SUBROUTINE TO CACULATE THE Bt & Abt
        subroutine btba(bx,by,bt,ba,thmax)
        DIMENSION bx(128,128),by(128,128)
        dimension bt(170,128),ba(170,128)
c
        thmax=0.0
        do j=1,128
         do i=1,128
	   li=i+20
          qq=bx(i,j)
          uu=by(i,j)
          tt=sqrt(qq**2+uu**2)
          bt(li,j)=tt
          if (tt.gt.thmax) thmax=tt
          if (qq.eq.0.00) then
           if (uu.gt.0.00) then
            ff=.5*180
           else
            ff=-.5*180
           end if
          else
           ff=atan2d(uu,qq)
          end if
        if (ff.lt.0.0) then
          ff=2*180+ff
        end if
          ba(li,j)=ff
c
         end do
        end do
        return
        END
                 

C------ASub.14.  Change the brightness or contrast for gray map
	SUBROUTINE GWHBC(A,B,IDIM,JDIM,WHITE,BLACK)
	REAL A(IDIM,JDIM),B(IDIM,JDIM),WHITE,BLACK
	character key
	real gy0,ratio,gymax,gymin
	
              WRITE(*,*)'--- Change the Contrast of Gray?<Y/N>'
              READ(*,*)KEY

              IF(KEY.EQ.'y') THEN
                 write(*,*)'input Bright=?[-128.(b)--+128.(w)] '
                 READ(*,*)gy0
                 write(*,*)'input Black-White Contrast=[0.0--1.0]'
                 READ(*,*)ratio
                 do i=1,512   
                 do j=1,512
                    B(j,i)=A(j,i)+gy0-128.0
                    gymin=min(B(j,i),gymin)
                    gymax=max(B(j,i),gymax)
                 end do
                 end do
                 WHITE=gymax*ratio
                 BLACK=gymin*ratio
	      ELSE 
	         do i=1,512
	         do j=1,512
		    B(j,i)=A(j,i)
		 enddo
	         enddo
              END IF
	    RETURN
	    END

C------- ASub.15. Display the gray image for FSub[labgy]----------
	SUBROUTINE  GGRAY(A, IDIM, JDIM, I1, I2, J1, J2,
     1             FG, BG, TR, BOX, NBOX)
	INTEGER  IDIM, JDIM, I1, I2, J1, J2, NBOX
	REAL FG, BG, TR(6)
	REAL A(IDIM,JDIM), PLOT(512,512), BOX(4,NBOX)
	INTEGER VISIBLE
 
              do M=J1, J2
              do N=I1, I2
                 visible=1
                 do k=1,nbox
                    if((N.ge.box(1,k)).and.(N.le.box(2,k)).and.
     1                 (M.ge.box(3,k)).and.(M.le.box(4,k))) then
                         visible=0
                    end if
                 end do
                 if(visible.eq.1) then
                    plot(N,M)=A(N,M)    
                 else
                    plot(N,M)=-1000.
                 end if
              end do
              end do
	   CALL PGGRAY(PLOT,IDIM,JDIM,I1,I2,J1,J2,FG,GB,TR)
	   RETURN
	   END


C      ASub.16 --- Find max or min value of array----
	SUBROUTINE FMAXMIN(ff,jdim,idim,fmax,fmin)
	real ff(jdim,idim),fmax,fmin
	integer idim,jdim

	fmax=-1.e6
	fmin=1.e6
	do i=1,idim
	do j=1,jdim
	   if(fmax.lt.ff(j,i)) fmax=ff(j,i)
	   if(fmin.gt.ff(j,i)) fmin=ff(j,i)
	end do
	end do
 	write(*,10)fmax,fmin
10      format('fmax== ',f10.3,'fmin== ',f10.3)
	return
	END

C-------ASub.17. Draw the Contour for FSub. [labcn] 
C---consist of five sub. modified by pgplot contour sources
C--they are [gcont,gplot,gconx,gcnsc,gcn01]
C-----------------------------------------------------
C*PGCONT -- contour map of a 2D data array (contour-following)
C+
	SUBROUTINE GCONT (A, IDIM, JDIM, I1, I2, J1, J2, C, NC, TR,box,nbox)
      INTEGER IDIM, JDIM, I1, J1, I2, J2, NC, nbox
      REAL A(IDIM,JDIM), C(*), TR(6)
      real box(4,nbox)
C
C Draw a contour map of an array.  The map is truncated if
C necessary at the boundaries of the viewport.  Each contour line
C is drawn with the current line attributes (color index, style, and
C width); except that if argument NC is positive (see below), the line
C style is set by PGCONT to 1 (solid) for positive contours or 2
C (dashed) for negative contours.
C
C Arguments:
C  A      (input) : data array.
C  IDIM   (input) : first dimension of A.
C  JDIM   (input) : second dimension of A.
C  I1, I2 (input) : range of first index to be contoured (inclusive).
C  J1, J2 (input) : range of second index to be contoured (inclusive).
C  C      (input) : array of NC contour levels; dimension at least NC.
C  NC     (input) : +/- number of contour levels (less than or equal
C                   to dimension of C). If NC is positive, it is the
C                   number of contour levels, and the line-style is
C                   chosen automatically as described above. If NC is
C                   negative, it is minus the number of contour
C                   levels, and the current setting of line-style is
C                   used for all the contours.
C  TR     (input) : array defining a transformation between the I,J
C                   grid of the array and the world coordinates.
C                   The world coordinates of the array point A(I,J)
C                   are given by:
C                     X = TR(1) + TR(2)*I + TR(3)*J
C                     Y = TR(4) + TR(5)*I + TR(6)*J
C                   Usually TR(3) and TR(5) are zero - unless the
C                   coordinate transformation involves a rotation or
C                   shear.
C--
C (7-Feb-1983)
C (24-Aug-1984) Revised to add the option of not automatically
C       setting the line-style. Sorry about the ugly way this is
C       done (negative NC); this is the least incompatible way of doing
C       it (TJP).
C (21-Sep-1989) Changed to call PGCONX instead of duplicating the code
C       [TJP].
C-----------------------------------------------------------------------
      INCLUDE  'pgplot.inc'
      INTEGER  I
      EXTERNAL GPLOT
C
C Save TRANS matrix.
C
      DO 10 I=1,6
          TRANS(I) = TR(I)
   10 CONTINUE
C
C Use PGCONX with external function PGCP, which applies the TRANS
C scaling.
C
	CALL GCONX (A, IDIM, JDIM, I1, I2, J1, J2, C, NC, GPLOT, box, nbox)
C
      END
C
      SUBROUTINE GPLOT (K, X, Y, Z,box,nbox)
C
C PGPLOT (internal routine): Draw one contour segment (for use by
C PGCNSC).
C
C Arguments:
C
C K (input, integer): if K=0, move the pen to (X,Y); if K=1, draw
C       a line from the current position to (X,Y); otherwise
C       do nothing.
C X (input, real): X world-coordinate of end point.
C Y (input, real): Y world-coordinate of end point.
C Z (input, real): the value of the contour level, not used by PGCP at
C       the moment.
C
C (7-Feb-1983)
C-----------------------------------------------------------------------
      INCLUDE  'pgplot.inc'
      INTEGER  K,VISBOX,nbox
      REAL     X,XX,Y,YY,Z,box(4,nbox)
C
      XX = TRANS(1) + TRANS(2)*X + TRANS(3)*Y
      YY = TRANS(4) + TRANS(5)*X + TRANS(6)*Y
      VISBOX=1
      IF (K.EQ.1) THEN
            DO I=1,NBOX
	 	IF((XX.GE.BOX(1,i)).AND.(XX.LE.BOX(2,i))
     1		.AND.(YY.GE.BOX(3,i)).AND.(YY.LE.BOX(4,i)))
     1		  THEN
	  	    VISBOX=0
		END IF
	    END DO
 	    IF(VISBOX.EQ.1) THEN
               CALL GRLINA(XX,YY)
	    ELSE 
	       CALL GRMOVA(XX,YY)
	    END IF
      ELSE IF (K.EQ.0) THEN
          CALL GRMOVA(XX,YY)
      END IF
      END
C*PGCONX -- contour map of a 2D data array (non-rectangular)
C+
	SUBROUTINE GCONX (A, IDIM, JDIM, I1, I2, J1, J2, C, NC, BPLOT,box,nbox)
      INTEGER  IDIM, JDIM, I1, J1, I2, J2, NC, nbox
      REAL     A(IDIM,JDIM), C(*), box(4,nbox)
      EXTERNAL BPLOT
c
C Draw a contour map of an array using a user-supplied plotting
C routine.  This routine should be used instead of PGCONT when the
C data are defined on a non-rectangular grid.  PGCONT permits only
C a linear transformation between the (I,J) grid of the array
C and the world coordinate system (x,y), but PGCONX permits any
C transformation to be used, the transformation being defined by a
C user-supplied subroutine. The nature of the contouring algorithm,
C however, dictates that the transformation should maintain the
C rectangular topology of the grid, although grid-points may be
C allowed to coalesce.  As an example of a deformed rectangular
C grid, consider data given on the polar grid theta=0.1n(pi/2),
C for n=0,1,...,10, and r=0.25m, for m=0,1,..,4. This grid
C contains 55 points, of which 11 are coincident at the origin.
C The input array for PGCONX should be dimensioned (11,5), and
C data values should be provided for all 55 elements.  PGCONX can
C also be used for special applications in which the height of the
C contour affects its appearance, e.g., stereoscopic views.
C
C The map is truncated if necessary at the boundaries of the viewport.
C Each contour line is drawn with the current line attributes (color
C index, style, and width); except that if argument NC is positive
C (see below), the line style is set by PGCONX to 1 (solid) for
C positive contours or 2 (dashed) for negative contours. Attributes
C for the contour lines can also be set in the user-supplied
C subroutine, if desired.
C
C Arguments:
C  A      (input) : data array.
C  IDIM   (input) : first dimension of A.
C  JDIM   (input) : second dimension of A.
C  I1, I2 (input) : range of first index to be contoured (inclusive).
C  J1, J2 (input) : range of second index to be contoured (inclusive).
C  C      (input) : array of NC contour levels; dimension at least NC.
C  NC     (input) : +/- number of contour levels (less than or equal
C                   to dimension of C). If NC is positive, it is the
C                   number of contour levels, and the line-style is
C                   chosen automatically as described above. If NC is
C                   negative, it is minus the number of contour
C                   levels, and the current setting of line-style is
C                   used for all the contours.
C  PLOT   (input) : the address (name) of a subroutine supplied by
C                   the user, which will be called by PGCONX to do
C                   the actual plotting. This must be declared
C                   EXTERNAL in the program unit calling PGCONX.
C
C The subroutine PLOT will be called with four arguments:
C      CALL PLOT(VISBLE,X,Y,Z)
C where X,Y (input) are real variables corresponding to
C I,J indices of the array A. If  VISBLE (input, integer) is 1,
C PLOT should draw a visible line from the current pen
C position to the world coordinate point corresponding to (X,Y);
C if it is 0, it should move the pen to (X,Y). Z is the value
C of the current contour level, and may be used by PLOT if desired.
C Example:
C       SUBROUTINE PLOT (VISBLE,X,Y,Z)
C       REAL X, Y, Z, XWORLD, YWORLD
C       INTEGER VISBLE
C       XWORLD = X*COS(Y) ! this is the user-defined
C       YWORLD = X*SIN(Y) ! transformation
C       IF (VISBLE.EQ.0) THEN
C           CALL PGMOVE (XWORLD, YWORLD)
C       ELSE
C           CALL PGDRAW (XWORLD, YWORLD)
C       END IF
C       END
C--
C 14-Nov-1985 - new routine [TJP].
C 12-Sep-1989 - correct documentation error [TJP].
C 22-Apr-1990 - corrected bug in panelling algorithm [TJP].
C-----------------------------------------------------------------------
      INTEGER  MAXEMX,MAXEMY
      PARAMETER (MAXEMX=50)
      PARAMETER (MAXEMY=50)
      INTEGER  I
      INTEGER  NNX,NNY, KX,KY, KI,KJ, IA,IB, JA,JB, LS, PX, PY
      LOGICAL  STYLE
C
C Check arguments.
C
c     IF (I1.LT.1 .OR. I2.GT.IDIM .OR. I1.GE.I2 .OR.
c    1    J1.LT.1 .OR. J2.GT.JDIM .OR. J1.GE.J2)
c    2          CALL PGERR('PGCONX: invalid range I1:I2, J1:J2')
      IF (NC.EQ.0) RETURN
      STYLE = NC.GT.0
      CALL PGQLS(LS)
      CALL PGBBUF
C
C Divide arrays into panels not exceeding MAXEMX by MAXEMY for
C contouring by PGCNSC.
C
CD    write (*,*) 'PGCONX window:',i1,i2,j1,j2
      NNX = I2-I1+1
      NNY = J2-J1+1
      KX = MAX(1,(NNX+MAXEMX-2)/(MAXEMX-1))
      KY = MAX(1,(NNY+MAXEMY-2)/(MAXEMY-1))
      PX = (NNX+KX-1)/KX
      PY = (NNY+KY-1)/KY
      DO 60 KI=1,KX
          IA = I1 + (KI-1)*PX
          IB = MIN(I2, IA + PX)
          DO 50 KJ=1,KY
              JA = J1 + (KJ-1)*PY
              JB = MIN(J2, JA + PY)
C
C             Draw the contours in one panel.
C
CD            write (*,*) 'PGCONX panel:',ia,ib,ja,jb
              IF (STYLE) CALL PGSLS(1)
              DO 40 I=1,ABS(NC)
                  IF (STYLE.AND.(C(I).LT.0.0)) CALL PGSLS(2)
 		    CALL GCNSC(A,IDIM,JDIM,IA,IB,JA,JB,C(I),BPLOT, box ,nbox)
                  IF (STYLE) CALL PGSLS(1)
   40         CONTINUE
   50     CONTINUE
   60 CONTINUE
C
      CALL PGSLS(LS)
      CALL PGEBUF
      END


	SUBROUTINE GCNSC (Z, MX, MY, IA, IB, JA, JB, Z0, BPLOT, box, nbox)
      INTEGER MX, MY, IA, IB, JA, JB, nbox
      REAL Z(MX,*), box(4,nbox)
      REAL Z0
      EXTERNAL BPLOT
C
C PGPLOT (internal routine): Draw a single contour.  This routine is
C called by PGCONT, but may be called directly by the user.
C
C Arguments:
C
C Z (real array dimension MX,MY, input): the array of function values.
C MX,MY (integer, input): actual declared dimension of Z(*,*).
C IA,IB (integer, input): inclusive range of the first index of Z to be
C       contoured.
C JA,JB (integer, input): inclusive range of the second index of Z to
C       be contoured.
C Z0 (real, input): the contour level sought.
C PLOT (the name of a subroutine declared EXTERNAL in the calling
C       routine): this routine is called by PGCNSC to do all graphical
C       output. The calling sequence is CALL PLOT(K,X,Y,Z) where Z is
C       the contour level, (X,Y) are the coordinates of a point (in the
C       inclusive range I1<X<I2, J1<Y<J2, and if K is 0, the routine is
C       to move then pen to (X,Y); if K is 1, it is to draw a line from
C       the current position to (X,Y).
C
C NOTE:  the intervals (IA,IB) and (JA,JB) must not exceed the
C dimensions of an internal array. These are currently
C set at 50.
C--
C 17-Sep-1989 - Completely rewritten [TJP]. The algorithm is my own,
C               but it is probably not original. It could probably be
C               coded more briefly, if not as clearly.
C-----------------------------------------------------------------------
      INTEGER UP, DOWN, LEFT, RIGHT
      PARAMETER (UP=1, DOWN=2, LEFT=3, RIGHT=4)
      INTEGER  MAXEMX, MAXEMY
      PARAMETER (MAXEMX=50, MAXEMY=50)
C
      LOGICAL FLAGS(MAXEMX,MAXEMY,2), RANGE
      INTEGER I, J, II, JJ
      REAL Z1, Z2, Z3, P, P1, P2
C
C The statement function RANGE decides whether a contour at level P
C crosses the line between two gridpoints with values P1 and P2. It is
C important that a contour cannot cross a line with equal endpoints.
C
      RANGE (P,P1,P2) = (P.GT.MIN(P1,P2)) .AND. (P.LE.MAX(P1,P2))
     1                  .AND. (P1.NE.P2)
C
C Check for errors.
C
      IF ( (IB-IA+1) .GT. MAXEMX .OR.  (JB-JA+1) .GT. MAXEMY ) THEN
          CALL GRWARN('PGCNSC - array index range exceeds'//
     1                ' built-in limit of 50')
          RETURN
      END IF
C
C Initialize the flags. The first flag for a gridpoint is set if
C the contour crosses the line segment to the right of the gridpoint
C (joining [I,J] to [I+1,J]); the second flag is set if if it crosses
C the line segment above the gridpoint (joining [I,J] to [I,J+1]).
C The top and right edges require special treatment. (For purposes
C of description only, we assume I increases horizontally to the right
C and J increases vertically upwards.)
C
      DO 20 I=IA,IB
          II = I-IA+1
          DO 10 J=JA,JB
              JJ = J-JA+1
              Z1 = Z(I,J)
              FLAGS(II,JJ,1) = .FALSE.
              FLAGS(II,JJ,2) = .FALSE.
              IF (I.LT.IB) THEN
                Z2 = Z(I+1,J)
                IF (RANGE(Z0,Z1,Z2)) FLAGS(II,JJ,1) = .TRUE.
              END IF
              IF (J.LT.JB) THEN
                Z3 = Z(I,J+1)
                IF (RANGE(Z0,Z1,Z3)) FLAGS(II,JJ,2) = .TRUE.
              END IF
   10     CONTINUE
   20 CONTINUE
C
C Search the edges of the array for the start of an unclosed contour.
C Note that (if the algorithm is implemented correctly) all unclosed
C contours must begin and end at the edge of the array. When one is
C found, call PGCN01 to draw the contour, telling it the correct
C starting direction so that it follows the contour into the array
C instead of out of it.
C
C Bottom edge.
C
      J = JA
      JJ = J-JA+1
      DO 26 I=IA,IB-1
          II = I-IA+1
          IF (FLAGS(II,JJ,1)) CALL GCN01(Z, MX, MY, IA, IB, JA, JB,
     1  		 Z0, BPLOT, FLAGS, I, J, UP, box, nbox)
   26 CONTINUE
C
C Right edge.
C
      I = IB
      II = I-IA+1
      DO 27 J=JA,JB-1
          JJ = J-JA+1
          IF (FLAGS(II,JJ,2)) CALL GCN01(Z, MX, MY, IA, IB, JA, JB,
     1    		Z0, BPLOT, FLAGS, I, J, LEFT, box, nbox)
   27 CONTINUE
C
C Top edge.
C
      J = JB
      JJ = J-JA+1
      DO 28 I=IB-1,IA,-1
          II = I-IA+1
          IF (FLAGS(II,JJ,1)) CALL GCN01(Z, MX, MY, IA, IB, JA, JB,
     1      		 Z0, BPLOT, FLAGS, I, J, DOWN, box, nbox)
   28 CONTINUE
C
C Left edge.
C
      I = IA
      II = I-IA+1
      DO 29 J=JB-1,JA,-1
          JJ = J-JA+1
          IF (FLAGS(II,JJ,2)) CALL GCN01(Z, MX, MY, IA, IB, JA, JB,
     1       		Z0, BPLOT, FLAGS, I, J, RIGHT, box, nbox)
   29 CONTINUE
C
C Now search the interior of the array for a crossing point, which will
C lie on a closed contour (because all unclosed contours have been
C eliminated). It is sufficient to search just the horizontal crossings
C (or the vertical ones); any closed contour must cross a horizontal
C and a vertical gridline. PGCN01 assumes that when it cannot proceed
C any further, it has reached the end of a closed contour. Thus all
C unclosed contours must be eliminated first. In this case the starting
C direction supplied to PGCN01 is not important. (It could be set, say,
C to circle maxima clockwise and minima anticlockwise, but that
C refinement is not needed here.)
C
      DO 40 I=IA+1,IB-1
          II = I-IA+1
          DO 30 J=JA+1,JB-1
              JJ = J-JA+1
              IF (FLAGS(II,JJ,1)) CALL GCN01(Z, MX, MY, IA, IB, JA, JB,
     1    			 Z0, BPLOT, FLAGS, I, J, UP, box, nbox)
   30     CONTINUE
   40 CONTINUE
C
C We didn't find any more crossing points: we're finished.
C
      RETURN
      END


	SUBROUTINE GCN01(Z, MX, MY, IA, IB, JA, JB, Z0, BPLOT,
     1 			FLAGS, IS, JS, SDIR, box ,nbox)
C
C Support routine for PGCNSC. This routine draws a continuous contour,
C starting at the specified point, until it either crosses the edge of
C the array or closes on itself.
C-----------------------------------------------------------------------
      INTEGER UP, DOWN, LEFT, RIGHT, nbox
      PARAMETER (UP=1, DOWN=2, LEFT=3, RIGHT=4)
      INTEGER  MAXEMX, MAXEMY
      PARAMETER (MAXEMX=50, MAXEMY=50)
      LOGICAL FLAGS(MAXEMX,MAXEMY,2)
      INTEGER MX, MY, IA, IB, JA, JB, IS, JS, I, J, II, JJ, DIR, SDIR
      REAL Z(MX,*), box(4,nbox)
      REAL Z0, X, Y, STARTX, STARTY
      EXTERNAL BPLOT
C
      I = IS
      J = JS
      DIR = SDIR
      II = 1+I-IA
      JJ = 1+J-JA
      IF (DIR.EQ.UP .OR. DIR.EQ.DOWN) THEN
          X = REAL(I) + (Z0-Z(I,J))/(Z(I+1,J)-Z(I,J))
          Y = REAL(J)
      ELSE
          X = REAL(I)
          Y = REAL(J) + (Z0-Z(I,J))/(Z(I,J+1)-Z(I,J))
      END IF
CD    WRITE (*,*) 'SEGMENT'
C
C Move to start of contour and record starting point.
C
      CALL BPLOT(0, X, Y, Z0, box, nbox)
      STARTX = X
      STARTY = Y
C
C We have reached grid-point (I,J) going in direction DIR (UP, DOWN,
C LEFT, or RIGHT). Look at the other three sides of the cell we are
C entering to decide where to go next. It is important to look to the
C two sides before looking straight ahead, in order to avoid self-
C intersecting contours. If all 3 sides have unused crossing-points,
C the cell is "degenerate" and we have to decide which of two possible 
C pairs of contour segments to draw; at present we make an arbitrary 
C choice. If we have reached the edge of the array, we have
C finished drawing an unclosed contour. If none of the other three
C sides of the cell have an unused crossing-point, we must have
C completed a closed contour, which requires a final segment back to
C the starting point.
C
  100 CONTINUE
CD    WRITE (*,*) I,J,DIR
      II = 1 + I - IA
      JJ = 1 + J - JA
      GOTO (110, 120, 130, 140), DIR
C
C DIR = UP
C
  110 CONTINUE
      FLAGS(II,JJ,1) = .FALSE.
      IF (J.EQ.JB) THEN
          RETURN
      ELSE IF (FLAGS(II,JJ,2)) THEN
          DIR = LEFT
          GOTO 200
      ELSE IF (FLAGS(II+1,JJ,2)) THEN
          DIR = RIGHT
          I = I+1
          GOTO 200
      ELSE IF (FLAGS(II,JJ+1,1)) THEN
C!        DIR = UP
          J = J+1
          GOTO 250
      ELSE
          GOTO 300
      END IF
C
C DIR = DOWN
C
  120 CONTINUE
      FLAGS(II,JJ,1) = .FALSE.
      IF (J.EQ.JA) THEN
          RETURN
      ELSE IF (FLAGS(II+1,JJ-1,2)) THEN
          DIR = RIGHT
          I = I+1
          J = J-1
          GOTO 200
      ELSE IF (FLAGS(II,JJ-1,2)) THEN
          DIR = LEFT
          J = J-1
          GOTO 200
      ELSE IF (FLAGS(II,JJ-1,1)) THEN
C!        DIR = DOWN
          J = J-1
          GOTO 250
      ELSE
          GOTO 300
      END IF
C
C DIR = LEFT
C
  130 CONTINUE
      FLAGS(II,JJ,2) = .FALSE.
      IF (I.EQ.IA) THEN
          RETURN
      ELSE IF (FLAGS(II-1,JJ,1)) THEN
          DIR = DOWN
          I = I-1
          GOTO 250
      ELSE IF (FLAGS(II-1,JJ+1,1)) THEN
          DIR = UP
          I = I-1
          J = J+1
          GOTO 250
      ELSE IF (FLAGS(II-1,JJ,2)) THEN
C!        DIR = LEFT
          I = I-1
          GOTO 200
      ELSE
          GOTO 300
      END IF
C
C DIR = RIGHT
C
  140 CONTINUE
      FLAGS(II,JJ,2) = .FALSE.
      IF (I.EQ.IB) THEN
          RETURN
      ELSE IF (FLAGS(II,JJ+1,1)) THEN
          DIR = UP
          J = J+1
          GOTO 250
      ELSE IF (FLAGS(II,JJ,1)) THEN
          DIR = DOWN
          GOTO 250
      ELSE IF (FLAGS(II+1,JJ,2)) THEN
C!        DIR = RIGHT
          I = I+1
          GOTO 200
      ELSE
          GOTO 300
      END IF
C
C Draw a segment of the contour.
C

  200 X = REAL(I)
      Y = REAL(J) + (Z0-Z(I,J))/(Z(I,J+1)-Z(I,J))
      CALL BPLOT(1,X,Y,Z0, box, nbox)
      GOTO 100
  250 X = REAL(I) + (Z0-Z(I,J))/(Z(I+1,J)-Z(I,J))
      Y = REAL(J)
      CALL BPLOT(1,X,Y,Z0, box, nbox)
      GOTO 100
C
C Close the contour and go look for another one.
C
  300 CALL BPLOT(1,STARTX,STARTY,Z0, box ,nbox)
      RETURN
      END


C_________________ASub.18. Set line style
C*PGSLS -- set line style
C+
      SUBROUTINE GSLS (LS)
      INTEGER  LS
C
C Set the line style attribute for subsequent plotting. This
C attribute affects line primitives only; it does not affect graph
C markers, text, or area fill.
C Five different line styles are available, with the following codes:
C 1 (full line), 2 (dashed), 3 (dot-dash-dot-dash), 4 (dotted),
C 5 (dash-dot-dot-dot). The default is 1 (normal full line).
C
C Argument:
C  LS     (input)  : the line-style code for subsequent plotting
C                    (in range 1-5).
C--
C  8-Aug-1985 - new routine, equivalent to GRSETLS [TJP].
C  3-Jun-1984 - add GMFILE device [TJP].
C-----------------------------------------------------------------------
      CALL GSETS(LS)
      END



C*GRSETLS -- set line style
C+
      SUBROUTINE GSETS (IS)
      INTEGER IS
C
C GRPCKG: Set the line style for subsequent plotting on the current
C device. The different line styles are generated in hardware on
C some devices and by GRPCKG software for the other devices. Five
C different line styles are available, with the following codes:
C 1 (full line), 2 (dashed), 3 (dot-dash-dot-dash), 4 (dotted),
C 5 (dash-dot-dot-dot). The default is 1 (normal full line). Line
C style is ignored when drawing characters, which are always drawn with
C a full line.
C
C Argument:
C
C IS (input, integer): the line-style code for subsequent plotting on
C       the current device (in range 1-5).
C--
C  9-Feb-1983 - [TJP].
C  3-Jun-1984 - add GMFILE device [TJP].
C  5-Aug-1986 - add GREXEC support [AFT].
C 21-Feb-1987 - If needed, calls begin picture [AFT].
C 19-Jan-1987 - fix bug in GREXEC call [TJP].
C 16-May-1989 - fix bug for hardware line dash [TJP].
C-----------------------------------------------------------------------
      INCLUDE 'pgplot.inc'
      INTEGER I, L, IDASH, NBUF,LCHR
      REAL    RBUF(6),TMP
      CHARACTER*10 CHR
      REAL PATERN(8,5)
C
      DATA PATERN/ 8*10.0,
     1             8*2.5,
     2             8.0, 6.0, 1.0, 6.0, 8.0, 6.0, 1.0, 6.0,
     3             1.0, 6.0, 1.0, 6.0, 1.0, 6.0, 1.0, 6.0,
     4             8.0, 6.0, 1.0, 6.0, 1.0, 6.0, 1.0, 6.0 /
C
      IF (GRCIDE.LT.1) THEN
          CALL GRWARN('GRSETLS - no graphics device is active.')
          RETURN
      END IF
C
      I = IS
      IF (I.LT.1 .OR. I.GT.5) THEN
          CALL GRWARN('GRSETLS - invalid line-style requested.')
          I = 1
      END IF
C
C Inquire if hardware dash is available.
C
      IDASH=0
      CALL GREXEC(GRGTYP, 4,RBUF,NBUF,CHR,LCHR)
      IF(CHR(3:3).EQ.'D') IDASH=1
C
C Set up for hardware dash.
C
      IF(IDASH.NE.0) THEN
          GRDASH(GRCIDE) = .FALSE.
          IF (GRPLTD(GRCIDE)) THEN
              RBUF(1)=I
              NBUF=1
              CALL GREXEC(GRGTYP,19,RBUF,NBUF,CHR,LCHR)
          END IF
C
C Set up for software dash.
C
      ELSE
          IF (I.EQ.1) THEN
              GRDASH(GRCIDE) = .FALSE.
          ELSE
              GRDASH(GRCIDE) = .TRUE.
              GRIPAT(GRCIDE) = 1
              GRPOFF(GRCIDE) = 0.0
              CALL GREXEC(GRGTYP, 6,RBUF,NBUF,CHR,LCHR)
              TMP=RBUF(4)/1000.
              DO 10 L=1,8
                  GRPATN(GRCIDE,L) = PATERN(L,I)*TMP
   10         CONTINUE
          END IF
      END IF
      GRSTYL(GRCIDE) = I
      END
