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 '
          read *, gamma
          print *, 'Enter beta2 '
          read *, 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'
          read *, ians

          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]'
            read 100, ptype
            
    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]'
            read 100, dtype
            
            if (dtype.eq.'z' .or. dtype.eq.'Z') then
              print 100, 'Graph plotting of pulse evolution'
              print 101, 'Filename 1: '
              read 100, filenm
              if (filenm.eq.' ') then
                 print *,'****no file name entered****'
                 stop
              endif
              print 101, 'Filename 2: '
              read 100, filenm1
              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 : '
              read 100, filenm
              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
                 read(20,*,end=5) preal(j), pimag(j)
                    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
                 read(20,*,end=7) preal(j),pimag(j),pabss(j),pabs(j)
              enddo
    C
    C real 2d
    C
            elseif (dtype.eq.'r'.or.dtype.eq.'R')then
              dtype = 'r'
              read(20,100) title
              read(20,*) nx,xmin,xmax,ymin,yinc
              do j=1,mdim
                read(20,*,end=7) (preal(i+(j-1)*nx),i=1,nx)
              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'
              read(20,100) title
              read(20,*) nx,xmin,xmax,ymin,yinc
              print *,'xmin = ',xmin
              print *,'xmax = ',xmax
              print *,'change x'
              read 100,choice
              if (choice.eq.'y' .or. choice.eq.'Y') then
                print *,'new xmin'
                read *,xminn
                print *,'new xmax'
                read *,xmaxn
    cccc        dummy read
    ccccccc            read(20,*) nn2,xm2,xx2
                xstep = (xmax-xmin)/nx
                nnx =  (xmaxn-xminn)/xstep 
                nnxmin = (xminn-xmin)/xstep
                nnxmax = nx - (xmax-xmaxn)/xstep
                do j=1,mdim
                  read(20,*,end=9) (ptemp(i),i=1,nx)
                  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
    cccc        dummy read
    ccccccc            read(20,*) nn2,xm2,xx2
                do j=1,mdim
                    read(20,*,end=10) (pcmplx(i+(j-1)*nx),i=1,nx)
                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'
              read(60,100) title
              read(60,*) nx,xmin,xmax,ymin,yinc
    cccc        dummy read
    ccccccc          read(60,*) nn2,xm2,xx2
                do j=1,mdim
                    read(60,*,end=88) (pcmplx(i+(j-1)*nx),i=1,nx)
                enddo
                print *,'****array out of bounds, increase mdim****'
                stop
    88        read(61,100) title
              read(61,*) nx,xmin,xmax,ymin,yinc
    cccc        dummy read
    cccccc          read(61,*) nn2,xm2,xx2
                do j=1,mdim
                    read(61,*,end=89) (pcmplx1(i+(j-1)*nx),i=1,nx)
                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 '
            read *, njoin
            print *,
            print *,'Display ?[a]'
            print *,'Front Face only...(a) '
            print *,'Both Faces........(b) '
            print *,'Reverse Face only.(c) '
            read 100, face
            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      read *,xn
    c      print *,' this is ',xn
            
            print *,'plotting routine '
            print *,'lines  (1)'
            print *,'curves (2)'
            print *,'spline (3)'
            read *, isect
            
    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)'
              read 100,rev
              print *,'slice number ?'
              print *,'nx ,ny ',nx,ny
              read *,nslice
            
              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)'
            read 100,dmp
            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) '
            read 100, gtype
            if (gtype.eq.'n' .or. gtype.eq.'N') then
               print *,'enter gridsize x x y'
               read *,kx,ky
            else
               kx = nx
               ky = ny
            endif
            
            print *,'default grid lines ?(y) '
            read 100, gtype
            if (gtype.eq.'n' .or. gtype.eq.'N') then
               print *,'enter steps at which contour lines are marked'
               read *,nxstep,nystep
            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)'
            read *, isect

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

            print *,'title ?[n]'
            read 100, ttype
            if (ttype.eq.'y' .or. ttype.eq.'Y') then
               read 100,title
            else
               title= ' '
            endif
            
            print *,'labelled axes ?[n]'
            read 100, ttype
            if (ttype.eq.'y' .or. ttype.eq.'Y') then
               print *,'x - title '
               read 100,xtit
               print *,'y - title '
               read 100,ytit
               print *,'z - title '
               read 100,ztit
            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)'
            read *,ncont
            
            print *,'title ?[n]'
            read 100, ttype
            if (ttype.eq.'y' .or. ttype.eq.'Y') then
               read 100,title
            else
               title=' '
               xtit = ' '
               ytit = ' '
            endif
            
            print *,'labelled axes ?[n]'
            read 100, ttype
            if (ttype.eq.'y' .or. ttype.eq.'Y') then
               print *,'x - title '
               read 100,xtit
               print *,'y - title '
               read 100,ytit
            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'

            read *, ichoice
            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)'
            read 100,choice
            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'
            read *, ichoice
            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]'
            read 100, ctype
      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
Advertisements

Leave a Comment »

No comments yet.

RSS feed for comments on this post. TrackBack URI

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

Create a free website or blog at WordPress.com.

%d bloggers like this: