home *** CD-ROM | disk | FTP | other *** search
/ ARM Club 3 / TheARMClub_PDCD3.iso / hensa / maths / pgplot_1 / Examples / f77 / PGDEMO8 < prev    next >
Text File  |  1996-05-14  |  8KB  |  314 lines
  1.       PROGRAM PGDEM8
  2. C
  3. C From: Philip Palmer <philip@maths.qmw.ac.uk>
  4. C Date: Mon, 26 Nov 90 12:27:22 GMT
  5. C
  6. C This program plots a 3d surface using PGSURF.
  7. C
  8.       INTEGER NCURVE
  9.       PARAMETER(NCURVE=21)
  10.       INTEGER I,J, PGBEG
  11.       REAL A(NCURVE,NCURVE),XLIMS(3,2)
  12.       REAL AA1,AA2,AA3,AA4
  13.       REAL X,Y,DX,DY,THD,PHD,Q
  14.       REAL POT
  15.       COMMON/COEFS/AA1,AA2,AA3,AA4
  16. C
  17.       THD=40.
  18.       PHD=-35.
  19.       Q=0.1
  20.       AA1=60./17.-24.*Q/17.
  21.       AA2=12./17.
  22.       AA3=52./17.
  23.       AA4=24./17.
  24.       XLIMS(1,1)=0.
  25.       XLIMS(1,2)=0.43
  26.       XLIMS(2,1)=-1.5
  27.       XLIMS(2,2)=1.5
  28.       XLIMS(3,1)=0.
  29.       XLIMS(3,2)=1.
  30.       DX=(XLIMS(1,2)-XLIMS(1,1))/(NCURVE-1)
  31.       DY=(XLIMS(2,2)-XLIMS(2,1))/(NCURVE-1)
  32.       DO 1 I=1,NCURVE
  33.           X=XLIMS(1,1)+(I-1)*DX
  34.           DO 2 J=1,NCURVE
  35.               Y=XLIMS(2,1)+(J-1)*DY
  36.               A(I,J)=POT(X,Y)
  37.     2     CONTINUE
  38.     1 CONTINUE
  39.       IF (PGBEG(0,'?',1,1) .NE. 1) STOP
  40.       CALL PGSURF(A,NCURVE,NCURVE,XLIMS,THD,PHD)
  41.       CALL PGEND
  42.       END
  43.  
  44.       REAL FUNCTION POT(X,Y)
  45.       REAL X,Y
  46. C
  47.       REAL X2,Y2,XX,YY
  48.       REAL AA1,AA2,AA3,AA4
  49.       COMMON/COEFS/AA1,AA2,AA3,AA4
  50. C
  51.       Y2=Y*Y
  52.       X2=12.*X*X
  53.       XX=.5*(Y2+X2)
  54.       YY=Y2/3.+Y*X2
  55.       POT=1.-AA1*XX-AA2*YY+AA3*XX*XX+AA4*XX*YY
  56.       RETURN
  57.       END
  58.  
  59.       SUBROUTINE PGSURF(A,NX,NY,XLIMS,THD,PHD)
  60.       INTEGER NX,NY
  61.       REAL A(NX,NY),XLIMS(3,2),THD,PHD
  62. C
  63. C This routine plots a 3d surface projected onto a 2D plane.
  64. C The underside of the surface appears dotted or in blue, on
  65. C clour terminals. This routine does the projection for you,
  66. C you just need to specify the viewing direction in terms of
  67. C spherical polar angles.  As a result, this routine calls 
  68. C pgwind for you.
  69. C
  70. C Arguments:
  71. C  a (input)     : data array, equally spaced in x and equally
  72. C          spaced in y, although steps in x and y may
  73. C          differ.
  74. C  nx (input)    : dimension of data array in x direction.
  75. C  ny (input)    : dimension of data array in y direction.
  76. C  xlims (input): array of min and max values in x, y and z
  77. C          respectively.
  78. C  thd (input)     : viewing direction given as spherical theta angle
  79. C          in degrees. 0<= thd <= 90.
  80. C  phd (input)    : viewing direction given as spherical phi angle
  81. C          in degrees. -180 <= phd <= 180.
  82. C
  83.       INTEGER I,J,K,N,IFLAG,I0,I1,NSTART,NLINE,ISGN,JSGN
  84.       INTEGER IMAX,JMAX,K0,K1
  85.       REAL DTOR, TH,PH,CTH,STH,CPH,SPH
  86.       REAL XMIN,XMAX,YMIN,YMAX,ZMIN,ZMAX,Z1,Z2,Q1,Q2,Q3
  87.       REAL A1,A2,A3,A4,B1,B2,B3,B4,XSTART,YSTART
  88.       REAL XPMIN,XPMAX,YPMIN,YPMAX,DX,DY,ZZ,T
  89.       REAL XA(2),YA(2),ZA(2)
  90.       COMMON/PGSRF1/CTH,STH,CPH,SPH,NLINE,NSTART,IFLAG
  91. C
  92.       IFLAG=0
  93.       CALL PGQCOL(I0,I1)
  94.       IF(I1.NE.1) IFLAG=1
  95.       CALL PGBBUF
  96.       DTOR=0.0174533
  97.       TH=THD*DTOR
  98.       PH=PHD*DTOR
  99.       CTH=COS(TH)
  100.       STH=SIN(TH)
  101.       CPH=COS(PH)
  102.       SPH=SIN(PH)
  103.       XMIN=XLIMS(1,1)
  104.       XMAX=XLIMS(1,2)
  105.       YMIN=XLIMS(2,1)
  106.       YMAX=XLIMS(2,2)
  107.       ZMIN=XLIMS(3,1)
  108.       ZMAX=XLIMS(3,2)
  109. C
  110. C Calculate plotting order from
  111. C viewing orientation
  112. C
  113.       Z1=ZMIN*STH
  114.       Z2=ZMAX*STH
  115.       IF(PHD.GT.0.) THEN
  116.         A1=-XMAX*SPH
  117.         A2=-XMIN*SPH
  118.         B3=-YMAX*SPH*CTH
  119.         B4=-YMIN*SPH*CTH
  120.         YSTART=YMAX
  121.         ISGN=-1
  122.       ELSE
  123.         A1=-XMIN*SPH
  124.         A2=-XMAX*SPH
  125.         B3=-YMIN*SPH*CTH
  126.         B4=-YMAX*SPH*CTH
  127.         YSTART=YMIN
  128.         ISGN=1
  129.       ENDIF
  130.       IF(ABS(PHD).LT.90.) THEN
  131.         B1=YMIN*CPH
  132.         B2=YMAX*CPH
  133.         A3=-XMAX*CPH*CTH
  134.         A4=-XMIN*CPH*CTH
  135.         XSTART=XMAX
  136.         JSGN=-1
  137.       ELSE
  138.         B1=YMAX*CPH
  139.         B2=YMIN*CPH
  140.         A3=-XMIN*CPH*CTH
  141.         A4=-XMAX*CPH*CTH
  142.         XSTART=XMIN
  143.         JSGN=1
  144.       ENDIF
  145.       XPMIN=MIN(A1,B1)
  146.       XPMAX=MAX(A2,B2)
  147.       YPMAX=MAX(A4,B4,Z2)
  148.       YPMIN=MIN(A3,B3,Z1)
  149.       CALL WINDOW(XPMIN,XPMAX,YPMIN,YPMAX)
  150. C
  151. C Draw coordinate axes
  152. C
  153.       NSTART=0
  154.       NLINE=1
  155.       XA(1)=0.
  156.       XA(2)=0.
  157.       YA(1)=YMIN
  158.       YA(2)=YMAX
  159.       CALL PROJ(XA,YA,XA,1,1)
  160.       NLINE=1
  161.       YA(1)=XMIN
  162.       YA(2)=XMAX
  163.       CALL PROJ(YA,XA,XA,1,1)
  164.       NLINE=1
  165.       YA(1)=ZMIN
  166.       YA(2)=ZMAX
  167.       CALL PROJ(XA,XA,YA,1,1)
  168. C
  169. C Draw curves stepped in x
  170. C
  171.       IMAX=NX-1
  172.       JMAX=NY-1
  173.       DX=JSGN*(XMAX-XMIN)/IMAX
  174.       DY=(YMAX-YMIN)/JMAX
  175.       Q1=DX*STH*SPH
  176.       Q2=DY*STH*CPH
  177.       Q3=DX*DY*CTH
  178.       IF(JSGN.EQ.1) THEN
  179.          K0=1
  180.       ELSE
  181.          K0=NX
  182.       ENDIF
  183.       DO 1 I=0,IMAX
  184.           NSTART=0
  185.           NLINE=1
  186.           XA(1)=XSTART+I*DX
  187.           XA(2)=XA(1)
  188.           K=K0+JSGN*I
  189.           DO 2 J=0,JMAX-1
  190.               YA(1)=YMIN+J*DY
  191.               YA(2)=YA(1)+DY
  192.               ZA(1)=A(K,J+1)
  193.               ZA(2)=A(K,J+2)
  194.               K1=K+JSGN
  195.               IF(K1.GE.1.AND.K1.LE.NX) THEN
  196.                   ZZ=A(K1,J+1)
  197.               ELSE
  198.                   K1=K-JSGN
  199.                   ZZ=2.*ZA(1)-A(K1,J+1)
  200.               ENDIF
  201.           T=Q3-Q2*(ZZ-ZA(1))-Q1*(ZA(2)-ZA(1))
  202.           T=JSGN*T
  203.           IF(T.GT.0.) THEN
  204.               N=1+IFLAG
  205.           ELSE
  206.               N=4
  207.           ENDIF
  208.           CALL PROJ(XA,YA,ZA,N,NX)
  209.     2     CONTINUE
  210.     1 CONTINUE
  211. C
  212. C Draw curves stepped in y
  213. C
  214.       DY=ISGN*DY
  215.       DX=JSGN*DX
  216.       Q1=JSGN*Q1
  217.       Q2=ISGN*Q2
  218.       Q3=ISGN*JSGN*Q3
  219.       IF(ISGN.EQ.1) THEN
  220.          K0=1
  221.       ELSE
  222.          K0=NY
  223.       ENDIF
  224.       DO 3 J=0,JMAX
  225.           NSTART=0
  226.           NLINE=1
  227.           YA(1)=YSTART+J*DY
  228.           YA(2)=YA(1)
  229.           K=K0+ISGN*J
  230.           DO 4 I=0,IMAX-1
  231.               XA(1)=XMIN+I*DX
  232.               XA(2)=XA(1)+DX
  233.               ZA(1)=A(I+1,K)
  234.               ZA(2)=A(I+2,K)
  235.               K1=K+ISGN
  236.               IF(K1.GE.1.AND.K1.LE.NY) THEN
  237.                   ZZ=A(I+1,K1)
  238.               ELSE
  239.                   K1=K-ISGN
  240.                   ZZ=2.*ZA(1)-A(I+1,K1)
  241.               ENDIF
  242.               T=Q3-Q1*(ZZ-ZA(1))-Q2*(ZA(2)-ZA(1))
  243.               T=ISGN*T
  244.               IF(T.GT.0.) THEN
  245.                   N=1+IFLAG
  246.               ELSE
  247.                   N=4
  248.               ENDIF
  249.               CALL PROJ(XA,YA,ZA,N,NY)
  250.     4     CONTINUE
  251.     3 CONTINUE
  252.       CALL PGEBUF
  253.       RETURN
  254.       END
  255.  
  256.       SUBROUTINE PROJ(X,Y,Z,N,NCURVE)
  257.       REAL X(2),Y(2),Z(2)
  258.       INTEGER N, NCURVE
  259. C
  260.       REAL XP(100),YP(100)
  261.       REAL CTH,STH,CPH,SPH
  262.       INTEGER NLINE,NSTART,IFLAG,M
  263.       COMMON/PGSRF1/CTH,STH,CPH,SPH,NLINE,NSTART,IFLAG
  264.       SAVE M
  265.       IF(NLINE.EQ.1) THEN
  266.           XP(1)=Y(1)*CPH-X(1)*SPH
  267.           YP(1)=Z(1)*STH-(X(1)*CPH+Y(1)*SPH)*CTH
  268.           XP(2)=Y(2)*CPH-X(2)*SPH
  269.           YP(2)=Z(2)*STH-(X(2)*CPH+Y(2)*SPH)*CTH
  270.           IF(N.NE.M) THEN
  271.               M=N
  272.               IF(IFLAG.EQ.0) THEN
  273.                   CALL PGSLS(N)
  274.               ELSE
  275.                   CALL PGSCI(N)
  276.               ENDIF
  277.           ENDIF
  278.           NLINE=2
  279.       ELSE
  280.           IF(N.NE.M) THEN
  281.               CALL PGLINE(NLINE,XP,YP)
  282.               NSTART=NSTART+NLINE-1
  283.               M=N
  284.               IF(IFLAG.EQ.0) THEN
  285.                   CALL PGSLS(N)
  286.               ELSE
  287.                   CALL PGSCI(N)
  288.               ENDIF
  289.               XP(1)=XP(NLINE)
  290.               YP(1)=YP(NLINE)
  291.               NLINE=1
  292.           ENDIF
  293.           NLINE=NLINE+1
  294.           XP(NLINE)=Y(2)*CPH-X(2)*SPH
  295.           YP(NLINE)=Z(2)*STH-(X(2)*CPH+Y(2)*SPH)*CTH
  296.       ENDIF
  297.       IF(NLINE+NSTART.GE.NCURVE) CALL PGLINE(NLINE,XP,YP)
  298.       RETURN
  299.       END
  300.  
  301.       SUBROUTINE WINDOW(XMIN,XMAX,YMIN,YMAX)
  302.       REAL XMIN, XMAX, YMIN, YMAX
  303. C
  304. C This subroutine sets up the standard pgwind, but with a
  305. C cream background. It can be ported to any program.
  306. C
  307.       CALL PGPAGE
  308.       CALL PGSCR(0,216/255.,216/255.,191/255.)
  309.       CALL PGERAS
  310.       CALL PGSWIN(XMIN,XMAX,YMIN,YMAX)
  311.       CALL PGSCI(1)
  312.       RETURN
  313.       END
  314.