# Rcjp's Weblog

## October 27, 1995

### Teaching MSc Numerical Methods

Filed under: physics — rcjp @ 11:51 am

Around this time I taught a few groups of MSc students some basic numerical methods – mainly to support another class in soliton theory.

One of the first programs was to plot optical reflectance against angle at a metal glass interface – this is interesting because at certain angles you will get a surface wave traveling in the metal film.

Mostly though this was just to get students comfortable using 2D plotting NAG routines…

``````
PROGRAM REFLECT
COMPLEX*16 ETA1, ETA2, ETA3, COSTH2, COSTH3
COMPLEX*16 RTEMP1, RTEMP2, R12, R23, R, ALPHA2, II
REAL*8 THETA1, COSTH1, SCALE, CONVERT
PARAMETER (NPTS=100, CONVERT=3.141592D0/180.0D0)
REAL*8 D, W, C, THETA(NPTS), RR(NPTS)

II   = DCMPLX(0.0D0, 1.0D0)
ETA1 = DCMPLX(1.5D0**2, 0.0D0)
ETA2 = (0.044D0, 2.4D0)**2
ETA3 = (1.0D0, 0.0D0)
D    = 40.0D-9
W    = 3.2D15
C    = 3.0D8

SCALE = 90.0D0/(NPTS-1)
DO 10 I=1, NPTS
THETA(I) = (I-1) * SCALE
THETA1 = THETA(I)*CONVERT
COSTH1 = DCOS(THETA1)
COSTH2 = SQRT( 1.0D0 - (ETA1 / ETA2) * DSIN(THETA1)**2 )
ths = acos(costh2)
COSTH3 = SQRT( 1.0D0 - (ETA1 / ETA3) * DSIN(THETA1)**2 )
RTEMP1 = SQRT(ETA2)*COSTH1
RTEMP2 = SQRT(ETA1)*COSTH2
R12 = (RTEMP1 - RTEMP2) / (RTEMP1 + RTEMP2)
RTEMP1 = SQRT(ETA3)*COSTH2
RTEMP2 = SQRT(ETA2)*COSTH3
R23 = (RTEMP1 - RTEMP2) / (RTEMP1 + RTEMP2)
ALPHA2 = SQRT(ETA2) * W / C
RTEMP1 = CDEXP( II*2.0D0*ALPHA2*COSTH2*D )
R = (R12 + R23 * RTEMP1) / (1.0D0 + R12 * R23 * RTEMP1)
RR(I) = CDABS(R)
10 CONTINUE
C      DO I = 1, 51
C      PRINT *, I, THETA(I), RR(I)
C      ENDDO
CALL J06VAF(1,6)
CALL XXXXXX
CALL J06WAF
CALL J06WBF(0.0D0, 90.0D0, 0.0D0, 1.0D0, 1)
CALL J06AHF('OPTICAL REFLECTANCE')
CALL J06AJF(1,'ANGLE')
CALL J06AJF(2,'REFLECTANCE')
CALL J06AAF
IFAIL = 0
CALL J06BAF(THETA, RR, NPTS, 1, 2, IFAIL)
CALL J06WZF
END

and in 3D...

PROGRAM REFL3D
PARAMETER(NPTS = 50, NSLICES = 50)
COMPLEX*16 ETA1, ETA2, ETA3, COSTH2, COSTH3
COMPLEX*16 RTEMP1, RTEMP2, R12, R23, R, ALPHA2
REAL*8 D, W, C, Z(NSLICES,NPTS), THETA1, COSTH1, CONVERT
REAL*8 FSTART, FRANGE
PARAMETER(CONVERT = 3.141592D0 / 180.0D0)
PARAMETER(FSTART = 1.0D12, FRANGE = 1.2D16)

ETA1 = DCMPLX(1.5D0**2, 0.0D0)
ETA2 = DCMPLX(0.044D0, 2.4D0)**2
ETA3 = DCMPLX(1.0D0, 0.0D0)
D    = 40.0D-9
C    = 3.0D8

DO 20 J = 1, NSLICES
W = FSTART +  FRANGE * (J-1) / (NSLICES-1)
DO 10 I=1, NPTS
THETA1 = CONVERT * 90.D0 * (I-1) / (NPTS-1)
COSTH1 = DCOS(THETA1)
COSTH2 = SQRT(1.0D0 -(ETA1/ETA2) * (DSIN(THETA1)**2))
COSTH3 = SQRT(1.0D0 -(ETA1/ETA3) * (DSIN(THETA1)**2))
RTEMP1 = SQRT(ETA2)*COSTH1
RTEMP2 = SQRT(ETA1)*COSTH2
R12    = (RTEMP1 - RTEMP2) / (RTEMP1 + RTEMP2)
RTEMP1 = SQRT(ETA3)*COSTH2
RTEMP2 = SQRT(ETA2)*COSTH3
R23    = (RTEMP1 - RTEMP2) / (RTEMP1 + RTEMP2)
ALPHA2 = SQRT(ETA2) * W / C
RTEMP1 = CDEXP( DCMPLX(0.0D0, 1.0D0) *
+                           2.0D0*ALPHA2*COSTH2* D)
R = (R12 + R23 * RTEMP1)/(1.0D0 + R12 * R23 * RTEMP1)
Z(J,I) = CDABS(R)
10       CONTINUE
20    CONTINUE

CALL J06VAF(1,6)
CALL XXXXXX
CALL J06WAF
CALL J06WBF(0.0D0, 1.0D0, 0.0D0, 1.0D0, 1)
CALL J06AHF('REFLECTANCE')
IFAIL = 0
IVIEW = 1
CALL J06HAF(Z, NSLICES, NPTS, IVIEW, IFAIL)
CALL J06WZF
END
``````

And then to implement a simple finite difference method to solve the Nonlinear Shrodinger equation. (I don’t think I wrote this though – I probably inherited
it from Hiren Mehta)

``````
program soliton
implicit real*8(a-h, o-z)
Parameter(m = 1000, n = 200, nsize = 50, pi = 3.14152)
real*8 U(nsize, nsize)
real*8 phi(nsize, nsize)
real*8 chirp(nsize, nsize)
real*8 R(0:m,0:n),S(0:m,0:n)
complex*16 ii,comu

ii = dcmplx(0.0,1.0)

RANGET = 12.0
vrange = 3.0
c      RANGEZ = pi / 2.0d0
RANGEZ = 0.2
Dt     = ranget/n
Dz     = rangez/m
print *,'dz is ' ,dz
print *,'dt is ' ,dt
print *,'dz/dt**2 is ', dz/dt**2

T0    = 1

print *, 'Enter gamma '
print *, 'Enter beta2 '

DO j=1,N
T      =(j-n/2.0)*Dt
comu   = exp(-t**2/(2.0*t0**2))
R(0,j) = dreal(comu)
S(0,j) = dimag(comu)
R(1,j) = R(0,j)
S(1,j) = S(0,j)
enddo

DO i=1,M-1
DO j=2,N-1
R(i+1,j) = R(i-1,j) +
&               beta2*Dz*( (S(i,j+1)-2.0*S(i,j)+S(i,j-1))/DT**2 )
&              -2.0*Dz*gamma*( R(i,j)**2+S(i,j)**2 ) * S(i,j)
S(i+1,j) = S(i-1,j) -
&               beta2*Dz*( (R(i,j+1)-2.0*R(i,j)+R(i,j-1))/DT**2 )
&              +2.0*Dz*gamma*( R(i,j)**2+S(i,j)**2 ) * R(i,j)
enddo
R(i+1,1) = R(i+1,n-1)
R(i+1,n) = R(i+1,2)
S(i+1,1) = S(i+1,n-1)
S(i+1,n) = S(i+1,2)
enddo

itskip = n / nsize
izskip = m / nsize
is = 75
itskip =1
do i = 1, nsize
do j = 1, nsize
u(i,j) = S(i*izskip,is+j*itskip)**2+R(i*izskip,is+j*itskip)**2
phi(i,j) = atan2(S(i*izskip,is+j*itskip),
&                     R(i*izskip,is+j*itskip))
enddo
do j = 2, nsize-1
chirp(i,j) = (phi(i,j+1)-phi(i,j-1))/(2.0*itskip*dt)
enddo
chirp(i,1) = chirp(i,2)
chirp(i,nsize) = chirp(i,nsize-1)
enddo

20    print *, 'Plot (1)intensity, (2)phase or (3) chirp: (4) to finish'

if (ians.eq.1) then
call plot(U, 'U(z,t)    ' , nsize)
elseif (ians.eq.2) then
call plot(phi, 'Phi(z,t)  ', nsize)
elseif (ians.eq.3) then
call plot(chirp, 'Chirp(z,t)', nsize)
elseif (ians.eq.4) then
stop
else
print *, ians, 'is not an option'
stop
endif
goto 20
end

subroutine plot(U, title, nsize)
real*8 u(nsize, nsize)
character*10 title

IFAIL=0
IVIEW = 1
call J06vaf(1,6)
Call XXXXXX
Call J06WAF
Call J06WBF(0.0D0,1.0D0,0.0D0,1.0D0,1)
Call J06AJF(1, 'z')
Call J06AJF(2, 't')
Call J06AHF(title)
CALL J06HAF(U, nsize, nsize, IVIEW, IFAIL)
Call J06WZF
end
``````

though I never use NAG for publication plotting as it looked a bit rough – I used GINO

``````
program gplot
parameter (nw=900000,mdim=22500)
real      preal(mdim),pimag(mdim),pabss(mdim),pabs(mdim)
real      w(nw),x(mdim),y(mdim),z(mdim),zz(mdim)
complex   pcmplx(mdim),pcmplx1(mdim)
integer   ifill(mdim),index(mdim)
character title*14,ptype*1,dtype*1

common /info/nx,ny,xmin,xmax,ymin,ymax,ptype,dtype
open (12,file='d:\dprn.dj')

call gino

call getdata(title,mdim,pcmplx,pcmplx1,
&                     preal,pimag,pabss,pabs)

10     print *, 'Enter plot type(s):'
print *, 'surface [s]'
print *, 'contour [c]'
print *, 'graph   [g]'
print *, 'sphere  [p]'

c
c     three dimensional graph-----------------------------------------
c
if (ptype.eq.'s' .or. ptype.eq.'S' .or. ptype.eq.' ') then
call selfunc(i)
if (i.eq.1) call splot(title,mdim,nw,preal,w,x,y,z,zz,ifill)
if (i.eq.2) call splot(title,mdim,nw,pimag,w,x,y,z,zz,ifill)
if (i.eq.3) call splot(title,mdim,nw,pabss,w,x,y,z,zz,ifill)
if (i.eq.4) call splot(title,mdim,nw,pabs,w,x,y,z,zz,ifill)
goto 10
endif
c
c     contour plot ---------------------------------------------------
c
if (ptype.eq.'c'.or.ptype.eq.'C') then
call selfunc(i)
if (i.eq.1) call conpt(title,mdim,nw,preal,w,x,y,z,zz,index)
if (i.eq.2) call conpt(title,mdim,nw,pimag,w,x,y,z,zz,index)
if (i.eq.3) call conpt(title,mdim,nw,pabss,w,x,y,z,zz,index)
if (i.eq.4) call conpt(title,mdim,nw,pabs,w,x,y,z,zz,index)
goto 10
endif
c
c     graph plot ---------------------------------------------------
c
if (ptype.eq.'g'.or.ptype.eq.'G') then
call selfunc(i)
if (i.eq.1) call grapt(title,mdim,pimag,preal)
if (i.eq.2) call grapt(title,mdim,pimag,pimag)
if (i.eq.3) call grapt(title,mdim,pimag,pabss)
if (i.eq.4) call grapt(title,mdim,pimag,pabs)
goto 10
endif
c
c     poincare sphere ----------------------------------------------
c
if (ptype.eq.'p'.or.ptype.eq.'P') then
call poin(mdim,preal,pimag,pabss,pabs)
goto 10
endif

call ginend
close(12)
100 format (a)
end

subroutine getdata(title,mdim,pcmplx,pcmplx1,
&    preal,pimag,pabss,pabs)
real      preal(mdim),pimag(mdim),pabss(mdim),pabs(mdim)
complex   pcmplx(mdim), pcmplx1(mdim),  ptemp(500)
character title*14, filenm*12, filenm1*12,ptype*1,dtype*1
character choice*1

common /info/nx,ny,xmin,xmax,ymin,ymax,ptype,dtype

print *, 'Enter data type(c):'
print *,
print *, 'complex 2d........[c]'
print *, 'real 2d...........[r]'
print *, 'complex 1 slice...[d]'
print *, 'real 1 slice......[s]'
print *, 'stokes parameters [p]'
print *, '2 x complex 2d....[z]'

if (dtype.eq.'z' .or. dtype.eq.'Z') then
print 100, 'Graph plotting of pulse evolution'
print 101, 'Filename 1: '
if (filenm.eq.' ') then
print *,'****no file name entered****'
stop
endif
print 101, 'Filename 2: '
if (filenm1.eq.' ') then
print *,'****no file name entered****'
stop
endif
open(60, file=filenm, status='old', err=999)
open(61, file=filenm1, status='old', err=999)
else
print 100, 'Graph plotting of pulse evolution'
print 101, 'Filename : '
if (filenm.eq.' ') then
print *,'****no file name entered****'
stop
endif
open(20, file=filenm, status='old', err=999)
endif
C
C real one slice
C
if (dtype.eq.'s' .or. dtype.eq.'S') then
dtype = 's'
do j=1,mdim
cccc this is x, y
pabss(j) = preal(j)
pabs(j) = preal(j)
enddo
print *,'****array out of bounds, increase mdim****'
stop
C
C stokes parameters
C
elseif (dtype.eq.'p'.or.dtype.eq.'P')then
dtype = 'p'
c next two values are irrelevent
ymin = 0.0
yinc = 1.0
do j = 1, mdim
enddo
C
C real 2d
C
elseif (dtype.eq.'r'.or.dtype.eq.'R')then
dtype = 'r'
do j=1,mdim
enddo
print *,'****array out of bounds, increase mdim****'
stop
C
C complex 2d
C
elseif (dtype.eq.'c'.or.dtype.eq.'C'.or.dtype.eq.' ')then
dtype = 'c'
print *,'xmin = ',xmin
print *,'xmax = ',xmax
print *,'change x'
if (choice.eq.'y' .or. choice.eq.'Y') then
print *,'new xmin'
print *,'new xmax'
xstep = (xmax-xmin)/nx
nnx =  (xmaxn-xminn)/xstep
nnxmin = (xminn-xmin)/xstep
nnxmax = nx - (xmax-xmaxn)/xstep
do j=1,mdim
iii = 1
do jjj = nnxmin, nnxmax
pcmplx(iii+(j-1)*nnx) = ptemp(jjj)
iii = iii + 1
enddo
enddo
print *,'****array out of bounds, increase mdim****'
stop
endif
do j=1,mdim
enddo
print *,'****array out of bounds, increase mdim****'
stop
C
C 2 x complex 2d
C
elseif (dtype.eq.'z'.or.dtype.eq.'Z')then
dtype = 'z'
do j=1,mdim
enddo
print *,'****array out of bounds, increase mdim****'
stop
do j=1,mdim
enddo
print *,'****array out of bounds, increase mdim****'
stop
endif

7  continue
ny = j - 1
ymax = ymin + yinc * ny
return

5  continue
nx = j - 1
ny = 1
ymin = 0.0
yinc = 0.0
c      do i=1,nx
c      pcmplx(i)=dcmplx(preal(i),pimag(i))
c      enddo
return
c
c  return from decreasing window
c
9     continue

do jjj = nnxmin, nnxmax
pcmplx(jjj+(j-1)*nnx) = ptemp(jjj)
enddo
nx = nnx
xmin = xminn
xmax = xmaxn
10 continue
ny = j - 1
cccc    ymax = ymin + yinc * ny
ymax = ymin + yinc * (ny-1)

do i=1,nx
do j=1,ny
preal(i+(j-1)*nx) = real( pcmplx(i+(j-1)*nx) )
pimag(i+(j-1)*nx) = aimag( pcmplx(i+(j-1)*nx) )
pabss(i+(j-1)*nx) = cabs( pcmplx(i+(j-1)*nx) )**2
pabs(i+(j-1)*nx) = cabs( pcmplx(i+(j-1)*nx) )
enddo
enddo
close(20)
return

89      continue
ny = j - 1
ymax = ymin + yinc * (ny-1)

do i=1,nx
do j=1,ny
preal(i+(j-1)*nx) = real( pcmplx(i+(j-1)*nx) )
&                               + real( pcmplx1(i+(j-1)*nx) )
pimag(i+(j-1)*nx) = aimag( pcmplx(i+(j-1)*nx) )
&                               + aimag( pcmplx1(i+(j-1)*nx) )
pabss(i+(j-1)*nx) = cabs( pcmplx(i+(j-1)*nx) )**2
&                               + cabs( pcmplx1(i+(j-1)*nx) )**2
pabs(i+(j-1)*nx) = cabs( pcmplx(i+(j-1)*nx) )
&                               + cabs( pcmplx1(i+(j-1)*nx) )
enddo
enddo

close(60)
close(61)
100 format(a)
101 format(a,\$)
return
999 print *,'****could not open file ',filenm,'****'
stop
end

subroutine poin(mdim,s0,s1,s2,s3)
real s0(mdim), s1(mdim), s2(mdim), s3(mdim)
character*1 ptype, dtype, face
real x1(4),y1(4),z1(4),x2(4),y2(4),z2(4),x3(4),y3(4),z3(4)
parameter(ab=30.0,ah=60.0)

common /info/nx,ny,xmin,xmax,ymin,ymax,ptype,dtype

data x1/0.0,  ah, -ah, 0.0/
data y1/-ab,  ab,  ab, -ab/
data z1/0.0, 0.0, 0.0, 0.0/

data x2/ab , -ab, -ab, ab /
data y2/0.0,  ah, -ah, 0.0/
data z2/0.0, 0.0, 0.0, 0.0/

data x3/ 0.0, 0.0, 0.0, 0.0/
data y3/-ab , ab,   ab, -ab/
data z3/ 0.0,  ah, -ah, 0.0/

cw = 19.5
ch = 19.5

print *, 'Join every ? points '
print *,
print *,'Display ?[a]'
print *,'Front Face only...(a) '
print *,'Both Faces........(b) '
print *,'Reverse Face only.(c) '
if (face.eq.'a' .or. face.eq.'A') then
face = 'a'
elseif (face.eq.'b' .or. face.eq.'B') then
face = 'b'
elseif (face.eq.'c' .or. face.eq.'C') then
face = 'c'
else
print *, '*****not an option*****'
stop
endif

call printq

qmax = s0(1)
do i = 2, ny
if (qmax.lt.s0(i)) qmax = s0(i)
enddo

scale = 1000.0 / qmax
do i = 1, ny
s0(i) = s0(i)*scale
s1(i) = s1(i)*scale
s2(i) = s2(i)*scale
s3(i) = s3(i)*scale
enddo
call chasiz(cw,ch)

call piccle
call fntwgt(-3)
call chafnt(15)
call linwid(2.0)

data vx, vy, vz/ -1.0, -0.3, -1.0/

call vspher(0.0, 0.0, 0.0, 1400., vx, vy, vz, 10000.0)
call vposit(500.0, 375.0)
call viewup
c
c draw sphere
c
call lincol(2)
call movto3(0.0, 1000.0, 0.0)
call arcto3(0.0, 0.0, 0.0, 0.0,1000.0, 0.0, -1.0, 0.0, 1.0)

call broken(1)
call movto3(1000.0, 0.0, 0.0)
call arcto3(0.0, 0.0, 0.0, 1000.0, 0.0, 0.0, 0.0, 0.0, -1.0)
call broken(0)
c
c draw axis
c
dis = 30.0
call movto3(-1070.0, 0.0, 0.0)
call linto3(1070.0, 0.0, 0.0)
call movto3(1190.0, 0.0, 0.0)
call chastr('s1')
call movto3(1070.0, 0.0, 0.0)
call pofby3(0,0,0,x1,y1,z1,4)

call movto3(0.0, -1070.0, 0.0)
call linto3(0.0, 1070.0, 0.0)
call movto3(-dis, 1150.0, dis)
call chastr('s2')
call movto3(0.0, 1070.0, 0.0)
call pofby3(0,0,0,x2,y2,z2,4)

call movto3(0.0, 0.0, -1070.0)
call linto3(0.0, 0.0, 1070.0)
call movto3(-dis, 0.0, 1290.0)
call chastr('s3')
call movto3(0.0, 0.0, 1070.0)
call pofby3(0,0,0,x3,y3,z3,4)
c
c plot stokes vectors
c
call lincol(1)
open (33, file = 'noplot')
eps = 0.2
noplot = 0
do i = 1, ny, njoin
do j = 1, njoin-1
call movto3(s1(i+j-1),s2(i+j-1),s3(i+j-1))
vp = vx*s1(i+j) + vy*s2(i+j) + vz*s3(i+j)
iflag = 0
if (abs(s1(i+j-1)-s1(i+j)).gt.eps) iflag = iflag + 1
if (abs(s2(i+j-1)-s2(i+j)).gt.eps) iflag = iflag + 1
if (abs(s3(i+j-1)-s3(i+j)).gt.eps) iflag = iflag + 1
if (iflag.gt.1) then
if (vp.lt.0.0) then
if (face.ne.'c')
&            call arcto3(0.0,0.0,0.0,s1(i+j),s2(i+j),s3(i+j),
&                                   s1(i+j),s2(i+j),s3(i+j))
elseif(face.eq.'b') then
call broken(1)
call arcto3(0.0,0.0,0.0,s1(i+j),s2(i+j),s3(i+j),
&                                   s1(i+j),s2(i+j),s3(i+j))
call broken(0)
endif
else
noplot = noplot + 1
write(33,99) s1(i+j-1), s2(i+j-1), s3(i+j-1)
write(33,99) s1(i+j), s2(i+j), s3(i+j)
write(33,99)
call dot(0.0)
endif
enddo
enddo
99      format (2x,3(g14.5,2x))

call get_key@(k)
call devend
if (noplot.gt.0) then
print *,'Number of points not plotted ', noplot
call get_key@(k)
endif
close(33)
100     format (a)
end

subroutine grapt(title,mdim,x,y)
character title*14, rev*1, ptype*1, dtype*1
character dmp*1
real x(mdim),y(mdim)
real p(2500), q(2500)

common /info/nx,ny,xmin,xmax,ymin,ymax,ptype,dtype
c      open (13,file='e:\hirdmp.dat')
open (13,file='hirdmp.dat')

c      print *,'number of plotting points '
c      print *,' this is ',xn

print *,'plotting routine '
print *,'lines  (1)'
print *,'curves (2)'
print *,'spline (3)'

c      step = nx/xn
c      nn = int (xn)
nn = nx
step = 1
if (dtype.eq.'d' .or. dtype.eq.'s') then
do i=1, nn
p(i) = x(int(i*step))
q(i) = y(int(i*step))
enddo
else
print *,'reverse surface ?(n)'
print *,'slice number ?'
print *,'nx ,ny ',nx,ny

if (rev.eq.'y'.or. rev.eq.'Y') then
nn = ny
if (nslice.gt.nx) then
print *,'slice number too big '
stop
endif
do i=1,nn
p(i)=ymin+(i-1)*(ymax-ymin)/(ny-1)
q(i)=y(nslice+(i-1)*nx)
enddo
elseif (rev.eq.'n'.or. rev.eq.'N' .or. rev.eq.' ') then
if (nslice.gt.ny) then
print *,'slice number too big '
stop
endif
do i=1,nn
p(i)=xmin+(i-1)*(xmax-xmin)/(nx-1)
q(i)=y(i+(nslice-1)*nx)
enddo
endif
endif

cw = 19.5
ch = 19.5

call printq

qmax = q(1)
qmin = q(1)
do i=1, nn
if (q(i).gt.qmax) qmax = q(i)
if (q(i).lt.qmin) qmin = q(i)
enddo

call chasiz(cw,ch)
qmin = qmin-(qmax-qmin)/5.0
call axisca(3,5,qmin,qmax,2)

call linwid(2.0)
call piccle
call fntwgt(-3)
call chafnt(15)
c        call chasiz(19.0,20.0)
call axnstr('time (nsec)',20.0,1,0)
call axnstr('intensity',70.0,2,0)
c      call graph (q,p,nn,0,isect,0)
call graph (p,q,nn,0,isect,0)

call get_key@(k)
call devend
print *,'dump slice ?(n)'
if (dmp.eq.'y'.or. dmp.eq.'Y') then
do i=1,nn
write(13,*) p(i), q(i)
enddo
endif
100   format(a)
close(13)
end

subroutine splot(title,mdim,nw,fdata,w,x,y,z,zz,ifill)
character title*14,centre@*14,gtype*1,ctype*1,ttype*1
character ptype*1, dtype*1
character*14 xtit,ytit,ztit
real fdata(mdim),x(mdim),y(mdim),z(mdim),zz(mdim),w(nw)
integer ifill(mdim)

common /info/nx,ny,xmin,xmax,ymin,ymax, ptype, dtype

call getcol(ctype)

print *,'default grid size ?(y) '
if (gtype.eq.'n' .or. gtype.eq.'N') then
print *,'enter gridsize x x y'
else
kx = nx
ky = ny
endif

print *,'default grid lines ?(y) '
if (gtype.eq.'n' .or. gtype.eq.'N') then
print *,'enter steps at which contour lines are marked'
else
nxstep = 1
nystep = 1
endif

print *,'type of section lines ? '
print *,'parallel to both axes (1)'
print *,'parallel to x axes (2)'
print *,'parallel to y axes (3)'

print *,'view number 0 - 3 ? '

if (iview.gt.3 .or. iview.lt.0) then
print *,'****funtion not an option****'
stop
endif

print *,'title ?[n]'
if (ttype.eq.'y' .or. ttype.eq.'Y') then
else
title= ' '
endif

print *,'labelled axes ?[n]'
if (ttype.eq.'y' .or. ttype.eq.'Y') then
print *,'x - title '
print *,'y - title '
print *,'z - title '
else
xtit = ' '
ytit = ' '
ztit = ' '
endif
c      xtit = 'distance (cm)'
c      ytit = 'time (nsec)'
c      ztit = 'intensity'

call printq
call setcol(ctype)

do i=1,nx
do j=1,ny
x(i+(j-1)*nx)=xmin+(i-1)*(xmax-xmin)/(nx-1)
y(i+(j-1)*nx)=ymin+(j-1)*(ymax-ymin)/(ny-1)
z(i+(j-1)*nx)=fdata(i+(j-1)*nx)
enddo
enddo

np = nx*ny
call rangrd(np,x,y,z,kx,xmin,xmax,ky,ymin,ymax,zz,nw,w)

zmin=zz(1)
zmax=zz(1)
do i=1,kx
do j=1,ky
if (zz(i+(j-1)*kx).gt.zmax) zmax=zz(i+(j-1)*kx)
if (zz(i+(j-1)*kx).lt.zmin) zmin=zz(i+(j-1)*kx)
enddo
enddo
c       zmax = 1.2

ioff = 11
call colinf(ncol,nint)
call colenq(icol)
if (ncol.le.1) ioff = 1001
scale = (ncol-11)/(zmax - zmin)
do i=1,kx-1
do j=1,ky-1
ifill(i+(j-1)*(kx-1)) = int((zz(i+(j-1)*kx)-zmin)*scale)+ioff
enddo
enddo

if (isect.eq.1) call prjsec(0)
if (isect.eq.2) call prjsec(1)
if (isect.eq.3) call prjsec(2)

call piccle
ccccccccccccccccccccccccccccccc
c        call chafnt(15)
c        call fntwgt(0)
call chafnt(100)
call fntwgt(3)
call linwid(2.0)

call chaenq(itype,width,height,nsize,slant,angle)
call prjgrd(nxstep,nystep)
if (ctype.eq.'n'.or. ctype.eq.'N') then
call isoprj(kx,xmin,xmax,ky,ymin,ymax,zz,iview,nw,w)
else
call isofil(kx,xmin,xmax,ky,ymin,ymax,zz,iview,ifill,nw,w)
endif

xmid = xmin + (xmax-xmin)/2.0
ymid = ymin + (ymax-ymin)/2.0
zmid = zmin + (zmax-zmin)/2.0
if (iview.eq.0) then
call isospa(xmax,ymax,zmax,xx,yy)
call chasiz(20.0,24.0)
call movto2(xx+30,yy+20.0)
call chajus(0)
call chastr(centre@(title,14))
call chasiz(8.0,10.0)
call isospa(xmin,ymax,zmid,xx,yy)
call movto2(xx-110.0,yy)
call chaang(90.0)
call chastr(centre@(ztit,14))
call isospa(xmid,ymin,zmin,xx,yy)
call movto2(xx+70.0,yy-50.0)
call chaang(30.0)
call chastr(centre@(xtit,14))
call isospa(xmin,ymid,zmin,xx,yy)
call movto2(xx-80.0,yy-60.0)
call chaang(-30.0)
call chastr(centre@(ytit,14))
endif
if (iview.eq.1) then
call isospa(xmin,ymin,zmax,x1,yy)
call isospa(xmax,ymax,zmax,x2,yy)
call isospa(xmin,ymax,zmax,xx,yy)
call chasiz(20.0,24.0)
call movto2(x1+0.5*(x2-x1)+20.0,yy+20.0)
call chajus(0)
call chastr(centre@(title,14))
call chasiz(8.0,10.0)
call isospa(xmin,ymin,zmid,xx,yy)
call movto2(xx-110.0,yy)
call chaang(90.0)
call chastr(centre@(ztit,14))
call isospa(xmax,ymid,zmin,xx,yy)
call movto2(xx+70.0,yy-50.0)
call chaang(30.0)
call chastr(centre@(ytit,14))
call isospa(xmid,ymin,zmin,xx,yy)
call movto2(xx-50.0,yy-60.0)
call chaang(-30.0)
call chastr(centre@(xtit,14))
endif
if (iview.eq.2) then
call isospa(xmin,ymin,zmax,xx,yy)
call chasiz(20.0,24.0)
call movto2(xx+30,yy+20.0)
call chajus(0)
call chastr(centre@(title,14))
cccccccccccc           call chasiz(8.0,10.0)
call chasiz(12.0,15.0)
call isospa(xmax,ymin,zmid,xx,yy)
call movto2(xx-80.0,yy)
call chaang(90.0)
call chastr(centre@(ztit,14))
call isospa(xmid,ymax,zmin,xx,yy)
call movto2(xx+60.0,yy-40.0)
call chaang(30.0)
call chastr(centre@(xtit,14))
call isospa(xmax,ymid,zmin,xx,yy)
call movto2(xx-60.0,yy-40.0)
call chajus(0)
call chaang(-30.0)
call chastr(centre@(ytit,14))
endif

call chajus(-1)
call chasiz(width,height)

call get_key@(k)
call devend
100 format(a)
end

subroutine conpt(title,mdim,nw,fdata,w,x,y,z,zz,index)
character title*14,centre@*14,ctype*1,ttype*1
character*14 xtit,ytit
character ptype*1,dtype*1
real fdata(mdim),x(mdim),y(mdim),z(mdim),zz(mdim),w(nw)
integer index(mdim)

common /info/nx,ny,xmin,xmax,ymin,ymax, ptype, dtype

call getcol(ctype)

ism = 0
print *,'number of contours ?(10)'

print *,'title ?[n]'
if (ttype.eq.'y' .or. ttype.eq.'Y') then
else
title=' '
xtit = ' '
ytit = ' '
endif

print *,'labelled axes ?[n]'
if (ttype.eq.'y' .or. ttype.eq.'Y') then
print *,'x - title '
print *,'y - title '
endif
c      ytit = 'distance (cm)'
c      xtit = 'time (nsec)'

call printq
call setcol(ctype)

call consca(0.0,0.0,0.0,1)
call chafnt(15)
call chaenq(itype,width,height,nsize,slant,angle)

ioff = 11
call colinf(ncol,nint)
call colenq(icol)
if (ncol.le.1) ioff = 1001
scale = (ncol-11)/(ncont+1)
do i=1,ncont+1
index(i)=int(i*scale)+ioff
enddo

do i=1,nx
do j=1,ny
x(i+(j-1)*nx)=xmin+(i-1)*(xmax-xmin)/(nx-1)
y(i+(j-1)*nx)=ymin+(j-1)*(ymax-ymin)/(ny-1)
z(i+(j-1)*nx)=fdata(i+(j-1)*nx)
enddo
enddo

zmin=0.0
zmax=0.0
do i=1,nx
do j=1,ny
if (z(i+(j-1)*nx).gt.zmax) zmax=z(i+(j-1)*nx)
if (z(i+(j-1)*nx).lt.zmin) zmin=z(i+(j-1)*nx)
enddo
enddo

call piccle
call levels(zmin,zmax)
if (ctype.eq.'n' .or. ctype.eq.'N') then
call dracon(nx,xmin,xmax,ny,ymin,ymax,fdata,ncont,ism,nw,w)
else
call filcon(nx,xmin,xmax,ny,ymin,ymax,fdata,ncont,ism,
*            index,1,nw,w)
endif

ymid = ymin + (ymax-ymin)/2.0
xmid = xmin + (xmax-xmin)/2.0

xth = xmid + xmid/2.7

call chasiz(20.0,24.0)
call conspa(xth,ymax,xx,yy)
call movto2(xx,yy+80.0)
call chajus(0)
c      call chastr(centre@(title,14))

call conspa(xmin,ymid,xx,yy)
call movto2(xx-100.0,yy)
call chaang(90.0)
call chajus(0)
call chasiz(12.0,14.0)
call chastr(centre@(xtit,14))
call chaang(0.0)
call conspa(xmid,ymin,xx,yy)
call movto2(xx,yy-50.0)
call chastr(centre@(ytit,14))
call chajus(-1)
call chasiz(width,height)

call get_key@(k)
call devend
100 format(a)
end

subroutine selfunc(ichoice)

print *,'enter function to plot '
print *,'1 - real'
print *,'2 - imaginary'
print *,'3 - absolute value squared'
print *,'4 - absolute value'

if (ichoice.gt.4 .or. ichoice.lt.1) then
print *,'****funtion not an option****'
stop
endif
end

subroutine printq

character choice*1
print *,'output to printer ?(n)'
if (choice.eq.'n' .or. choice.eq.'N' .or. choice.eq.' ') then
call vga
call papenq(x,y,ipapty)
call setvp2(0.0,1000.0,0.0,750.0,0.0,x,0.0,y)
c           call rgbdef(0,0.0,0.0,0.0)
return
endif
print *,'which printer ?'
print *,'0 - epson'
print *,'1 - hp deskjet colour'
print *,'2 - hp deskjet black'
print *,'3 - hp 7550 plotter '
print *,'4 - ps'
if (ichoice.gt.4 .or. ichoice.lt.0) then
print *,'****funtion not an option****'
stop
endif

if (ichoice.eq.0) then
c          call fx80
call lq850
call papenq(xpap,ypap,i)
call device(12,0)
c          call device(1,1)
call devpap(ypap,xpap,0)
call papenq(xpap,ypap,i)
call setvp2(0.0,1000.0,0.0,750.0,0.0,xpap,0.0,ypap)
call vptswi(0)
call colenq(icol)
return
endif
if (ichoice.eq.1) then
call dj550c
call papenq(xpap,ypap,i)
call device(12,0)
call devpap(ypap,xpap,0)
call papenq(xpap,ypap,i)
call setvp2(0.0,1000.0,0.0,750.0,0.0,xpap,0.0,ypap)
call vptswi(0)
call colenq(icol)
return
endif
if (ichoice.eq.2) then
call dj500m
call papenq(xpap,ypap,i)
call device(12,0)
call devpap(ypap,xpap,0)
call papenq(xpap,ypap,i)
call setvp2(0.0,1000.0,0.0,750.0,0.0,xpap,0.0,ypap)
call vptswi(0)
call colenq(icol)
return
endif

if (ichoice.eq.3) then
c         call hp7550
call papenq(xpap,ypap,i)
call setvp2(0.0,1000.0,0.0,750.0,0.0,xpap,0.0,ypap)
call device(12,0)
endif
if (ichoice.eq.4) then

c         call eps(0,10.0,10.0,1000.0,750.0,1000.0,750.0)
c          call papenq(xpap,ypap,i)
c          call setvp2(0.0,1000.0,0.0,750.0,0.0,xpap,0.0,ypap)

call eps(0,0.0,0.0,1000.0,750.0,50.0,50.0)
call papenq(x,y,i)
factor = min(y/1000.0, x/750.0)
call devpap(x,y,0)
call units(factor)
c         call setvp2(0.0,1000.0,0.0,750.0,0.0,xpap,0.0,ypap)
call device(12,0)
endif
100 format(a)
end

subroutine getcol(ctype)
character ctype*1
print *,'colour map ?(a) '
print *,'red[r]'
print *,'blue[b]'
print *,'all available[a]'
print *,'no colour map[n]'
100 format(a)
end

subroutine setcol(ctype)
character ctype*1
call colinf(ncol,nint)
call colenq(icol)
if (ctype.eq.'a' .or. ctype.eq.'A') then
do i=11,ncol
z = 30.0 + 330.0/(ncol-11)*(i-11)
call hsvdef(i,z,1.0,1.0)
enddo
endif

if (ctype.eq.'r' .or. ctype.eq.'R') then
do i=11,ncol
z = 120.0/(ncol-11)*(i-11)
call hsvdef(i,z,1.0,1.0)
enddo
endif

if (ctype.eq.'b' .or. ctype.eq.'B') then
do i=11,ncol
z = 320.0-160.0/(ncol-11)*(i-11)
call hsvdef(i,z,1.0,1.0)
enddo
endif
call hsvdef(10,0.0,0.0,0.0)
end
``````