c************************************************************ c subroutine file << sub2a.f >> * c * c update ; 01/09/12 /home/g3/ogino/VRML/jst2k1 * c************************************************************ c c --------------------- c subroutine initvrml c --------------------- subroutine initvrml c initialization of VRML write(10,10) 10 format("#VRML V2.0 utf8") return end c --------------------- c subroutine point3d c --------------------- subroutine point3d(nx,ny,nz,x0,y0,z0,xl,yl,zl,ico,u) c c nx: number in x-direction c ny: number in y-direction c nz: number in z-direction c x0: x-value at the left and bottom point c y0: y-value at the left and bottom point c z0: z-value at the left and bottom point c xl: length in x-direction c yl: length in y-direction c zl: length in z-direction c ico: color code c u(nx*ny): input data [0.0,1.0] c dimension u(1) c hx=xl/float(nx-1) hy=yl/float(ny-1) hz=zl/float(nz-1) c call colorb call colorpb(ico,nx,ny,nz,u) call colore c call pointb call pointp(nx,ny,nz,x0,y0,z0,hx,hy,hz) call pointe c return end c --------------------- c subroutine colorb c --------------------- subroutine colorb c write(10,*) " Shape {" write(10,*) " geometry PointSet {" write(10,*) " color Color {" write(10,*) " color [" c return end c --------------------- c subroutine colore c --------------------- subroutine colore write(10,*) " ]" write(10,*) " }" c return end c --------------------- c subroutine pointb c --------------------- subroutine pointb write(10,*) " coord Coordinate {" write(10,*) " point [" return end subroutine pointe write(10,*) " ]" write(10,*) " }" write(10,*) " }" write(10,*) " }" c return end c --------------------- c subroutine colorpb c --------------------- subroutine colorpb(ico,nx,ny,nz,u) dimension u(1) c n1=nx n2=n1*ny n3=n2*nz last=nx*ny*nz c do 10 k=1,nz do 10 j=1,ny do 10 i=1,nx x=hx*float(i-1) y=hy*float(j-1) z=hz*float(k-1) ii=i+n1*(j-1)+n2*(k-1) c x1=u(ii) c a0=1.20 a0=1.00 if(ico.eq.2) go to 32 if(ico.eq.3) go to 33 if(ico.eq.4) go to 34 if(ico.eq.5) go to 35 if(ico.eq.6) go to 36 if(ico.eq.7) go to 37 r=a0-(a0-0.5)*x1 g=a0-(a0-0.5)*x1 b=a0-(a0-0.5)*x1 go to 40 32 continue r=a0-a0*x1 g=1.0 b=a0-a0*x1 go to 40 33 continue r=1.0 g=a0-a0*x1 b=a0-a0*x1 go to 40 34 continue a0=2.0 y1=x1*(a0+4.0) r=a0-y1 if(y1.ge.a0) r=y1-(a0+1.0) g=a0+3.0-y1 b=a0+1.0-y1 if(y1.ge.a0+2.0) b=y1-(a0+3.0) go to 40 35 continue a0=2.0 y1=x1*(a0+4.0) b=a0-y1 if(y1.ge.a0) b=y1-(a0+1.0) r=a0+3.0-y1 g=a0+1.0-y1 if(y1.ge.a0+2.0) g=y1-(a0+3.0) go to 40 36 continue a0=2.0 y1=x1*(a0+4.0) r=a0-y1 if(y1.ge.a0) r=y1-(a0+1.0) b=a0+3.0-y1 g=a0+1.0-y1 if(y1.ge.a0+2.0) g=y1-(a0+3.0) go to 40 37 continue x1=10.0*x1 b=x1 if(x1.ge.3.0) b=4.0-x1 if(x1.ge.6.0) b=0.25*(x1-6.0) g=x1-2.0 if(x1.ge.5.0) g=6.0-x1 if(x1.ge.6.0) g=0.25*(x1-6.0) r=x1 if(x1.ge.1.0) r=2.0-x1 if(x1.ge.4.0) r=x1-4.0 go to 40 c 40 continue r=amax1(r,0.0) g=amax1(g,0.0) b=amax1(b,0.0) r=amin1(r,1.0) g=amin1(g,1.0) b=amin1(b,1.0) c if(ii.lt.last) then write(10,20) r,g,b 20 format(1h ,f3.1,1x,f3.1,1x,f3.1,1h,) else write(10,22) r,g,b 22 format(1h ,f3.1,1x,f3.1,1x,f3.1) end if 10 continue c return end c --------------------- c subroutine pointp c --------------------- subroutine pointp(nx,ny,nz,x0,y0,z0,hx,hy,hz) c n1=nx n2=n1*ny n3=n2*nz last=nx*ny*nz c do 10 k=1,nz z=z0+hz*float(k-1) do 10 j=1,ny y=y0+hy*float(j-1) do 10 i=1,nx x=x0+hx*float(i-1) ii=i+n1*(j-1)+n2*(k-1) if(ii.lt.last) then write(10,20) x,z,y 20 format(1h ,f5.2,1x,f5.2,1x,f5.2,1h,) else write(10,22) x,z,y 22 format(1h ,f5.2,1x,f5.2,1x,f5.2) end if 10 continue return end c --------------------- c subroutine line3d c --------------------- subroutine line3d(lasl,r,g,b,po) c c lasl: number of points on line c r: red color [0.0,1.0] c g: green color [0.0,1.0] c b: blue color [0.0,1.0] c po(3*lasl): x,y,z-values in position c dimension po(1) c call linebe(r,g,b) c do 40 ii=1,lasl lin=ii li1=3*(lin-1)+1 li2=3*(lin-1)+2 li3=3*lin x=po(li1) y=po(li2) z=po(li3) call linep(x,y,z,ii,lasl) 40 continue c call lineen(lasl) return end c --------------------- c subroutine symbvb c --------------------- subroutine symbvb write(10,10) 10 format("Viewpoint { position 10 4 8 }") return end c --------------------- c subroutine symblv c --------------------- subroutine symblv(x,y,z,h,chr,n) c c x: x-position of start point of character c y: y-position of start point of character c z: z-position of start point of character c h: height of character c chr: character c n: number of character c character chr*80,ica*80,ich(80)*1 equivalence (ica,ich) ica=chr write(10,10) c 10 format("Viewpoint { position 10 4 8 }"/"Transform {") 10 format("Transform {") write(10,12) x,y,z 12 format("translation ",3(1x,f5.1)) write(10,14) 14 format("children ["/"Shape {"/"geometry Text {") c write(10,16) (ich(i),i=1,n) c 16 format("string ",'"',80a1,'"') write(10,*) "string ",'"',(ich(i),i=1,n),'"' write(10,18) 18 format("fontStyle FontStyle {") write(10,20) h 20 format("size ",f5.1) write(10,22) 22 format('family "SERIF"'/'style "ROMAN"') write(10,24) 24 format("}"/"}"/"}"/"]"/"}") return end c --------------------- c subroutine symblvc c --------------------- subroutine symblvc(x,y,z,h,r,g,b,chr,n) c c x: x-position of start point of character c y: y-position of start point of character c z: z-position of start point of character c h: height of character c r: red [0.0,1.0] c g: green [0.0,1.0] c b: blue [0.0,1.0] c chr: character c n: number of character c character chr*80,ica*80,ich(80)*1 equivalence (ica,ich) ica=chr write(10,10) c 10 format("Viewpoint { position 10 4 8 }"/"Transform {") 10 format("Transform {") write(10,12) x,y,z 12 format("translation ",3(1x,f5.1)) write(10,14) r,g,b 14 format("children ["/"Shape {"/"appearance Appearance {" 1 /"material Material {"/"diffuseColor ",3(1x,f3.1)/ 2 "}"/"}"/"geometry Text {") c write(10,16) (ich(i),i=1,n) c 16 format("string ",'"',80a1,'"') write(10,*) "string ",'"',(ich(i),i=1,n),'"' write(10,18) 18 format("fontStyle FontStyle {") write(10,20) h 20 format("size ",f5.1) write(10,22) 22 format('family "SERIF"'/'style "ROMAN"') write(10,24) 24 format("}"/"}"/"}"/"]"/"}") return end c --------------------- c subroutine image1 c --------------------- subroutine image1(nx,ny,xb,yb,xl,yl,ico,icc,zcc,u,p) c c nx: number in x-direction c nx: number in x-direction c x0: x-value at the left and bottom point c y0: y-value at the left and bottom point c xl: length in x-direction c yl: length in y-direction c ico: color code c icc: selection of plane, i c : 1 for xy plane, s2 for xz plane, 3 for yz plane c zcc: position of selected plane c u(nx*ny): input data c p: work array c dimension p(1),u(1) ny0=ny-1 nx2=nx*2 c do 100 ii=1,ny0 c call trans2(nx,ny,ii,u,p) call facebe c call colorp(nx2,p,ico) c call coordbe c c call coordp(nx2,p) call coordp(ii,nx,ny,xb,yb,xl,yl,icc,zcc,p) c call coorden(nx) call faceen c 100 continue c return end c --------------------- c subroutine image2 c --------------------- subroutine image2(nx,ny,xb,yb,xl,yl,ico,icc,zcc,v,u,p) c c nx: number in x-direction c nx: number in x-direction c x0: x-value at the left and bottom point c y0: y-value at the left and bottom point c xl: length in x-direction c yl: length in y-direction c ico: color code c icc: selection of plane, i c : 1 for xy plane, s2 for xz plane, 3 for yz plane c zcc: position of selected plane c v: (nx-1)*(ny-1) additional input data, unnecessary c u(nx*ny): input data c p: work array c dimension p(1),u(1),v(1) ny0=ny-1 nx2=nx*2 c do 100 ii=1,ny0 c call trans2(nx,ny,ii,u,p) call facebe c call colorp(nx2,p,ico) c call coordbe c c call coordp(nx2,p) call coordp(ii,nx,ny,xb,yb,xl,yl,icc,zcc,p) c call coorden(nx) call faceen c 100 continue c return end c --------------------- c subroutine trans1 c --------------------- subroutine trans1(nx,ny,v,u) dimension v(1),u(1) c nx0=nx-1 ny0=ny-1 nxa=nx0-1 nya=ny0-1 c do 10 j=1,nya do 10 i=1,nxa i1=i+1+nx*j j1=i+nx0*(j-1) u(i1)=0.25*(v(j1)+v(j1+1)+v(j1+nx0)+v(j1+nx)) 10 continue c do 20 i=1,nxa i1=i+1 i2=i1+nx*ny0 j1=i+1 j2=i+1+nx0*nya u(i1)=0.75*(v(j1-1)+v(j1))-0.25*(v(j1-1+nx0)+v(j1+nx0)) u(i2)=0.75*(v(j2-1)+v(j2))-0.25*(v(j2-1-nx0)+v(j2-nx0)) 20 continue c do 30 j=1,nya i1=1+nx*j i2=i1+nx0 j1=1+nx0*j j2=j1+nxa u(i1)=0.75*(v(j1-nx0)+v(j1))-0.25*(v(j1-nx0+1)+v(j1+1)) u(i2)=0.75*(v(j2-nx0)+v(j2))-0.25*(v(j2-nx0-1)+v(j2-1)) 30 continue c i1=1 i2=i1+nx0 i3=1+nx*ny0 i4=i3+nx0 j1=1 j2=j1+nxa j3=1+nx0*nya j4=j3+nxa u(i1)=0.5*(3.0*v(j1)-v(j1+nx)) u(i2)=0.5*(3.0*v(j2)-v(j2+nxa)) u(i3)=0.5*(3.0*v(j3)-v(j3-nxa)) u(i4)=0.5*(3.0*v(j4)-v(j4-nx)) c c write(6,*) i1,i2,i3,i4 c write(6,*) j1,j2,j3,j4 c write(6,*) nx,nx0,nxa,ny,ny0,nya c return end c --------------------- c subroutine trans2 c --------------------- subroutine trans2(nx,ny,ii,u,p) dimension p(1),u(1) do 10 i=1,nx j1=i+nx*(ii-1) j2=j1+nx i1=1+2*(i-1) i2=i1+1 p(i1)=u(j1) p(i2)=u(j2) 10 continue return end c --------------------- c subroutine colorp c --------------------- subroutine colorp(nx,p,ico) dimension p(1) r=1.0 g=1.0 b=1.0 do 100 i=1,nx x1=p(i) c a0=1.20 a0=1.00 if(ico.eq.2) go to 32 if(ico.eq.3) go to 33 if(ico.eq.4) go to 34 if(ico.eq.5) go to 35 if(ico.eq.6) go to 36 if(ico.eq.7) go to 37 r=a0-(a0-0.5)*x1 g=a0-(a0-0.5)*x1 b=a0-(a0-0.5)*x1 go to 40 32 continue r=a0-a0*x1 g=1.0 b=a0-a0*x1 go to 40 33 continue r=1.0 g=a0-a0*x1 b=a0-a0*x1 go to 40 34 continue a0=2.0 y1=x1*(a0+4.0) r=a0-y1 if(y1.ge.a0) r=y1-(a0+1.0) g=a0+3.0-y1 b=a0+1.0-y1 if(y1.ge.a0+2.0) b=y1-(a0+3.0) go to 40 35 continue a0=2.0 y1=x1*(a0+4.0) b=a0-y1 if(y1.ge.a0) b=y1-(a0+1.0) r=a0+3.0-y1 g=a0+1.0-y1 if(y1.ge.a0+2.0) g=y1-(a0+3.0) go to 40 36 continue a0=2.0 y1=x1*(a0+4.0) r=a0-y1 if(y1.ge.a0) r=y1-(a0+1.0) b=a0+3.0-y1 g=a0+1.0-y1 if(y1.ge.a0+2.0) g=y1-(a0+3.0) go to 40 37 continue x1=10.0*x1 b=x1 if(x1.ge.3.0) b=4.0-x1 if(x1.ge.6.0) b=0.25*(x1-6.0) g=x1-2.0 if(x1.ge.5.0) g=6.0-x1 if(x1.ge.6.0) g=0.25*(x1-6.0) r=x1 if(x1.ge.1.0) r=2.0-x1 if(x1.ge.4.0) r=x1-4.0 go to 40 c 40 continue r=amax1(r,0.0) g=amax1(g,0.0) b=amax1(b,0.0) r=amin1(r,1.0) g=amin1(g,1.0) b=amin1(b,1.0) c if(i.ne.nx) write(10,20) r,g,b if(i.eq.nx) write(10,30) r,g,b 20 format(f4.2,1x,f4.2,1x,f4.2,",") 30 format(f4.2,1x,f4.2,1x,f4.2) 100 continue return end c --------------------- c subroutine coordp c --------------------- subroutine coordp(ii,nx,ny,xb,yb,xl,yl,icc,zcc,p) dimension p(1) x=1.0 y=1.0 z=1.0 c xb=0.0 c yb=0.0 c dx=1.0 c dy=1.0 nxy=nx*ny dx=xl/float(nx-1) dy=yl/float(ny-1) c do 100 j1=1,2 do 100 i=1,nx do 100 j1=1,2 j=j1-1+ii ixy=i+nx*(j-1) c if (icc.eq.2) go to 52 if (icc.eq.3) go to 53 c x=xb+dx*float(i-1) y=yb+dy*float(j-1) z=zcc go to 60 52 continue y=xb+dx*float(i-1) z=yb+dy*float(j-1) x=zcc go to 60 53 continue x=xb+dx*float(i-1) z=yb+dy*float(j-1) y=zcc 60 continue c if(i.ne.nx) write(10,20) x,y,z if(i.eq.nx) write(10,30) x,y,z 20 format(f4.2,1x,f4.2,1x,f4.2,",") 30 format(f4.2,1x,f4.2,1x,f4.2) 100 continue return end c --------------------- c subroutine facebe c --------------------- subroutine facebe write(10,21) 21 format("Shape {") write(10,22) 22 format("geometry IndexedFaceSet {") write(10,23) 23 format("color Color {") write(10,24) 24 format("color [") return end c --------------------- c subroutine coordbe c --------------------- subroutine coordbe write(10,21) 21 format("]") write(10,22) 22 format("}") write(10,23) 23 format("coord Coordinate {") write(10,24) 24 format("point [") return end c --------------------- c subroutine coorden c --------------------- subroutine coorden(nx) dimension ia(2000) write(10,21) 21 format("]") write(10,22) 22 format("}") write(10,23) 23 format("colorPerVertex TRUE") write(10,24) 24 format("coordIndex [") c do 10 i=1,nx i1=1+5*(i-1) i2=2+5*(i-1) i3=3+5*(i-1) i4=4+5*(i-1) i5=5+5*(i-1) ia(i1)=0+2*(i-1) ia(i2)=1+2*(i-1) ia(i3)=3+2*(i-1) ia(i4)=2+2*(i-1) ia(i5)=-1 10 continue c na2=5*(nx-1)-1 write(10,102) (ia(i),i=1,na2) 102 format(4(4i2,i3),2i2,3i3,220i3,2i3,2i4,i3, 1 449(4i4,i3),2i4,2i5,i3,3000(4i5,i3)) c return end c --------------------- c subroutine faceen c --------------------- subroutine faceen write(10,21) 21 format("]") write(10,22) 22 format("solid FALSE") write(10,23) write(10,23) 23 format("}") return end c --------------------- c subroutine linep c --------------------- subroutine linep(x,y,z,ii,lasl) c if(ii.lt.lasl) then write(10,10) x,z,y c 10 format(f5.1,1x,f5.1,1x,f5.1,1h,) c 10 format(f5.2,1x,f5.2,1x,f5.2,1h,) 10 format(f6.2,1x,f6.2,1x,f6.2,1h,) else write(10,20) x,z,y c 20 format(f5.1,1x,f5.1,1x,f5.1) c 20 format(f5.2,1x,f5.2,1x,f5.2) 20 format(f6.2,1x,f6.2,1x,f6.2) end if return end c --------------------- c subroutine linep2 c --------------------- subroutine linep2(x,y,z,ii,lasl) c if(ii.lt.lasl) then write(10,10) x,z,y 10 format(f5.1,1x,f5.1,1x,f5.1,1h,) else write(10,20) x,z,y 20 format(f5.1,1x,f5.1,1x,f5.1) end if return end c --------------------- c subroutine linebe c --------------------- subroutine linebe(r,g,b) c write(10,21) 21 format("Shape {") write(10,22) 22 format("geometry IndexedLineSet {") write(10,23) 23 format("color Color {") write(10,20) r,g,b 20 format("color [ ",f3.1,1x,f3.1,1x,f3.1," ]") write(10,25) 25 format("}") write(10,26) 26 format("coord Coordinate {") write(10,27) 27 format("point [") return end c --------------------- c subroutine lineen c --------------------- subroutine lineen(lasl) integer ic(9000) c do 10 i=1,9000 ic(i)=i-1 10 continue write(10,21) 21 format("]") write(10,22) 22 format("}") write(10,23) 23 format("colorIndex [ 0 ]") write(10,24) 24 format("coordIndex [") c write(10,20) (ic(i),i=1,lasl) 20 format(10i2,90i3,900i4/(2000i5)) c write(10,25) 25 format("]") write(10,26) 26 format("colorPerVertex FALSE") write(10,27) 27 format("}") write(10,28) 28 format("}") return end c --------------------- c subroutine rect1 c --------------------- subroutine rect1(icc,zcc,r,g,b,x1,y1,x2,y2) lasl=5 call linebe(r,g,b) c if (icc.eq.2) go to 52 if (icc.eq.3) go to 53 c z=zcc call linep(x1,y1,z,1,lasl) call linep(x2,y1,z,2,lasl) call linep(x2,y2,z,3,lasl) call linep(x1,y2,z,4,lasl) call linep(x1,y1,z,5,lasl) go to 60 c 52 continue x=zcc call linep(x,x1,y1,1,lasl) call linep(x,x2,y1,2,lasl) call linep(x,x2,y2,3,lasl) call linep(x,x1,y2,4,lasl) call linep(x,x1,y1,5,lasl) go to 60 c 53 continue y=zcc call linep(x1,y,y1,1,lasl) call linep(x2,y,y1,2,lasl) call linep(x2,y,y2,3,lasl) call linep(x1,y,y2,4,lasl) call linep(x1,y,y1,5,lasl) 60 continue call lineen(lasl) c return end c zcsub1.f subroutine program c --------------------- c subroutine zcoordp c --------------------- subroutine zcoordp(nx,p) dimension p(1) c do 100 i=1,nx i1=1+3*(i-1) i2=i1+1 i3=i2+1 x=p(i1) y=p(i2) z=p(i3) if(i.ne.nx) write(10,20) x,y,z if(i.eq.nx) write(10,30) x,y,z c 20 format(f3.1,1x,f3.1,1x,f3.1,",") c 30 format(f3.1,1x,f3.1,1x,f3.1) 20 format(f5.3,1x,f5.3,1x,f5.3,",") 30 format(f5.3,1x,f5.3,1x,f5.3) 100 continue return end c --------------------- c subroutine zcoorden c --------------------- subroutine zcoorden(nx,ib) dimension ib(1) write(10,21) 21 format("]") write(10,22) 22 format("}") write(10,23) 23 format("colorPerVertex TRUE") write(10,24) 24 format("coordIndex [") c c do 10 i=1,nx c ib(i)=i-1 c 10 continue c write(10,102) (ib(i),i=1,nx) 102 format(10i5) c return end c c zsub12.f cube subroutine program c --------------------- c subroutine csub2 c --------------------- subroutine csub2(ic,ib,aa,p,pt,ct) parameter(nol=12,nop=3*nol) c dimension pp(nop),p(8),nf(nol,6),ii(6),jj(6) dimension pp(nop),p(8),nf(nol,6),ii(10),jj(10) dimension ic(1),ib(1),aa(1),pt(1),ct(1) c c ico=ic(1) ua=aa(1) x1=aa(2) x2=aa(3) y1=aa(4) y2=aa(5) z1=aa(6) z2=aa(7) c do 10 j1=1,6 do 10 i1=1,nol nf(i1,j1)=0 10 continue c do 12 i1=1,nop pp(i1)=-1.0 12 continue c u1=(ua-p(1))*(ua-p(2)) if(u1.lt.0.0) then nf(1,1)=1 nf(1,3)=1 x=(x1*(p(2)-ua)+x2*(ua-p(1)))/(p(2)-p(1)) y=y1 z=z1 i1=1 i2=i1+1 i3=i2+1 pp(i1)=x pp(i2)=y pp(i3)=z end if c u1=(ua-p(2))*(ua-p(3)) if(u1.lt.0.0) then nf(2,1)=1 nf(2,4)=1 x=x2 y=(y1*(p(3)-ua)+y2*(ua-p(2)))/(p(3)-p(2)) z=z1 i1=1+3*1 i2=i1+1 i3=i2+1 pp(i1)=x pp(i2)=y pp(i3)=z end if ccc c do 242 j1=1,12 c i1=1+3*(j1-1) c i2=3*j1 c write(6,22) j1,(pp(ia),ia=i1,i2) c 242 continue c u1=(ua-p(3))*(ua-p(4)) if(u1.lt.0.0) then nf(3,1)=1 nf(3,5)=1 x=(x2*(p(4)-ua)+x1*(ua-p(3)))/(p(4)-p(3)) y=y2 z=z1 i1=1+3*2 i2=i1+1 i3=i2+1 pp(i1)=x pp(i2)=y pp(i3)=z end if c u1=(ua-p(4))*(ua-p(1)) if(u1.lt.0.0) then nf(4,1)=1 nf(4,6)=1 x=x1 y=(y2*(p(1)-ua)+y1*(ua-p(4)))/(p(1)-p(4)) z=z1 i1=1+3*3 i2=i1+1 i3=i2+1 pp(i1)=x pp(i2)=y pp(i3)=z end if c c u1=(ua-p(5))*(ua-p(6)) if(u1.lt.0.0) then nf(5,2)=1 nf(5,3)=1 x=(x1*(p(6)-ua)+x2*(ua-p(5)))/(p(6)-p(5)) y=y1 z=z2 i1=1+3*4 i2=i1+1 i3=i2+1 pp(i1)=x pp(i2)=y pp(i3)=z end if c u1=(ua-p(6))*(ua-p(7)) if(u1.lt.0.0) then nf(6,2)=1 nf(6,4)=1 x=x2 y=(y1*(p(7)-ua)+y2*(ua-p(6)))/(p(7)-p(6)) z=z2 i1=1+3*5 i2=i1+1 i3=i2+1 pp(i1)=x pp(i2)=y pp(i3)=z end if c u1=(ua-p(7))*(ua-p(8)) if(u1.lt.0.0) then nf(7,2)=1 nf(7,5)=1 x=(x2*(p(8)-ua)+x1*(ua-p(7)))/(p(8)-p(7)) y=y2 z=z2 i1=1+3*6 i2=i1+1 i3=i2+1 pp(i1)=x pp(i2)=y pp(i3)=z end if c u1=(ua-p(8))*(ua-p(5)) if(u1.lt.0.0) then nf(8,2)=1 nf(8,6)=1 x=x1 y=(y2*(p(5)-ua)+y1*(ua-p(8)))/(p(5)-p(8)) z=z2 i1=1+3*7 i2=i1+1 i3=i2+1 pp(i1)=x pp(i2)=y pp(i3)=z end if c c u1=(ua-p(1))*(ua-p(5)) if(u1.lt.0.0) then nf(9,3)=1 nf(9,6)=1 x=x1 y=y1 z=(z1*(p(5)-ua)+z2*(ua-p(1)))/(p(5)-p(1)) i1=1+3*8 i2=i1+1 i3=i2+1 pp(i1)=x pp(i2)=y pp(i3)=z end if c u1=(ua-p(2))*(ua-p(6)) if(u1.lt.0.0) then nf(10,3)=1 nf(10,4)=1 x=x2 y=y1 z=(z1*(p(6)-ua)+z2*(ua-p(2)))/(p(6)-p(2)) i1=1+3*9 i2=i1+1 i3=i2+1 pp(i1)=x pp(i2)=y pp(i3)=z end if c u1=(ua-p(3))*(ua-p(7)) if(u1.lt.0.0) then nf(11,4)=1 nf(11,5)=1 x=x2 y=y2 z=(z1*(p(7)-ua)+z2*(ua-p(3)))/(p(7)-p(3)) i1=1+3*10 i2=i1+1 i3=i2+1 pp(i1)=x pp(i2)=y pp(i3)=z end if c u1=(ua-p(4))*(ua-p(8)) if(u1.lt.0.0) then nf(12,5)=1 nf(12,6)=1 x=x1 y=y2 z=(z1*(p(8)-ua)+z2*(ua-p(4)))/(p(8)-p(4)) i1=1+3*11 i2=i1+1 i3=i2+1 pp(i1)=x pp(i2)=y pp(i3)=z end if c c c write(6,20) (p(i1),i1=1,8) 20 format(1h ,4f10.4) do 24 j1=1,12 i1=1+3*(j1-1) i2=3*j1 c write(6,22) j1,(pp(ia),ia=i1,i2) 22 format(1h ,i5,3f10.4) 24 continue do 28 j1=1,12 c write(6,26) j1,(nf(j1,ia),ia=1,6) 26 format(1h ,i5,6i5) 28 continue c c call csub1(ic,ib,p,pp,nf,pt,ct) c c 100 continue return end c c c c cube subroutine program c --------------------- c subroutine csub1 c --------------------- subroutine csub1(ic,ib,p,pp,nf,pt,ct) parameter(nol=12,nop=3*nol) c dimension p(8),pp(nop),nf(nol,6),ii(6),jj(6),ic(1) dimension p(8),pp(nop),nf(nol,6),ii(10),jj(10),ic(1) dimension ib(1),pa(nop),pt(1),ct(1) c ico=ic(1) iin=0 do 10 i=1,4 if(nf(i,1).eq.1) then iin=iin+1 ii(iin)=i jj(iin)=1 end if 10 continue c do 12 i=5,8 if(nf(i,2).eq.1) then iin=iin+1 ii(iin)=i jj(iin)=2 end if 12 continue cc cc if(iin.eq.0) go to 100 if(iin.eq.2.and.ii(2).le.4) go to 200 if(iin.eq.2.and.ii(2).ge.5) go to 300 if(iin.eq.4) go to 400 go to 990 c 100 continue do 110 i=9,12 j=i-6 if(nf(i,j).eq.1) then iin=iin+1 ii(iin)=i jj(iin)=j end if 110 continue go to 1000 c c 200 continue i1=ii(1) i2=ii(2) j1=jj(1) j2=jj(2) c do 210 j=3,6 if(nf(i2,j).eq.1) then j3=j end if 210 continue do 212 i=9,12 if(i.eq.i2) go to 212 if(nf(i,j3).eq.1) then iin=iin+1 i3=i ii(iin)=i3 jj(iin)=j3 end if 212 continue c i4=0 do 220 j=3,6 if(j.eq.j3) go to 220 if(nf(i3,j).eq.1) then j4=j end if 220 continue do 222 i=9,12 if(i.eq.i2) go to 222 if(i.eq.i3) go to 222 if(nf(i,j4).eq.1) then iin=iin+1 i4=i ii(iin)=i4 jj(iin)=j4 end if 222 continue if(i4.eq.0) go to 1000 c do 230 j=3,6 if(j.eq.j3) go to 230 if(j.eq.j4) go to 230 if(nf(i4,j).eq.1) then j5=j end if 230 continue do 232 i=9,12 if(i.eq.i2) go to 232 if(i.eq.i3) go to 232 if(i.eq.i4) go to 232 if(nf(i,j5).eq.1) then iin=iin+1 i5=i ii(iin)=i5 jj(iin)=j5 end if 232 continue go to 1000 c c 300 continue i1=ii(1) i2=ii(2) j1=jj(1) j2=jj(2) c do 310 j=3,6 if(nf(i2,j).eq.1) then j3=j end if 310 continue do 312 i=9,12 if(i.eq.i2) go to 312 if(nf(i,j3).eq.1) then iin=iin+1 i3=i ii(iin)=i3 jj(iin)=j3 end if 312 continue c i4=0 do 320 j=3,6 if(j.eq.j3) go to 320 if(nf(i3,j).eq.1) then j4=j end if 320 continue do 322 i=9,12 if(i.eq.i2) go to 322 if(i.eq.i3) go to 322 if(nf(i,j4).eq.1) then iin=iin+1 i4=i ii(iin)=i4 jj(iin)=j4 end if 322 continue if(i4.eq.0) go to 1000 c do 330 j=3,6 if(j.eq.j3) go to 330 if(j.eq.j4) go to 330 if(nf(i4,j).eq.1) then j5=j end if 330 continue do 332 i=9,12 if(i.eq.i2) go to 332 if(i.eq.i3) go to 332 if(i.eq.i4) go to 332 if(nf(i,j5).eq.1) then iin=iin+1 i5=i ii(iin)=i5 jj(iin)=j5 end if 332 continue go to 1000 c c 400 continue iin=2 i1=ii(1) i2=ii(2) j1=jj(1) j2=jj(2) c do 410 j=3,6 if(nf(i2,j).eq.1) then j3=j end if 410 continue do 412 i=5,12 if(i.eq.i2) go to 412 if(nf(i,j3).eq.1) then iin=iin+1 i3=i ii(iin)=i3 jj(iin)=j3 end if 412 continue c i4=0 do 420 j=2,6 if(j.eq.j3) go to 420 if(nf(i3,j).eq.1) then j4=j end if 420 continue do 422 i=5,12 if(i.eq.i2) go to 422 if(i.eq.i3) go to 422 if(nf(i,j4).eq.1) then iin=iin+1 i4=i ii(iin)=i4 jj(iin)=j4 end if 422 continue if(i4.eq.0) go to 1000 c i5=0 do 430 j=2,6 if(j.eq.j3) go to 430 if(j.eq.j4) go to 430 if(nf(i4,j).eq.1) then j5=j end if 430 continue do 432 i=5,12 if(i.eq.i2) go to 432 if(i.eq.i3) go to 432 if(i.eq.i4) go to 432 if(nf(i,j5).eq.1) then iin=iin+1 i5=i ii(iin)=i5 jj(iin)=j5 end if 432 continue if(i5.eq.0) go to 1000 c do 440 j=2,6 if(j.eq.j3) go to 440 if(j.eq.j4) go to 440 if(j.eq.j5) go to 440 if(nf(i5,j).eq.1) then j6=j end if 440 continue do 442 i=5,12 if(i.eq.i2) go to 442 if(i.eq.i3) go to 442 if(i.eq.i4) go to 442 if(i.eq.i5) go to 442 if(nf(i,j6).eq.1) then iin=iin+1 i6=i ii(iin)=i6 jj(iin)=j6 end if 442 continue c c 1000 continue do 1010 j=1,iin ii1=ii(j) jj1=jj(j) i1=1+3*(ii1-1) i2=i1+1 i3=i2+1 j1=1+3*(j-1) j2=j1+1 j3=j2+1 c ic(5)=ic(5)+1 ic(6)=ic(6)+1 k=ic(5) ka=ic(6) k1=1+3*(k-1) k2=k1+1 k3=k2+1 ct(k)=p(j) ib(ka)=k-1 c x=pp(i1) y=pp(i2) z=pp(i3) pt(k1)=x pt(k2)=y pt(k3)=z c write(6,1020) j,ii1,jj1,x,y,z 1020 format(1h ,3i4,3f10.4) 1010 continue c if(iin.eq.0) go to 990 ic(6)=ic(6)+1 ka=ic(6) ib(ka)=-1 c c 990 continue return end c cube subroutime program c --------------------- c subroutine csuba c --------------------- subroutine csuba(ic,aa,u) c isosurface plot c c ic(10): c ic(1)=ico,ic(2)=nx,ic(3)=ny,ic(4)=nz c aa(20): c aa(1)=ua,aa(2)=x1,aa(3)=x2,aa(4)=y1,aa(5)=y2 c aa(6)=z1,aa(7)=z2,aa(8)=xl,aa(9)=yl,aa(10)=zl c aa(11)=x0,aa(12)=y0,aa(13)=z0 c u(nx*ny*nz): input data c parameter(ipt=9999,ipp=ipt*3) parameter(nol=12,nop=3*nol) c dimension pp(nop),p(8),nf(nol,6) dimension pp(nop),p(8),nf(nol,6) dimension ic(1),aa(1),u(1) dimension ib(ipt),ct(ipt),pt(ipp) c c c dimension ic(1),aa(1),u(1),pp(nop),nf(nol,6) c dimension p(8) c ico=ic(1) nx=ic(2) ny=ic(3) nz=ic(4) ua=aa(1) x1=aa(2) x2=aa(3) y1=aa(4) y2=aa(5) z1=aa(6) z2=aa(7) xl=aa(8) yl=aa(9) zl=aa(10) x0=aa(11) y0=aa(12) z0=aa(13) c n1=nx n2=n1*ny n3=n2*nz nxs=nx-1 nys=ny-1 nzs=nz-1 dx=xl/float(nxs) dy=yl/float(nys) dz=zl/float(nzs) c c do 200 k=1,nzs z1=z0+dz*float(k-1) z2=z1+dz ic(5)=0 ic(6)=0 c do 100 j=1,nys y1=y0+dy*float(j-1) y2=y1+dy do 100 i=1,nxs x1=x0+dx*float(i-1) x2=x1+dx c i1=i+n1*(j-1)+n2*(k-1) i2=i1+1 i3=i2+nx i4=i1+nx i5=i1+n2 i6=i5+1 i7=i6+nx i8=i5+nx p(1)=u(i1) p(2)=u(i2) p(3)=u(i3) p(4)=u(i4) p(5)=u(i5) p(6)=u(i6) p(7)=u(i7) p(8)=u(i8) c c aa(1)=ua aa(2)=x1 aa(3)=x2 aa(4)=y1 aa(5)=y2 aa(6)=z1 aa(7)=z2 c c c write(6,12) (ic(i1),i1=1,10) c write(6,14) (aa(i1),i1=1,20) c write(6,16) (p(i1),i1=1,8) 12 format(1h ,10i5) 14 format(1h ,5f10.5) 16 format(1h ,4f10.5) c call csub2(ic,ib,aa,p,pt,ct) c 100 continue c iin1=ic(5) iin2=ic(6) c write(6,12) k,iin1,iin2 c if(iin1.eq.0) go to 200 c call facebe c call colorp(iin1,ct,ico) c call coordbe c call zcoordp(iin1,pt) c iin21=iin2-1 call zcoorden(iin21,ib) call faceen c 200 continue return end c c c c************************************************************ c subroutine file << msubrtn2.f >> * c * c update ; 01/03/24 /home/g3/inayoshi/VRML/subrtn4 * c************************************************************ c------------------------------------- subroutine triangm1(nx,ny,u,p) c------------------------------------- dimension p(1),u(1) ny1=ny-1 c call facebe c do 10 j=1,ny1 call colorp2(nx,j,u) 10 continue c call coordbe c do 20 j=1,ny1 call coordp2(nx,j,p) 20 continue c call coorden5(nx,ny1) c call ifaceen c return end c------------------------------------- subroutine triangm(nx,ny,u,p) c------------------------------------- dimension p(1),u(1) ny1=ny-1 c do 10 j=1,ny1 call triangb(nx,j,u,p) 10 continue c return end c------------------------------------- subroutine triangb(nx,ny,u,p) c------------------------------------- dimension p(1),u(1) c call facebe call colorp2(nx,ny,u) call coordbe call coordp2(nx,ny,p) call coorden1(nx) call ifaceen c return end c------------------------------------------- subroutine linemesh(p,n1,rd,r,g,b) c------------------------------------------- dimension p(1) dimension p1(3*2) dimension q(3*12*2) c n11=n1-2 do 20 i=1,n11 do 10 j=1,2 c i1=3*(i-1)+1+3*(j-1) i2=i1+1 i3=i1+2 c---<< chusin ten set >> p1(1+(j-1)*3)=p(i1) p1(2+(j-1)*3)=p(i2) p1(3+(j-1)*3)=p(i3) c xa=abs(p(i1)-p(i1+3)) ya=abs(p(i2)-p(i2+3)) za=abs(p(i3)-p(i3+3)) amx=amax1(xa,ya,za) c x=p(i1) y=p(i2) z=p(i3) c if(amx.eq.xa) then k=1 else if(amx.eq.ya) then k=2 else if(amx.eq.za) then k=3 end if end if end if c c---<< ensyu no point wo motomeru >> c call ensyu(x,y,z,rd,q,j,k) 10 continue c call linemeshb(12,1,q,rd,r,g,b,k,p1) 20 continue c---<< cone >> rr=1.57 if(k.eq.1) then x1=0.0 y1=0.0 z1=-1.0 end if if(k.eq.3) then x1=-1.0 y1=0.0 z1=0.0 end if c rd2=rd*2 h=rd*3 x=p1(4) y=p1(5) z=p1(6) call cone1(rd2,h,r,g,b,x,y,z,x1,y1,z1,rr) c return end c-------------------------------------------------- subroutine linemeshb(nx,ny,q,rd,r,g,b,k,p1) c * 1 color * c-------------------------------------------------- dimension q(1) c call facebe c c----( color ) write(10,10) r,g,b 10 format(f9.2,1x,f9.2,1x,f9.2) c call coordbe1(nx) call coordp4(nx,ny,q,rd,k,p1) call coorden3(nx) call ifaceen c return end c--------------------------------------- subroutine ensyu(x,y,z,rd,q,j,k) c--------------------------------------- dimension q(1) c x1=rd/2*1.73205 y1=rd/2 j2=36*(j-1) c--[ max=xa ] if(k.eq.1) then c (x) do 10 i=1,12 q(1+(i-1)*3+j2)=x 10 continue c (y) q(2+j2)=y+rd q(5+j2)=y+x1 q(8+j2)=y+y1 q(11+j2)=y q(14+j2)=y-y1 q(17+j2)=y-x1 q(20+j2)=y-rd q(23+j2)=y-x1 q(26+j2)=y-y1 q(29+j2)=y q(32+j2)=y+y1 q(35+j2)=y+x1 c (z) q(3+j2)=z q(6+j2)=z+y1 q(9+j2)=z+x1 q(12+j2)=z+rd q(15+j2)=z+x1 q(18+j2)=z+y1 q(21+j2)=z q(24+j2)=z-y1 q(27+j2)=z-x1 q(30+j2)=z-rd q(33+j2)=z-x1 q(36+j2)=z-y1 go to 99 end if c--[ max=ya ] if(k.eq.2) then c (x) q(1+j2)=x q(4+j2)=x+y1 q(7+j2)=x+x1 q(10+j2)=x+rd q(13+j2)=x+x1 q(16+j2)=x+y1 q(19+j2)=x q(22+j2)=x-y1 q(25+j2)=x-x1 q(28+j2)=x-rd q(31+j2)=x-x1 q(34+j2)=x-y1 c (y) do 20 i=1,12 q(2+(i-1)*3+j2)=y 20 continue c (z) q(2+j2)=z+rd q(5+j2)=z+x1 q(8+j2)=z+y1 q(11+j2)=z q(14+j2)=z-y1 q(17+j2)=z-x1 q(20+j2)=z-rd q(23+j2)=z-x1 q(26+j2)=z-y1 q(29+j2)=z q(32+j2)=z+y1 q(35+j2)=z+x1 go to 99 end if c--[ max=za ] if(k.eq.3) then c (x) q(1+j2)=x q(4+j2)=x+y1 q(7+j2)=x+x1 q(10+j2)=x+rd q(13+j2)=x+x1 q(16+j2)=x+y1 q(19+j2)=x q(22+j2)=x-y1 q(25+j2)=x-x1 q(28+j2)=x-rd q(31+j2)=x-x1 q(34+j2)=x-y1 c (y) q(2+j2)=y+rd q(5+j2)=y+x1 q(8+j2)=y+y1 q(11+j2)=y q(14+j2)=y-y1 q(17+j2)=y-x1 q(20+j2)=y-rd q(23+j2)=y-x1 q(26+j2)=y-y1 q(29+j2)=y q(32+j2)=y+y1 q(35+j2)=y+x1 c (z) do 30 i=1,12 q(3*i+j2)=z 30 continue end if c 99 continue return end c-------------------------------------------------- subroutine cylinderm(nx,ny,u,p,rd,br,bg,bb) c * gradation color * c-------------------------------------------------- dimension p(1),u(1) ny1=ny-1 c do 10 j=1,ny1 call cylinderb(nx,j,u,p,rd,br,bg,bb) 10 continue c return end c-------------------------------------------------- subroutine cylinderb(nx,ny,u,p,rd,br,bg,bb) c * gradation color * c-------------------------------------------------- dimension p(1),u(1) c call facebe call colorp3(nx,ny,u,br,bg,bb) call coordbe call coordp3(nx,ny,p,rd) call coorden2(nx) call ifaceen c return end c-------------------------------------------------- subroutine cylinderm1(nx,ny,p,rd,r,g,b) c * 1 color * c-------------------------------------------------- dimension p(1) ny1=ny-1 c do 10 j=1,ny1 call cylinderb1(nx,j,p,rd,r,g,b) 10 continue c return end c-------------------------------------------------- subroutine cylinderb1(nx,ny,p,rd,r,g,b) c * 1 color * c-------------------------------------------------- dimension p(1) c call facebe c----( color ) write(10,10) r,g,b 10 format(f8.2,1x,f8.2,1x,f8.2) c call coordbe1(nx) call coordp3(nx,ny,p,rd) call coorden3(nx) call ifaceen c return end c---------------------------------------------------------------------- subroutine arrow(nx,ny,u,p,rd1,rd2,h,cx,cy,cz,br,bg,bb,cr,cg,cb) c---------------------------------------------------------------------- dimension p(1),u(1) ny1=ny-1 c do 10 j=1,ny1 call cylinderb(nx,j,u,p,rd1,br,bg,bb) 10 continue c x1=0.0 y1=0.0 z1=1.0 rr=1.57 c call cone1(rd2,h,cr,cg,cb,cx,cy,cz,x1,y1,z1,rr) return end c---------------------------------------------------------------------- subroutine arrow2(nx,ny,u,p,rd1,rd2,h,cx,cy,cz,br,bg,bb,cr,cg,cb) c---------------------------------------------------------------------- dimension p(1),u(1) ny1=ny-1 c do 10 j=1,ny1 call cylinderb(nx,j,u,p,rd1,br,bg,bb) 10 continue c jj=ny1 i=nx j=2 j1=j+jj-1 i1=3*(i-1)+1+3*nx*(j1-1) i2=i1+1 i3=i1+2 cx=p(i1)+0.8*rd2 cy=p(i2) cz=p(i3)+rd1 c x1=0.0 y1=0.0 z1=-1.0 rr=1.57 c call cone1(rd2,h,cr,cg,cb,cx,cy,cz,x1,y1,z1,rr) return end c----------------------------------------------------- subroutine cone1(rd,h,r,g,b,x,y,z,x1,y1,z1,rr) c----------------------------------------------------- write(10,20) x,y,z 20 format("Transform {",/"translation ",f7.2,f7.2,f7.2) write(10,21) x1,y1,z1,rr 21 format("rotation ",f7.2,f7.2,f7.2,f7.2) write(10,22) 22 format("children [") write(10,23) 23 format("Shape {") write(10,24) 24 format("appearance Appearance {") write(10,25) 25 format("material Material {") write(10,26) r,g,b 26 format("diffuseColor",f5.2,f5.2,f5.2) write(10,27) 27 format("}",/"}") write(10,28) 28 format("geometry Cone {") write(10,29) rd 29 format("bottomRadius ",f5.2) write(10,30) h 30 format("height ",f5.2) write(10,31) 31 format("}",/"}",/"]",/"}") c return end c------------------------------- subroutine coordbe1(nx) c------------------------------- write(10,21) 21 format("]") write(10,22) 22 format("}") c----( colorindexed ) nx4=nx*4 write(10,20) 20 format("colorIndex [ ") write(10,30) (0,i=1,nx4) 30 format(5000i3) write(10,40) 40 format(" ]") c write(10,23) 23 format("coord Coordinate {") write(10,24) 24 format("point [") return end c-------------------------------------- subroutine colorp1(nx,ny,u) c-------------------------------------- dimension u(1) c do 100 i=1,nx c i1=3*(i-1)+1+3*nx*(ny-1) i2=i1+1 i3=i1+2 r=u(i1) g=u(i2) b=u(i3) c r=amax1(r,0.0) g=amax1(g,0.0) b=amax1(b,0.0) r=amin1(r,1.0) g=amin1(g,1.0) b=amin1(b,1.0) c if(i.eq.nx) then write(10,10) r,g,b 10 format(f6.2,1x,f6.2,1x,f6.2) else write(10,20) r,g,b 20 format(f6.2,1x,f6.2,1x,f6.2,",") end if c 100 continue return end c-------------------------------------- subroutine colorp2(nx,ny,u) c-------------------------------------- dimension u(1) c do 100 i=1,nx do 100 j=1,2 j1=j+ny-1 c i1=3*(i-1)+1+3*nx*(j1-1) i2=i1+1 i3=i1+2 r=u(i1) g=u(i2) b=u(i3) c r=amax1(r,0.0) g=amax1(g,0.0) b=amax1(b,0.0) r=amin1(r,1.0) g=amin1(g,1.0) b=amin1(b,1.0) c if(i.eq.nx.and.j.eq.2) then write(10,10) r,g,b 10 format(f6.2,1x,f6.2,1x,f6.2) else write(10,20) r,g,b 20 format(f6.2,1x,f6.2,1x,f6.2,",") end if c 100 continue return end c-------------------------------------------- subroutine colorp3(nx,ny,u,br,bg,bb) c * bottom tuki * c-------------------------------------------- dimension u(1) c do 100 i=1,nx do 100 j=1,2 j1=j+ny-1 c i1=3*(i-1)+1+3*nx*(j1-1) i2=i1+1 i3=i1+2 r=u(i1) g=u(i2) b=u(i3) c r=amax1(r,0.0) g=amax1(g,0.0) b=amax1(b,0.0) r=amin1(r,1.0) g=amin1(g,1.0) b=amin1(b,1.0) c write(10,10) r,g,b 10 format(f8.2,1x,f8.2,1x,f8.2,",") c 100 continue c----( bottom chusin color ) write(10,120) br,bg,bb 120 format(f8.2,1x,f8.2,1x,f8.2) c return end c------------------------------------- subroutine coordp1(nx,ny,p) c------------------------------------- dimension p(1) c do 100 i=1,nx c i1=3*(i-1)+1+3*nx*(ny-1) i2=i1+1 i3=i1+2 x=p(i1) y=p(i2) z=p(i3) c if(i.eq.nx) then write(10,10) x,y,z 10 format(f9.2,1x,f9.2,1x,f9.2) else write(10,20) x,y,z 20 format(f9.2,1x,f9.2,1x,f9.2,",") end if c 100 continue return end c------------------------------------- subroutine coordp2(nx,ny,p) c------------------------------------- dimension p(1) c do 100 i=1,nx do 100 j=1,2 j1=j+ny-1 c i1=3*(i-1)+1+3*nx*(j1-1) i2=i1+1 i3=i1+2 x=p(i1) y=p(i2) z=p(i3) c if(i.eq.nx.and.j.eq.2) then write(10,10) x,y,z 10 format(f9.2,1x,f9.2,1x,f9.2) else write(10,20) x,y,z 20 format(f9.2,1x,f9.2,1x,f9.2,",") end if c 100 continue return end c------------------------------------------- subroutine coordp3(nx,ny,p,rd) c * bottom tuki * c------------------------------------------- dimension p(1) c do 100 i=1,nx do 100 j=1,2 j1=j+ny-1 c i1=3*(i-1)+1+3*nx*(j1-1) i2=i1+1 i3=i1+2 x=p(i1) y=p(i2) z=p(i3) c write(10,20) x,y,z 20 format(f9.2,1x,f9.2,1x,f9.2,",") c 100 continue c c----( bottom no chusin ) c i1=3*nx*2+1 i2=i1+1 i3=i1+2 x=p(i1) y=p(i2) z=p(i3) c write(10,110) x,y,z 110 format(f9.2,1x,f9.2,1x,f9.2) c return end c------------------------------------------- subroutine coordp4(nx,ny,q,rd,k,p1) c * bottom tuki * c------------------------------------------- dimension q(1) dimension p1(1) c do 100 i=1,nx do 100 j=1,2 c i1=3*(i-1)+1+3*nx*(j-1) i2=i1+1 i3=i1+2 x=q(i1) y=q(i2) z=q(i3) c write(10,20) x,y,z 20 format(f9.2,1x,f9.2,1x,f9.2,",") c 100 continue c c----( bottom no chusin ) x=p1(1) y=p1(2) z=p1(3) c write(10,30) x,y,z 30 format(f9.2,1x,f9.2,1x,f9.2,",") c c----( cap no chusin ) x=p1(4) y=p1(5) z=p1(6) c write(10,40) x,y,z 40 format(f9.2,1x,f9.2,1x,f9.2) c return end c------------------------------- subroutine coorden1(nx) c------------------------------- dimension ia(9000) write(10,21) 21 format("]") write(10,22) 22 format("}") write(10,23) 23 format("colorPerVertex TRUE") write(10,24) 24 format("coordIndex [") c nx2=nx*2 do 100 i=1,nx2 i1=1+4*(i-1) i2=2+4*(i-1) i3=3+4*(i-1) i4=4+4*(i-1) ia(i1)=0+i-1 ia(i2)=1+i-1 ia(i3)=2+i-1 ia(i4)=-1 100 continue c jj=4*(nx-1)*2-1 write(10,102) (ia(i),i=1,jj) 102 format(8(3i2,i3),2i2,2i3,i2,3i3,88i3,2i3,i4,i3,i3,2i4,i3, 1 398(3i4,i3),/500(3i4,i3),/2i4,i5,i3,i4,2i5,i3, 2 /500(3i5,i3),/500(3i5,i3),/500(3i5,i3),/500(3i5,i3)) return end c----------------------------------- subroutine coorden2(nx) c * bottom tuki * c----------------------------------- dimension ia(9000) write(10,21) 21 format("]") write(10,22) 22 format("}") write(10,23) 23 format("colorPerVertex TRUE") write(10,24) 24 format("coordIndex [") c nx2=nx*2 do 100 i=1,nx2 i1=1+4*(i-1) i2=2+4*(i-1) i3=3+4*(i-1) i4=4+4*(i-1) ia(i1)=0+i-1 ia(i2)=1+i-1 ia(i3)=2+i-1 ia(i4)=-1 100 continue c jj=4*(nx-1)*2 write(10,102) (ia(i),i=1,jj) 102 format(8(3i2,i3),2i2,2i3,i2,3i3,88i3,2i3,i4,i3,i3,2i4,i3, 1 398(3i4,i3),/500(3i4,i3),/2i4,i5,i3,i4,2i5,i3, 2 /500(3i5,i3),/500(3i5,i3),/500(3i5,i3),/500(3i5,i3)) c write(10,103) ia(jj-2),ia(jj-1),ia(1) 103 format(3i5," -1") write(10,104) ia(jj-1),ia(1),ia(2) 104 format(3i5," -1") c c---- bottom nx2=nx*2 do 200 i=1,nx i1=1+4*(i-1) i2=2+4*(i-1) i3=3+4*(i-1) i4=4+4*(i-1) ia(i1)=0+2*(i-1) ia(i2)=ia(i1)+2 ia(i3)=nx2 ia(i4)=-1 200 continue c jj=4*(nx-1) write(10,202) (ia(i),i=1,jj) 202 format(500(3i4,i3),/500(3i4,i3),/500(3i5,i3),/500(3i5,i3)) c write(10,203) ia(jj-2),ia(1),nx2 203 format(3i4," -1") c c---- cap nx2=nx*2 do 300 i=1,nx i1=1+4*(i-1) i2=2+4*(i-1) i3=3+4*(i-1) i4=4+4*(i-1) ia(i1)=1+2*(i-1) ia(i2)=ia(i1)+2 ia(i3)=nx2 ia(i4)=-1 300 continue c jj=4*(nx-1) write(10,302) (ia(i),i=1,jj) 302 format(500(3i4,i3),/500(3i4,i3),/500(3i5,i3),/500(3i5,i3)) c write(10,303) ia(jj-2),ia(1),nx2 303 format(3i4," -1") c return end c----------------------------------- subroutine coorden3(nx) c * bottom tuki * c * 1 color * c----------------------------------- dimension ia(9000) write(10,21) 21 format("]") write(10,22) 22 format("}") write(10,23) 23 format("colorPerVertex FALSE") write(10,24) 24 format("coordIndex [") c nx2=nx*2 do 100 i=1,nx2 i1=1+4*(i-1) i2=2+4*(i-1) i3=3+4*(i-1) i4=4+4*(i-1) ia(i1)=0+i-1 ia(i2)=1+i-1 ia(i3)=2+i-1 ia(i4)=-1 100 continue c jj=4*(nx-1)*2 write(10,102) (ia(i),i=1,jj) 102 format(8(3i2,i3),2i2,2i3,i2,3i3,88i3,2i3,i4,i3,i3,2i4,i3, 1 398(3i4,i3),/500(3i4,i3),/2i4,i5,i3,i4,2i5,i3, 2 /500(3i5,i3),/500(3i5,i3),/500(3i5,i3),/500(3i5,i3)) c write(10,103) ia(jj-2),ia(jj-1),ia(1) 103 format(3i5," -1") write(10,104) ia(jj-1),ia(1),ia(2) 104 format(3i5," -1") c c---- bottom nx2=nx*2 do 200 i=1,nx i1=1+4*(i-1) i2=2+4*(i-1) i3=3+4*(i-1) i4=4+4*(i-1) ia(i1)=0+2*(i-1) ia(i2)=ia(i1)+2 ia(i3)=nx2 ia(i4)=-1 200 continue c jj=4*(nx-1) write(10,202) (ia(i),i=1,jj) 202 format(500(3i4,i3),/500(3i4,i3),/500(3i5,i3),/500(3i5,i3)) c write(10,203) ia(jj-2),ia(1),nx2 203 format(3i4," -1") c c---- cap nx21=nx*2+1 do 300 i=1,nx i1=1+4*(i-1) i2=2+4*(i-1) i3=3+4*(i-1) i4=4+4*(i-1) ia(i1)=1+2*(i-1) ia(i2)=ia(i1)+2 ia(i3)=nx21 ia(i4)=-1 300 continue c jj=4*(nx-1) write(10,302) (ia(i),i=1,jj) 302 format(500(3i4,i3),/500(3i4,i3),/500(3i5,i3),/500(3i5,i3)) c write(10,303) ia(jj-2),ia(1),nx21 303 format(3i4," -1") c return end c---------------------------------- subroutine coorden5(nx,ny1) c---------------------------------- dimension ia(9000) write(10,21) 21 format("]") write(10,22) 22 format("}") write(10,23) 23 format("colorPerVertex TRUE") write(10,24) 24 format("coordIndex [") c cc nx2=(nx-1)*2 nx2=nx*2 do 200 j=1,ny1 c do 100 i=1,nx2 i1=1+4*(i-1) i2=2+4*(i-1) i3=3+4*(i-1) i4=4+4*(i-1) ia(i1)=(j-1)*2*nx+0+i-1 ia(i2)=(j-1)*2*nx+1+i-1 ia(i3)=(j-1)*2*nx+2+i-1 ia(i4)=-1 100 continue c jj=4*(nx-1)*2-1 write(10,102) (ia(i),i=1,jj) c 102 format(8(3i2,i3),2i2,2i3,i2,3i3,88i3,2i3,i4,i3,i3,2i4,i3, c 1 398(3i4,i3),/500(3i4,i3),/2i4,i5,i3,i4,2i5,i3, 102 format(250(3i7,i3),/250(3i7,i3),/250(3i7,i3),/250(3i7,i3), 1 /250(3i7,i3),/250(3i7,i3),/250(3i7,i3),/250(3i7,i3), 2 /250(3i7,i3),/250(3i7,i3),/250(3i7,i3),/250(3i7,i3), 3 /250(3i7,i3),/250(3i7,i3),/250(3i7,i3),/250(3i7,i3), 4 /250(3i7,i3),/250(3i7,i3),/250(3i7,i3),/250(3i7,i3), 5 /250(3i7,i3),/250(3i7,i3),/250(3i7,i3),/250(3i7,i3), 6 /250(3i7,i3),/250(3i7,i3),/250(3i7,i3),/250(3i7,i3), 7 /250(3i7,i3),/250(3i7,i3),/250(3i7,i3),/250(3i7,i3), 8 /250(3i7,i3),/250(3i7,i3),/250(3i7,i3),/250(3i7,i3)) c 200 continue return end c------------------------ subroutine ifaceen c------------------------ write(10,20) 20 format("]") write(10,21) 21 format("solid FALSE",/"}") write(10,22) 22 format("appearance Appearance {") write(10,23) 23 format("material Material {") write(10,24) 24 format("transparency 0.000000",/"}",/"}") write(10,25) 25 format("}") return end c-------------------------- subroutine faceen2 c-------------------------- write(10,20) 20 format("]") write(10,21) 21 format("solid FALSE",/"}") write(10,25) 25 format("}") return end c-------------------------------------- subroutine pointlight(px,py,pz) c-------------------------------------- write(10,20) 20 format("PointLight {") write(10,21) px,py,pz 21 format("location ",f8.2,f8.2,f8.2) write(10,22) 22 format("}") return end c----------------------------------------- subroutine directlight(px,py,pz) c----------------------------------------- write(10,20) 20 format("DirectionalLight {") write(10,21) px,py,pz 21 format("direction ",f8.2,f8.2,f8.2) write(10,22) 22 format("}") return end c------------------------------------- subroutine lineset1(nx,ny,u,p) c------------------------------------- dimension p(1),u(1) ny1=ny-1 c do 10 j=1,ny call line1(nx,j,u,p) 10 continue c return end c------------------------------------- subroutine line1(nx,ny,u,p) c------------------------------------- dimension p(1),u(1) c call ilinebe call colorp1(nx,ny,u) call coordbe call coordp1(nx,ny,p) call ilineen(nx) c return end c------------------------ subroutine ilinebe c------------------------ write(10,21) 21 format("Shape {") write(10,22) 22 format("geometry IndexedLineSet {") write(10,23) 23 format("color Color {") write(10,24) 24 format("color [") return end c------------------------------- subroutine ilineen(nx) c------------------------------- dimension ia(9000) c do 10 i=1,9000 ia(i)=i-1 10 continue c write(10,21) 21 format("]") write(10,22) 22 format("}") write(10,24) 24 format("coordIndex [") c write(10,25) (ia(i),i=1,nx) 25 format(10i2,90i3,900i4/(2000i5)) c write(10,27) 27 format("]") write(10,29) write(10,29) 29 format("}") return end c---------------------------------------- subroutine imagetx(nx,ny,p,ch1,n) c---------------------------------------- dimension p(1),u(1) ny1=ny-1 c do 10 j=1,ny1 call imagea(nx,j,p,ch1,n) 10 continue c return end c--------------------------------------- subroutine imagea(nx,ny,p,ch1,n) c--------------------------------------- dimension p(1),u(1) c call imagebe(ch1,n) call coordp2(nx,ny,p) call coorden1(nx) call ifaceen c return end c------------------------------- subroutine imagebe(chr,n) c------------------------------- character chr*15,ica*15,ich(15)*1 equivalence (ica,ich) ica=chr write(10,21) 21 format("Shape {") write(10,22) 22 format("appearance Appearance {") write(10,23) 23 format("texture ImageTexture {") 24 write(10,*) "url",'"',(ich(i),i=1,n),'"' write(10,25) 25 format("}",/"}") write(10,26) 26 format("geometry IndexedFaceSet {") write(10,27) 27 format("coord Coordinate {") write(10,28) 28 format("point [") return end c---------------------------------------- subroutine movietx(nx,ny,p,ch1,n) c---------------------------------------- dimension p(1),u(1) ny1=ny-1 c do 10 j=1,ny1 call movie1(nx,j,p,ch1,n) 10 continue c return end c--------------------------------------- subroutine movie1(nx,ny,p,ch1,n) c--------------------------------------- dimension p(1),u(1) c call moviebe(ch1,n) call coordp2(nx,ny,p) call coorden1(nx) call ifaceen c return end c------------------------------- subroutine moviebe(chr,n) c------------------------------- character chr*15,ica*15,ich(15)*1 equivalence (ica,ich) ica=chr write(10,21) 21 format("Shape {") write(10,22) 22 format("appearance Appearance {") write(10,23) 23 format("texture MovieTexture {") 24 write(10,*) "url",'"',(ich(i),i=1,n),'"' write(10,30) cc 30 format("loop TRUE",/"startTime 0") 30 format("loop FALSE",/"startTime 0",/"stopTime 0") write(10,25) 25 format("}",/"}") write(10,26) 26 format("geometry IndexedFaceSet {") write(10,27) 27 format("coord Coordinate {") write(10,28) 28 format("point [") return end c----------------------------------------------------------- subroutine defusem(nx,ny,u,p,x1,y1,z1,x2,y2,z2,r,n) c---------------------------------------------------------- dimension p(1),u(1) ny1=ny-1 c do 10 j=1,ny1 call defusem1(nx,j,u,p,x1,y1,z1,x2,y2,z2,r,n) 10 continue c return end c----------------------------------------------------------- subroutine defusem1(nx,ny,u,p,x1,y1,z1,x2,y2,z2,r,n) c----------------------------------------------------------- dimension p(1),u(1) c call defface call colorp2(nx,ny,u) call coordbe call coordp2(nx,ny,p) call coorden1(nx) call ifaceen c do 10 i=1,n call useface(x1,y1,z1,x2,y2,z2,r,i) 10 continue c return end c-------------------------- subroutine defface c-------------------------- write(10,21) 21 format("DEF Face Shape {") write(10,22) 22 format("geometry IndexedFaceSet {") write(10,23) 23 format("color Color {") write(10,24) 24 format("color [") return end c------------------------------------------------ subroutine useface(x1,y1,z1,x2,y2,z2,r,i) c------------------------------------------------ write(10,21) 21 format("Transform {") write(10,22) x1,y1,z1,r 22 format("rotation ",f6.2,f6.2,f6.2,f6.2) write(10,23) x2,y2*i,z2 23 format("translation ",f9.2,f9.2,f9.2) write(10,24) 24 format("children USE Face",/"}") return end c----------------------------------------- subroutine sphere2(p,n1,rd,r,g,b) c----------------------------------------- dimension p(1) c do 20 i=1,n1 c i1=3*(i-1)+1 i2=i1+1 i3=i1+2 x=p(i1) y=p(i2) z=p(i3) c call sphere(x,y,z,r,g,b,rd) c 20 continue c return end c---------------------------------------- subroutine sphere(x,y,z,r,g,b,rd) c---------------------------------------- write(10,20) 20 format("Transform {") write(10,21) x,y,z 21 format("translation ",f9.2,f9.2,f9.2) write(10,22) 22 format("children [") write(10,23) 23 format("Shape {") write(10,24) 24 format("appearance Appearance {") write(10,25) 25 format("material Material {") write(10,26) 26 format("transparency 0.000000") write(10,27) r,g,b 27 format("diffuseColor ",f5.2,f5.2,f5.2) write(10,28) 28 format("}",/"}") write(10,29) 29 format("geometry Sphere {") write(10,30) rd 30 format("radius ",f6.2) write(10,31) 31 format("}",/"}",/"]",/"}") c return end c------------------------------------------------------- subroutine cylinder1(p,n1,rx,ry,rz,dd,rd,h,r,g,b) c------------------------------------------------------- dimension p(1) c do 20 i=1,n1 c i1=3*(i-1)+1 i2=i1+1 i3=i1+2 x=p(i1) y=p(i2) z=p(i3) c call cylinder(x,y,z,rx,ry,rz,dd,rd,h,r,g,b) c 20 continue c return end c----------------------------------------------- subroutine cylinderm2(p,n1,rd,r,g,b,h) c----------------------------------------------- dimension p(1) rx=0.0 ry=0.0 rz=0.0 dd=0.78539 c do 20 i=1,n1 c i1=3*(i-1)+1 i2=i1+1 i3=i1+2 x=p(i1) y=p(i2) z=p(i3) c x1=p(i1+3) y1=p(i2+3) z1=p(i3+3) c x2=abs(x1-x) y2=abs(y1-y) z2=abs(z1-z) c c---<< x houkou >> if(x2.gt.y2.and.x2.gt.z2) then rz=-1 if(z2.gt.y2) then rx=1 end if end if c---<< y houkou >> c if(y2.gt.z2.and.y2.gt.x2) then c rx=-1 c else c rx=1 c end if c end if c---<< z houkou >> if(z2.gt.x2.and.z2.gt.y2) then if(x2.lt.y2) then rx=-1 else rx=1 end if end if c call cylinder(x,y,z,rx,ry,rz,dd,rd,h,r,g,b) c call def1(x,y,z,rx,ry,rz,dd,rd,h,r,g,b) c 20 continue c 99 continue c return end c----------------------------------------------------------- subroutine cylinder(cx,cy,cz,rx,ry,rz,dd,rd,t,r,g,b) c----------------------------------------------------------- write(10,20) 20 format("Transform {") write(10,21) rx,ry,rz,dd 21 format("rotation ",f7.2,f7.2,f7.2,f7.4) write(10,22) cx,cy,cz 22 format("translation ",f7.2,f7.2,f7.2) write(10,23) 23 format("children Shape {") write(10,24) 24 format("appearance Appearance {") write(10,25) 25 format("material Material {") write(10,26) 26 format("transparency 0.000000") write(10,27) r,g,b 27 format("diffuseColor",f5.2,f5.2,f5.2) write(10,28) 28 format("}",/"}") write(10,29) 29 format("geometry Cylinder {") write(10,30) rd 30 format("radius",f7.2) write(10,31) t 31 format("height",f7.2) write(10,32) 32 format("}",/"}",/"}") c return end c------------------------------------------------------- subroutine defl(cx,cy,cz,rx,ry,rz,dd,rd,t,r,g,b) c------------------------------------------------------- write(10,20) 20 format("Transform {") write(10,21) rx,ry,rz,dd 21 format("rotation ",f9.2,f9.2,f9.2,f9.2) write(10,22) cx,cy,cz 22 format("translation ",f9.2,f9.2,f9.2) write(10,23) 23 format("children DEF def1 Shape {") write(10,24) 24 format("appearance Appearance {") write(10,25) 25 format("material Material {") write(10,26) 26 format("transparency 0.000000") write(10,27) r,g,b 27 format("diffuseColor",f5.2,f5.2,f5.2) write(10,28) 28 format("}",/"}") write(10,29) 29 format("geometry Cylinder {") write(10,30) rd 30 format("radius",f8.2) write(10,31) t 31 format("height",f8.2) write(10,32) 32 format("}",/"}",/"}") c return end c------------------------------------------------ subroutine use1(cx,cy,cz,rx,ry,rz,dd) c------------------------------------------------ write(10,20) 20 format("Transform {") write(10,21) rx,ry,rz,dd 21 format("rotation ",f9.2,f9.2,f9.2,f9.2) write(10,22) cx,cy,cz 22 format("translation ",f9.2,f9.2,f9.2) write(10,23) 23 format("children USE def1") write(10,32) 32 format("}") c return end c --------------------- c subroutine pixel1 c --------------------- subroutine pixel1(nx,ny,xb,yb,xl,yl,ipx0,ico,icc,zcc, 1 vmin,vmax,u) c c nx: number in x-direction c nx: number in x-direction c x0: x-value at the left and bottom point c y0: y-value at the left and bottom point c xl: length in x-direction c yl: length in y-direction c ipx0: 3 for color, 4 for transparency c ico: color code c icc: selection of plane, i c : 1 for xy plane, s2 for xz plane, 3 for yz plane c zcc: position of selected plane c u(nx*ny): input data c dimension u(1) nxy=nx*ny c do 10 j=1,ny do 10 i=1,nx i1=i+nx*(j-1) u(i1)=(u(i1)-vmin)/(vmax-vmin) 10 continue c call pixebe(nx,ny,ipx0) if(ipx0.eq.3) call pixeco(nx,ny,ico,u) if(ipx0.eq.4) call pixect(nx,ny,ico,u) call pixeen(xb,yb,xl,yl,icc,zcc) c return end c --------------------- c subroutine pixel2 c --------------------- subroutine pixel2(nx,ny,xb,yb,xl,yl,ipx0,ico,icc,zcc,u) c c nx: number in x-direction c nx: number in x-direction c x0: x-value at the left and bottom point c y0: y-value at the left and bottom point c xl: length in x-direction c yl: length in y-direction c ipx0: 3 for color, 4 for transparency c ico: color code c icc: selection of plane, i c : 1 for xy plane, s2 for xz plane, 3 for yz plane c zcc: position of selected plane c u(nx*ny): input data c dimension u(1) nxy=nx*ny c call pixebe(nx,ny,ipx0) if(ipx0.eq.3) call pixeco(nx,ny,ico,u) if(ipx0.eq.4) call pixect(nx,ny,ico,u) call pixeen(xb,yb,xl,yl,icc,zcc) c return end c --------------------- c subroutine pixebe c --------------------- subroutine pixebe(nx,ny,ipx0) write(10,21) 21 format("Shape {") write(10,22) 22 format("appearance Appearance {") write(10,23) 23 format("texture PixelTexture {") write(10,24) nx,ny,ipx0 24 format("image ",2i5,i2) return end c --------------------- c subroutine pixeco c --------------------- subroutine pixeco(nx,ny,ico,u) dimension u(1),ipx(6) r=1.0 g=1.0 b=1.0 do 100 j=1,ny do 100 i=1,nx i1=i+nx*(j-1) x1=u(i1) c a0=1.00 a0=1.20 if(ico.eq.2) go to 32 if(ico.eq.3) go to 33 if(ico.eq.4) go to 34 if(ico.eq.5) go to 35 if(ico.eq.6) go to 36 if(ico.eq.7) go to 37 c x1=2.0*x1 r=x1 g=2.0-x1 b=0.0 go to 40 32 continue r=x1 g=x1 b=x1 go to 40 33 continue x1=7.0*x1 b=x1 if(x1.ge.2.0) b=3.0-x1 if(x1.ge.5.0) b=0.5*(x1-5.0) g=x1-1.0 if(x1.ge.4.0) g=5.0-x1 if(x1.ge.5.0) g=0.5*(x1-5.0) r=x1-3.0 go to 40 34 continue x1=10.0*x1 b=x1 if(x1.ge.3.0) b=4.0-x1 if(x1.ge.6.0) b=0.25*(x1-6.0) g=x1-2.0 if(x1.ge.5.0) g=6.0-x1 if(x1.ge.6.0) g=0.25*(x1-6.0) r=x1 if(x1.ge.1.0) r=2.0-x1 if(x1.ge.4.0) r=x1-4.0 go to 40 35 continue a0=2.0 y1=x1*(a0+4.0) b=a0-y1 if(y1.ge.a0) b=y1-(a0+1.0) r=a0+3.0-y1 g=a0+1.0-y1 if(y1.ge.a0+2.0) g=y1-(a0+3.0) go to 40 36 continue a0=2.0 y1=x1*(a0+4.0) r=a0-y1 if(y1.ge.a0) r=y1-(a0+1.0) b=a0+3.0-y1 g=a0+1.0-y1 if(y1.ge.a0+2.0) g=y1-(a0+3.0) go to 40 37 continue x1=10.0*x1 b=x1 if(x1.ge.3.0) b=4.0-x1 if(x1.ge.6.0) b=0.25*(x1-6.0) g=x1-2.0 if(x1.ge.5.0) g=6.0-x1 if(x1.ge.6.0) g=0.25*(x1-6.0) r=x1 if(x1.ge.1.0) r=2.0-x1 if(x1.ge.4.0) r=x1-4.0 go to 40 c 40 continue r=amax1(r,0.0) g=amax1(g,0.0) b=amax1(b,0.0) r=amin1(r,1.0) g=amin1(g,1.0) b=amin1(b,1.0) c ir=r*256.0 ig=g*256.0 ib=b*256.0 ir=max(0,ir) ir=min(255,ir) ig=max(0,ig) ig=min(255,ig) ib=max(0,ib) ib=min(255,ib) ir1=ir/16 ir2=ir-16*ir1 ig1=ig/16 ig2=ig-16*ig1 ib1=ib/16 ib2=ib-16*ib1 ipx(1)=ir1 ipx(2)=ir2 ipx(3)=ig1 ipx(4)=ig2 ipx(5)=ib1 ipx(6)=ib2 c write(10,20) ipx c 20 format("0x",6z1,1x) 20 format(" 0x",6z1) c 100 continue return end c --------------------- c subroutine pixect c --------------------- subroutine pixect(nx,ny,ico,u) dimension u(1),ipx(8) r=1.0 g=1.0 b=1.0 do 100 j=1,ny do 100 i=1,nx i1=i+nx*(j-1) x1=u(i1) c a0=1.00 a0=1.20 if(ico.eq.2) go to 32 if(ico.eq.3) go to 33 if(ico.eq.4) go to 34 if(ico.eq.5) go to 35 if(ico.eq.6) go to 36 if(ico.eq.7) go to 37 c x1=2.0*x1 r=x1 g=2.0-x1 b=0.0 go to 40 32 continue r=x1 g=x1 b=x1 go to 40 33 continue x1=7.0*x1 b=x1 if(x1.ge.2.0) b=3.0-x1 if(x1.ge.5.0) b=0.5*(x1-5.0) g=x1-1.0 if(x1.ge.4.0) g=5.0-x1 if(x1.ge.5.0) g=0.5*(x1-5.0) r=x1-3.0 go to 40 34 continue x1=10.0*x1 b=x1 if(x1.ge.3.0) b=4.0-x1 if(x1.ge.6.0) b=0.25*(x1-6.0) g=x1-2.0 if(x1.ge.5.0) g=6.0-x1 if(x1.ge.6.0) g=0.25*(x1-6.0) r=x1 if(x1.ge.1.0) r=2.0-x1 if(x1.ge.4.0) r=x1-4.0 go to 40 35 continue a0=2.0 y1=x1*(a0+4.0) b=a0-y1 if(y1.ge.a0) b=y1-(a0+1.0) r=a0+3.0-y1 g=a0+1.0-y1 if(y1.ge.a0+2.0) g=y1-(a0+3.0) go to 40 36 continue a0=2.0 y1=x1*(a0+4.0) r=a0-y1 if(y1.ge.a0) r=y1-(a0+1.0) b=a0+3.0-y1 g=a0+1.0-y1 if(y1.ge.a0+2.0) g=y1-(a0+3.0) go to 40 37 continue x1=10.0*x1 b=x1 if(x1.ge.3.0) b=4.0-x1 if(x1.ge.6.0) b=0.25*(x1-6.0) g=x1-2.0 if(x1.ge.5.0) g=6.0-x1 if(x1.ge.6.0) g=0.25*(x1-6.0) r=x1 if(x1.ge.1.0) r=2.0-x1 if(x1.ge.4.0) r=x1-4.0 go to 40 c 40 continue r=amax1(r,0.0) g=amax1(g,0.0) b=amax1(b,0.0) r=amin1(r,1.0) g=amin1(g,1.0) b=amin1(b,1.0) c ir=r*256.0 ig=g*256.0 ib=b*256.0 ir=max(0,ir) ir=min(255,ir) ig=max(0,ig) ig=min(255,ig) ib=max(0,ib) ib=min(255,ib) ir1=ir/16 ir2=ir-16*ir1 ig1=ig/16 ig2=ig-16*ig1 ib1=ib/16 ib2=ib-16*ib1 ipx(1)=ir1 ipx(2)=ir2 ipx(3)=ig1 ipx(4)=ig2 ipx(5)=ib1 ipx(6)=ib2 c it=u(i1)*256.0 t=max(0,it) it=min(255,it) it1=it/16 it2=it-16*it1 ipx(7)=it1 ipx(8)=it2 c write(10,20) ipx c 20 format("0x",8z1,1x) 20 format(" 0x",8z1) c 100 continue return end c --------------------- c subroutine pixeen c --------------------- subroutine pixeen(xb,yb,xl,yl,icc,zcc) x1=xb x2=xb+xl y1=yb y2=yb+yl c write(10,21) 21 format(" }}") write(10,22) 22 format("geometry IndexedFaceSet {") write(10,23) 23 format("coord Coordinate { point [") c write(10,24) nx,ny,ipx0 c 24 format("image ",2i5,i2) c if (icc.eq.2) go to 52 if (icc.eq.3) go to 53 c z=zcc write(10,31) x1,y1,z write(10,31) x2,y1,z write(10,31) x2,y2,z write(10,32) x1,y2,z 31 format(3f6.2,",") 32 format(3f6.2) go to 60 c 52 continue x=zcc write(10,31) x,x1,y1 write(10,31) x,x2,y1 write(10,31) x,x2,y2 write(10,32) x,x1,y2 go to 60 c 53 continue y=zcc write(10,31) x1,y,y1 write(10,31) x2,y,y1 write(10,31) x2,y,y2 write(10,32) x1,y,y2 60 continue c write(10,41) 41 format(" ]}") write(10,42) 42 format("coordIndex [ 0 1 2 3 ]") write(10,43) 43 format("solid FALSE"/"texCoord TextureCoordinate {") write(10,44) 44 format("point ["/"0.0 0.0,"/"1.0 0.0,"/"1.0 1.0,") write(10,45) 45 format("0.0 1.0"/"]}"/"texCoordIndex [ 0 1 2 3 ]") write(10,46) 46 format("}}") c return end c --------------------- SUBROUTINE QUANT1(X,Y,Z,HX,HY,HZ,MM,NX,NY,NZ,NXP,AA,F,P) c --------------------- DIMENSION F(1),P(1),AA(1),B(3) C xo=x yo=y zo=z if(yo.lt.0.0) then y=-y z=-z endif NX1=NX+1 NX2=NX+2 NY1=NY+1 NZ1=NZ+1 NY2=NY+2 NZ2=NZ+2 N1=NX+2 N2=N1*(NY+2) N3=N2*(NZ+2) N11=N1+1 N12=N2+1 N21=N2+N1 N22=N21+1 ARU=AA(2) AR1=AA(3) B0=AA(8) C XI=X/HX+0.5*FLOAT(NX2+1-2*NXP) c XJ=Y/HY+0.5*FLOAT(ny2+1) xj=y/hy+1.5 c XK=Z/HZ+1.5 XK=Z/HZ+0.5*FLOAT(nz2+1) I=IFIX(XI) J=IFIX(XJ) K=IFIX(XK) IF(I.LT.1) I=1 IF(J.LT.1) J=1 IF(K.LT.1) K=1 IF(I.GT.NX1) I=NX1 IF(J.GT.NY1) J=NY1 IF(K.GT.NZ1) K=NZ1 XI=XI-FLOAT(I) XJ=XJ-FLOAT(J) XK=XK-FLOAT(K) XI1=1.0-XI XJ1=1.0-XJ XK1=1.0-XK I1=I+N1*(J-1)+N2*(K-1)-N3 L1=8*(MM-1) DO 10 M=1,3 I1=I1+N3 I2=M+L1 P(I2)=XI*XJ*XK*F(I1+N22)+XI*XJ*XK1*F(I1+N11) & +XI*XJ1*XK*F(I1+N12)+XI1*XJ*XK*F(I1+N21) & +XI*XJ1*XK1*F(I1+1)+XI1*XJ*XK1*F(I1+N1) & +XI1*XJ1*XK*F(I1+N2)+XI1*XJ1*XK1*F(I1) P(I2+5)=P(I2) 10 CONTINUE C AR2=AR1*AR1 RO2=X*X+Y*Y+Z*Z RO5=RO2*RO2*SQRT(RO2) X2=0.0 IF(RO2.GT.AR2) X2=RO2/AR2-1.0 X2=ARU*X2*X2 X2=X2/(X2+1.0) B(1)=-3.0*B0*X*Z/RO5 B(2)=-3.0*B0*Y*Z/RO5 B(3)=B0*(RO2-3.0*Z*Z)/RO5 L1=5+8*(MM-1) DO 20 M=1,3 I1=M+L1 P(I1)=P(I1)*X2+B(M)*(1.0-X2) 20 CONTINUE if(yo.lt.0.0) then i1=1+L1 p(i1)=-p(i1) endif x=xo y=yo z=zo RETURN END c --------------------- SUBROUTINE QUANT2(X,Y,Z,HX,HY,HZ,MM,NX,NY,NZ,NXP,AA,F,P) c --------------------- DIMENSION F(1),P(1),AA(1),B(3) C xo=x yo=y zo=z if(yo.lt.0.0) then y=-y z=-z endif NX1=NX+1 NX2=NX+2 NY1=NY+1 ny2=ny+2 NZ1=NZ+1 nz2=nz+2 N1=NX+2 N2=N1*(NY+2) N3=N2*(NZ+2) N11=N1+1 N12=N2+1 N21=N2+N1 N22=N21+1 ARU=AA(2) AR1=AA(3) B0=AA(8) C XI=X/HX+0.5*FLOAT(NX2+1-2*NXP) XJ=Y/HY+1.5 c XK=Z/HZ+1.5 xk=z/hz+0.5*float(nz2+1) I=IFIX(XI) J=IFIX(XJ) K=IFIX(XK) IF(I.LT.1) I=1 IF(J.LT.1) J=1 IF(K.LT.1) K=1 IF(I.GT.NX1) I=NX1 IF(J.GT.NY1) J=NY1 IF(K.GT.NZ1) K=NZ1 XI=XI-FLOAT(I) XJ=XJ-FLOAT(J) XK=XK-FLOAT(K) XI1=1.0-XI XJ1=1.0-XJ XK1=1.0-XK I1=I+N1*(J-1)+N2*(K-1)-N3 L1=8*(MM-1) DO 10 M=1,3 I1=I1+N3 I2=M+L1 P(I2)=XI*XJ*XK*F(I1+N22)+XI*XJ*XK1*F(I1+N11) & +XI*XJ1*XK*F(I1+N12)+XI1*XJ*XK*F(I1+N21) & +XI*XJ1*XK1*F(I1+1)+XI1*XJ*XK1*F(I1+N1) & +XI1*XJ1*XK*F(I1+N2)+XI1*XJ1*XK1*F(I1) P(I2+5)=P(I2) 10 CONTINUE C AR2=AR1*AR1 RO2=X*X+Y*Y+Z*Z RO5=RO2*RO2*SQRT(RO2) X2=0.0 IF(RO2.GT.AR2) X2=RO2/AR2-1.0 X2=ARU*X2*X2 X2=X2/(X2+1.0) B(1)=-3.0*B0*X*Z/RO5 B(2)=-3.0*B0*Y*Z/RO5 B(3)=B0*(RO2-3.0*Z*Z)/RO5 L1=5+8*(MM-1) DO 20 M=1,3 I1=M+L1 P(I1)=P(I1)*X2+B(M)*(1.0-X2) 20 CONTINUE if(yo.lt.0.0) then i1=1+L1 p(i1)=-p(i1) endif x=xo y=yo z=zo RETURN END c --------------------- SUBROUTINE AINTE1(IA,AA,F,P) c --------------------- parameter(nnx=102,nny=102,nnz=2) PARAMETER(nline=9900,nlt=nline*3) DIMENSION IA(1),AA(1),F(1),P(1),Q(16),po(nlt) common app(nnx,nny,nnz),ie(nnx,nny,nnz) C IA 1 2 3 4 5 6 7 8 9 10 C NX NY MX MY MZ NXP MI MO C AA 1 2 3 4 5 6 7 8 9 10 C TH0 ARU AR1 ARB XXL YYL ZZL B0 DL SIGN C 11 12 13 14 15 16 17 18 19 20 C DIRE BOUN GX0 GY0 GXL GYL GTH GXMI GXMA EP1 C P(1)=P(1) PI=3.1415926 NX=IA(1) NY=IA(2) MX=IA(3) MY=IA(4) MZ=IA(5) NXP=IA(6) MI=IA(7) MO=IA(8) TH0=AA(1)*PI/180.0 ARU=AA(2) AR1=AA(3) ARB=AA(4) XXL=AA(5) YYL=AA(6) ZZL=AA(7) B0=AA(8) GX0=AA(13) GY0=AA(14) GXL=AA(15) GYL=AA(16) GTH=AA(17)*PI/180.0 GXMI=AA(18) GXMA=AA(19) GYMA=AA(21) EP1=AA(20) ccc waku kaku tokoro C x1=gxmi x2=0.0 x3=gxma y1=-gyma y2=0.0 y3=gyma z1=-gyma z2=0.0 z3=gyma c r=1.0 g=1.0 b=1.0 r=0.0 g=0.5 b=1.0 lasl=5 c call linebe(r,g,b) call linep2(x1,y2,z1,1,lasl) call linep2(x1,y3,z1,2,lasl) call linep2(x1,y3,z3,3,lasl) call linep2(x1,y2,z3,4,lasl) call linep2(x1,y2,z1,5,lasl) call lineen(lasl) call linebe(r,g,b) call linep2(x2,y2,z1,1,lasl) call linep2(x2,y3,z1,2,lasl) call linep2(x2,y3,z3,3,lasl) call linep2(x2,y2,z3,4,lasl) call linep2(x2,y2,z1,5,lasl) call lineen(lasl) call linebe(r,g,b) call linep2(x3,y2,z1,1,lasl) call linep2(x3,y3,z1,2,lasl) call linep2(x3,y3,z3,3,lasl) call linep2(x3,y2,z3,4,lasl) call linep2(x3,y2,z1,5,lasl) call lineen(lasl) c call linebe(r,g,b) call linep2(x1,y2,z1,1,lasl) call linep2(x1,y3,z1,1,lasl) call linep2(x3,y3,z1,1,lasl) call linep2(x3,y2,z1,1,lasl) call linep2(x1,y2,z1,5,lasl) call lineen(lasl) call linebe(r,g,b) call linep2(x1,y2,z2,1,lasl) call linep2(x1,y3,z2,1,lasl) call linep2(x3,y3,z2,1,lasl) call linep2(x3,y2,z2,1,lasl) call linep2(x1,y2,z2,5,lasl) call lineen(lasl) call linebe(r,g,b) call linep2(x1,y2,z3,1,lasl) call linep2(x1,y3,z3,1,lasl) call linep2(x3,y3,z3,1,lasl) call linep2(x3,y2,z3,1,lasl) call linep2(x1,y2,z3,5,lasl) call lineen(lasl) c call linebe(r,g,b) call linep2(x1,y2,z2,1,lasl) call linep2(x1,y3,z2,1,lasl) call linep2(x2,y3,z2,1,lasl) call linep2(x2,y2,z2,1,lasl) call linep2(x1,y2,z2,5,lasl) call lineen(lasl) c c call linebe(r,g,b) call linep2(x1,-y2,z1,1,lasl) call linep2(x1,-y3,z1,2,lasl) call linep2(x1,-y3,z3,3,lasl) call linep2(x1,-y2,z3,4,lasl) call linep2(x1,-y2,z1,5,lasl) call lineen(lasl) call linebe(r,g,b) call linep2(x2,-y2,z1,1,lasl) call linep2(x2,-y3,z1,2,lasl) call linep2(x2,-y3,z3,3,lasl) call linep2(x2,-y2,z3,4,lasl) call linep2(x2,-y2,z1,5,lasl) call lineen(lasl) call linebe(r,g,b) call linep2(x3,-y2,z1,1,lasl) call linep2(x3,-y3,z1,2,lasl) call linep2(x3,-y3,z3,3,lasl) call linep2(x3,-y2,z3,4,lasl) call linep2(x3,-y2,z1,5,lasl) call lineen(lasl) c call linebe(r,g,b) call linep2(x1,-y2,z1,1,lasl) call linep2(x1,-y3,z1,1,lasl) call linep2(x3,-y3,z1,1,lasl) call linep2(x3,-y2,z1,1,lasl) call linep2(x1,-y2,z1,5,lasl) call lineen(lasl) call linebe(r,g,b) call linep2(x1,-y2,z2,1,lasl) call linep2(x1,-y3,z2,1,lasl) call linep2(x3,-y3,z2,1,lasl) call linep2(x3,-y2,z2,1,lasl) call linep2(x1,-y2,z2,5,lasl) call lineen(lasl) call linebe(r,g,b) call linep2(x1,-y2,z3,1,lasl) call linep2(x1,-y3,z3,1,lasl) call linep2(x3,-y3,z3,1,lasl) call linep2(x3,-y2,z3,1,lasl) call linep2(x1,-y2,z3,5,lasl) call lineen(lasl) c call linebe(r,g,b) call linep2(x1,-y2,z2,1,lasl) call linep2(x1,-y3,z2,1,lasl) call linep2(x2,-y3,z2,1,lasl) call linep2(x2,-y2,z2,1,lasl) call linep2(x1,-y2,z2,5,lasl) call lineen(lasl) c c XX=999.0 XM=0.075 YM=0.075 X1=GX0 Y1=GY0-2.0*YM C CALL CHARA1(X1,Y1,0,0,-1,0,MI,GXMI,GXMA) C XC1=COS(GTH) XS1=SIN(GTH) Y=GY0+0.5*GYL c CALL PLOT(GX0,Y,3) X=GX0 Y=GY0+GYL c CALL PLOT(X,Y,2) X=GX0+GXL c CALL PLOT(X,Y,2) Y=GY0+0.5*GYL c CALL PLOT(X,Y,2) X=GX0 c CALL PLOT(X,Y,2) X=GX0+0.5*GYL*XC1 Y=GY0+0.5*GYL*XS1+0.5*GYL c CALL PLOT(X,Y,2) X=GX0+0.5*GYL*XC1+GXL c CALL PLOT(X,Y,2) X=GX0+GXL Y=GY0+0.5*GYL c CALL PLOT(X,Y,2) C y=gy0 c call plot(x,y,2) X=GX0 c CALL PLOT(X,Y,2) Y=GY0+0.5*GYL c CALL PLOT(X,Y,2) c XDD1=GXL*(0.0-GXMI)/(GXMA-GXMI) X=GX0+XDD1 Y=GY0+GYL c CALL PLOT(X,Y,3) Y=GY0+0.5*GYL c CALL PLOT(X,Y,2) X=X+0.5*GYL*XC1 Y=Y+0.5*GYL*XS1 c CALL PLOT(X,Y,2) c X=GX0+xdd1 Y=GY0+GYL*0.5 c CALL PLOT(X,Y,2) Y=GY0 c CALL PLOT(X,Y,2) c c call linee ccc waku kakiowatta! C ccc boundary?? condition ?? XL0=COS(TH0) XL01=XL0*1.0 XL=2.0*XL0 YL=XL NX1=NX+1 NX2=NX+2 NY1=NY+1 NY2=NY+2 N1=NX2 N2=N1*NY2 AR2=AR1*AR1 DX=XL/FLOAT(NX1) DY=YL/FLOAT(NY1) HX=XXL/FLOAT(MX+1) HY=YYL/FLOAT(MY+1) HZ=ZZL/FLOAT(MZ+1) XLMI=0.5*HX*FLOAT(2*NXP-MX-1) XLMA=0.5*HX*FLOAT(2*NXP+MX+1) DL0=HX*AA(9) DL=DL0 LAST=IFIX(XXL/DL) GYMI=-0.1*HY GYMA=AA(21) AMPY=0.5*(GXMA-GXMI)*GYL/(GYMA*GXL) c GZMI=-0.1*HZ gymi=-gyma gzmi=-gyma c GZMI GZMA ?? GZMA=GYMA AMP=GXL/(GXMA-GXMI) M1=MI M2=M1+1 M3=M2+1 M4=MI+8 M5=M4+1 M6=M5+1 c arcc0=hx*0.5 ardr0=hx*5.0 C cc polar cap kara start suru basyo wo kimeru JB1=NY2/2+1 do 100 jj=1,2 DO 100 J=JB1,NY2 cc DO 100 J=1,NY2 DO 100 I=1,NX2 X=0.5*DX*(2*I-NX2-1) Y=0.5*DY*(2*J-NY2-1) R=SQRT(X*X+Y*Y) IF(R.GT.XL01) GO TO 100 c if(x.gt.0.0) go to 100 c z ga + sika naiyouna.... Z=SQRT(1.0-R*R) R=ARB*SQRT(ARB) X1=X*R Y1=Y*R Z1=ARB*ARB-X1*X1-Y1*Y1 c zga 0.0 yori tiisaitoki 100 he IF(Z1.LE.0.0) GO TO 100 c IF(Z1.LE.0.0) GO TO 1030 if(jj.eq.1) Z1=SQRT(Z1) if(jj.eq.2) Z1=-SQRT(Z1) 1040 CALL QUANT1(X1,Y1,Z1,HX,HY,HZ,1,MX,MY,MZ,NXP,AA,F,Q) B1=SQRT(Q(M1)*Q(M1)+Q(M2)*Q(M2)+Q(M3)*Q(M3)) IF(B1.LT.1.0E-20) GO TO 100 AA(11)=1.0 X2=Q(M1)*X1+Q(M2)*Y1+Q(M3)*Z1 IF(X2.LT.0.0) AA(11)=-1.0 C X2=0.0 Y2=0.0 Z2=0.0 X3=0.0 Y3=0.0 Z3=0.0 GX1=GX0+XDD1+AMP*(X1+AMPY*Y1*XC1) cc GY1=GY0+0.5*GYL+AMP*AMPY*(Z1+Y1*XS1) GY1=GY0+0.5*GYL+AMP*AMPY*(Z1+Y1*XS1) c CALL PLOT(GX1,GY1,3) lin=1 li1=3*(lin-1)+1 li2=3*(lin-1)+2 li3=3*lin po(li1)=X1 po(li2)=Y1 po(li3)=Z1 DX2=0.0 DY2=0.0 DZ2=0.0 DR2=SQRT(DX2*DX2+DY2*DY2+DZ2*DZ2) DR2=AMAX1(dr2,1.0e-8) arcc=0.0 drar=0.0 c**** DO 30 III=1,LAST R=SQRT((X3-X2)**2+(Y3-Y2)**2+(Z3-Z2)**2) DL=R/EP1 IF(DL.GT.0.5*HX) DL=0.5*HX IF(DL.LT.DL0) DL=DL0 DL1=DL*AA(11) X2=X1+DL1*Q(M1)/B1 Y2=Y1+DL1*Q(M2)/B1 Z2=Z1+DL1*Q(M3)/B1 R2=SQRT(X2*X2+Y2*Y2+Z2*Z2) IF(R2.LT.0.5*ARB) GO TO 210 IF(X2.LT.GXMI.OR.X2.GT.GXMA) GO TO 200 IF(Y2.LT.GYMI.OR.Y2.GT.GYMA) GO TO 200 IF(Z2.LT.GZMI.OR.Z2.GT.GZMA) GO TO 200 CALL QUANT1(X2,Y2,Z2,HX,HY,HZ,2,MX,MY,MZ,NXP,AA,F,Q) B2=SQRT(Q(M4)*Q(M4)+Q(M5)*Q(M5)+Q(M6)*Q(M6)) IF(B2.LT.1.0E-20) GO TO 100 DO 32 II=1,1 DX3=0.5*DL1*(Q(M1)/B1+Q(M4)/B2) DY3=0.5*DL1*(Q(M2)/B1+Q(M5)/B2) DZ3=0.5*DL1*(Q(M3)/B1+Q(M6)/B2) DR3=SQRT(DX3*DX3+DY3*DY3+DZ3*DZ3) X3=X1+DX3 Y3=Y1+DY3 Z3=Z1+DZ3 R2=SQRT(X3*X3+Y3*Y3+Z3*Z3) IF(R2.LT.0.5*ARB) GO TO 210 IF(X3.LT.GXMI.OR.X3.GT.GXMA) GO TO 200 IF(Y3.LT.GYMI.OR.Y3.GT.GYMA) GO TO 200 IF(Z3.LT.GZMI.OR.Z3.GT.GZMA) GO TO 200 CALL QUANT1(X3,Y3,Z3,HX,HY,HZ,2,MX,MY,MZ,NXP,AA,F,Q) B2=SQRT(Q(M4)*Q(M4)+Q(M5)*Q(M5)+Q(M6)*Q(M6)) IF(B2.LT.1.0E-20) GO TO 100 32 CONTINUE C GX1=GX0+XDD1+AMP*(X3+AMPY*Y3*XC1) cc GY1=GY0+0.5*GYL+AMP*AMPY*(Z3+Y3*XS1) GY1=GY0+0.5*GYL+AMP*AMPY*(Z3+Y3*XS1) c CALL PLOT(GX1,GY1,2) c arc1=sqrt((dx2/dr2-dx3/dr3)**2+(dy2/dr2-dy3/dr3)**2+ 1 (dz2/dr2-dz3/dr3)**2) arcc=arcc+arc1 ardr=ardr+dr3 DX2=DX3 DY2=DY3 DZ2=DZ3 DR2=DR3 if(arcc.gt.arcc0.or.ardr.gt.ardr0) go to 326 go to 327 326 continue c arcc=0.0 ardr=0.0 lin=lin+1 c lin=iii+1 li1=3*(lin-1)+1 li2=3*(lin-1)+2 li3=3*lin po(li1)=X3 po(li2)=Y3 po(li3)=Z3 327 continue X1=X3 Y1=Y3 Z1=Z3 B1=B2 DO 34 II=1,8 Q(II)=Q(II+8) 34 CONTINUE c if(mod(iii,100).eq.0) then c call linee c call plot(gx1,gy1,3) c end if 30 CONTINUE c**** 200 continue if(lin.eq.1) go to 100 r=0.0 g=0.2 b=1.0 go to 300 210 continue if(lin.eq.1) go to 100 r=0.0 g=1.0 b=0.0 300 continue c iie1=ie(i,j,jj) if(iie1.le.0) go to 100 if(iie1.ge.8) go to 100 c if(iie1.le.1) go to 100 c if(iie1.ge.7) go to 100 c lasl=lin c do 42 kkk=1,2 c call linebe(r,g,b) c do 40 ii=1,lasl lin=ii li1=3*(lin-1)+1 li2=3*(lin-1)+2 li3=3*lin if(kkk.eq.1) then X1=po(li1) Y1=po(li2) Z1=po(li3) else X1=po(li1) Y1=-po(li2) Z1=-po(li3) endif c call linep(x1,y1,z1,ii,lasl) c 40 continue c call lineen(lasl) 42 continue c 100 CONTINUE RETURN END c --------------------- SUBROUTINE AINTE21(IA,AA,F,P) c --------------------- PARAMETER(nline=6000,nlt=nline*3) DIMENSION IA(1),AA(1),F(1),P(1),Q(16),po(nlt,2),ipo(3,2) C IA 1 2 3 4 5 6 7 8 9 10 C NX NY MX MY MZ NXP MI MO C AA 1 2 3 4 5 6 7 8 9 10 C TH0 ARU AR1 ARB XXL YYL ZZL B0 DL SIGN C 11 12 13 14 15 16 17 18 19 20 C DIRE BOUN GX0 GY0 GXL GYL GTH GXMI GXMA EP1 C P(1)=P(1) PI=3.1415926 NX=IA(1) NY=IA(2) MX=IA(3) MY=IA(4) MZ=IA(5) NXP=IA(6) MI=IA(7) MO=IA(8) TH0=AA(1)*PI/180.0 ARU=AA(2) AR1=AA(3) ARB=AA(4) XXL=AA(5) YYL=AA(6) ZZL=AA(7) B0=AA(8) GX0=AA(13) GY0=AA(14) GXL=AA(15) GYL=AA(16) GTH=AA(17)*PI/180.0 GXMI=AA(18) GXMA=AA(19) EP1=AA(20) XLIM=AA(23) XLI2=AA(24) C XX=999.0 XM=0.075 YM=0.075 X1=GX0 Y1=GY0-2.0*YM C CALL CHARA1(X1,Y1,0,0,-1,0,MI,GXMI,GXMA) C XC1=COS(GTH) XS1=SIN(GTH) Y=GY0+0.5*GYL c CALL PLOT(GX0,Y,3) X=GX0 Y=GY0+GYL c CALL PLOT(X,Y,2) X=GX0+GXL c CALL PLOT(X,Y,2) Y=GY0+0.5*GYL c CALL PLOT(X,Y,2) X=GX0 c CALL PLOT(X,Y,2) X=GX0+0.5*GYL*XC1 Y=GY0+0.5*GYL*XS1+0.5*GYL c CALL PLOT(X,Y,2) X=GX0+0.5*GYL*XC1+GXL c CALL PLOT(X,Y,2) X=GX0+GXL Y=GY0+0.5*GYL c CALL PLOT(X,Y,2) C C y=gy0 c call plot(x,y,2) X=GX0 c CALL PLOT(X,Y,2) Y=GY0+0.5*GYL c CALL PLOT(X,Y,2) c XDD1=GXL*(0.0-GXMI)/(GXMA-GXMI) X=GX0+XDD1 Y=GY0+GYL c CALL PLOT(X,Y,3) Y=GY0+0.5*GYL c CALL PLOT(X,Y,2) X=X+0.5*GYL*XC1 Y=Y+0.5*GYL*XS1 c CALL PLOT(X,Y,2) c X=GX0+xdd1 Y=GY0+GYL*0.5 c CALL PLOT(X,Y,2) Y=GY0 c CALL PLOT(X,Y,2) c c call linee C XL0=COS(TH0) XL01=XL0*1.0 XL=2.0*XL0 YL=XL NX1=NX+1 NX2=NX+2 NY1=NY+1 NY2=NY+2 N1=NX2 N2=N1*NY2 AR2=AR1*AR1 DX=XXL/FLOAT(NX1) DY=YYL/FLOAT(NY1) HX=XXL/FLOAT(MX+1) HY=YYL/FLOAT(MY+1) HZ=ZZL/FLOAT(MZ+1) XLMI=0.5*HX*FLOAT(2*NXP-MX-1) XLMA=0.5*HX*FLOAT(2*NXP+MX+1) DL0=HX*AA(9) DL=DL0 LAST=IFIX(XXL/DL) GYMI=-0.1*HY GYMA=AA(21) gymi=-gyma AMPY=0.5*(GXMA-GXMI)*GYL/(GYMA*GXL) c GZMI=-0.1*HZ GZMI=-GYMA GZMA=GYMA AMP=GXL/(GXMA-GXMI) M1=MI M2=M1+1 M3=M2+1 M4=MI+8 M5=M4+1 M6=M5+1 c arcc0=hx*0.5 ardr0=hx*5.0 C JB1=NY2/2+1 c do 100 jj=1,2 DO 100 J=JB1,NY2 C DO 100 J=1,NY1 c DO 100 J=1,NY2 DO 100 I=1,NX2 do 300 jj=1,2 X1=0.5*DX*(2*I-NX2-1)+HX*FLOAT(NXP) Y1=0.5*DY*(2*J-NY2-1) C Y1=0.5*DY*(J) C y1=0.5*(2*J-1) Z1=0.0 c c rlim=sqrt(x1*x1+y1*y1) if(rlim.le.8.0) go to 100 c xbou=x1-(-100.0+4.0*y1*y1)/30.0 x1po=8.0 y1po=10.0 y1po=5.0 y1po=10.0 xbou=x1/x1po+1.0-(y1/y1po)**2 if(xbou.le.0.0) go to 100 c c CALL QUANT1(X1,Y1,Z1,HX,HY,HZ,1,MX,MY,MZ,NXP,AA,F,Q) B1=SQRT(Q(M1)*Q(M1)+Q(M2)*Q(M2)+Q(M3)*Q(M3)) IF(B1.LT.1.0E-20) GO TO 300 AA(11)=1.0 X2=Q(M3) IF(X2.LT.0.0) AA(11)=-1.0 if(jj.eq.2) aa(11)=-aa(11) C X2=0.0 Y2=0.0 Z2=0.0 X3=0.0 Y3=0.0 Z3=0.0 GX1=GX0+XDD1+AMP*(X1+AMPY*Y1*XC1) GY1=GY0+0.5*GYL+AMP*AMPY*(Z1+Y1*XS1) c CALL PLOT(GX1,GY1,3) lin=1 li1=3*(lin-1)+1 li2=3*(lin-1)+2 li3=3*lin po(li1,jj)=X1 po(li2,jj)=Y1 po(li3,jj)=Z1 DX2=0.0 DY2=0.0 DZ2=0.0 DR2=SQRT(DX2*DX2+DY2*DY2+DZ2*DZ2) DR2=AMAX1(dr2,1.0e-8) arcc=0.0 drar=0.0 c**** DO 30 III=1,LAST R=SQRT((X3-X2)**2+(Y3-Y2)**2+(Z3-Z2)**2) DL=R/EP1 IF(DL.GT.0.5*HX) DL=0.5*HX IF(DL.LT.DL0) DL=DL0 DL1=DL*AA(11) X2=X1+DL1*Q(M1)/B1 Y2=Y1+DL1*Q(M2)/B1 Z2=Z1+DL1*Q(M3)/B1 R2=SQRT(X2*X2+Y2*Y2+Z2*Z2) IF(R2.LT.0.5*ARB) GO TO 210 c XLIM=0.0 IF(X2.LT.XLIM.AND.Y2.LT.GYMI) GO TO 100 IF(X2.LT.XLIM.AND.Y2.GT.GYMA) GO TO 100 IF(X2.LT.XLIM.AND.Z2.LT.GZMI) GO TO 100 IF(X2.LT.XLIM.AND.Z2.GT.GZMA) GO TO 100 IF(X2.LT.XLI2.AND.Y2.LT.GYMI) GO TO 220 IF(X2.LT.XLI2.AND.Y2.GT.GYMA) GO TO 220 IF(X2.LT.XLI2.AND.Z2.LT.GZMI) GO TO 220 IF(X2.LT.XLI2.AND.Z2.GT.GZMA) GO TO 220 IF(X2.LT.GXMI.OR.X2.GT.GXMA) GO TO 200 IF(Y2.LT.GYMI.OR.Y2.GT.GYMA) GO TO 200 IF(Z2.LT.GZMI.OR.Z2.GT.GZMA) GO TO 200 CALL QUANT1(X2,Y2,Z2,HX,HY,HZ,2,MX,MY,MZ,NXP,AA,F,Q) B2=SQRT(Q(M4)*Q(M4)+Q(M5)*Q(M5)+Q(M6)*Q(M6)) IF(B2.LT.1.0E-20) GO TO 300 DO 32 II=1,1 DX3=0.5*DL1*(Q(M1)/B1+Q(M4)/B2) DY3=0.5*DL1*(Q(M2)/B1+Q(M5)/B2) DZ3=0.5*DL1*(Q(M3)/B1+Q(M6)/B2) DR3=SQRT(DX3*DX3+DY3*DY3+DZ3*DZ3) X3=X1+DX3 Y3=Y1+DY3 Z3=Z1+DZ3 R2=SQRT(X3*X3+Y3*Y3+Z3*Z3) IF(R2.LT.0.5*ARB) GO TO 210 IF(X3.LT.XLIM.AND.Y3.LT.GYMI) GO TO 100 IF(X3.LT.XLIM.AND.Y3.GT.GYMA) GO TO 100 IF(X3.LT.XLIM.AND.Z3.LT.GZMI) GO TO 100 IF(X3.LT.XLIM.AND.Z3.GT.GZMA) GO TO 100 IF(X3.LT.XLI2.AND.Y3.LT.GYMI) GO TO 220 IF(X3.LT.XLI2.AND.Y3.GT.GYMA) GO TO 220 IF(X3.LT.XLI2.AND.Z3.LT.GZMI) GO TO 220 IF(X3.LT.XLI2.AND.Z3.GT.GZMA) GO TO 220 IF(X3.LT.GXMI.OR.X3.GT.GXMA) GO TO 200 IF(Y3.LT.GYMI.OR.Y3.GT.GYMA) GO TO 200 IF(Z3.LT.GZMI.OR.Z3.GT.GZMA) GO TO 200 CALL QUANT1(X3,Y3,Z3,HX,HY,HZ,2,MX,MY,MZ,NXP,AA,F,Q) B2=SQRT(Q(M4)*Q(M4)+Q(M5)*Q(M5)+Q(M6)*Q(M6)) IF(B2.LT.1.0E-20) GO TO 300 32 CONTINUE C GX1=GX0+XDD1+AMP*(X3+AMPY*Y3*XC1) GY1=GY0+0.5*GYL+AMP*AMPY*(Z3+Y3*XS1) c CALL PLOT(GX1,GY1,2) c arc1=sqrt((dx2/dr2-dx3/dr3)**2+(dy2/dr2-dy3/dr3)**2+ 1 (dz2/dr2-dz3/dr3)**2) arcc=arcc+arc1 ardr=ardr+dr3 DX2=DX3 DY2=DY3 DZ2=DZ3 DR2=DR3 if(arcc.gt.arcc0.or.ardr.gt.ardr0) go to 326 go to 327 326 continue c arcc=0.0 ardr=0.0 lin=lin+1 c lin=iii+1 li1=3*(lin-1)+1 li2=3*(lin-1)+2 li3=3*lin po(li1,jj)=X3 po(li2,jj)=Y3 po(li3,jj)=Z3 327 continue X1=X3 Y1=Y3 Z1=Z3 B1=B2 DO 34 II=1,8 Q(II)=Q(II+8) 34 CONTINUE c if(mod(iii,300).eq.0) then c call linee c call plot(gx1,gy1,3) c end if 30 CONTINUE c go to 100 c**** 200 continue ipo(1,jj)=lin ipo(2,jj)=1 c if(lin.eq.1) go to 100 c call newrgb(1.0,0.0,0.0) r=1.0 g=0.0 b=0.0 go to 300 220 continue ipo(1,jj)=lin ipo(2,jj)=2 c if(lin.eq.1) go to 100 c call newrgb(1.0,0.8,0.0) r=1.0 g=0.5 b=0.0 go to 300 210 continue ipo(1,jj)=lin ipo(2,jj)=3 c if(lin.eq.1) go to 100 c call newrgb(0.0,1.0,0.0) r=0.0 g=1.0 b=0.0 300 continue c line1=ipo(1,1) line2=ipo(1,2) col1=ipo(2,1) col2=ipo(2,2) line1=min(line1,line2) if(line1.eq.1) go to 100 if(col1.eq.1.and.col2.eq.1) color=1 if(col1.eq.1.and.col2.eq.2) color=2 if(col1.eq.1.and.col2.eq.3) color=4 if(col1.eq.2.and.col2.eq.1) color=2 if(col1.eq.2.and.col2.eq.2) color=2 if(col1.eq.2.and.col2.eq.3) color=4 if(col1.eq.3.and.col2.eq.1) color=4 if(col1.eq.3.and.col2.eq.2) color=4 if(col1.eq.3.and.col2.eq.3) color=3 ipo(2,1)=color ipo(2,2)=color c do 44 kkk=1,2 do 42 jj=1,2 lin=ipo(1,jj) if(lin.eq.1) go to 42 icol=ipo(2,jj) c if(icol.eq.1) call newrgb(1.0,0.0,0.0) c if(icol.eq.2) call newrgb(1.0,0.8,0.0) c if(icol.eq.3) call newrgb(0.0,1.0,0.0) c if(icol.eq.4) call newrgb(0.0,0.5,1.0) c r=0.0 g=0.0 b=0.0 if(icol.eq.1) r=1.0 if(icol.eq.2) r=1.0 if(icol.eq.2) g=0.8 if(icol.eq.3) g=1.0 if(icol.eq.4) g=0.2 if(icol.eq.4) b=1.0 cccc if(icol.ge.2.0) go to 42 cccc c lasl=lin c call linebe(r,g,b) c do 40 ii=1,lasl lin=ii li1=3*(lin-1)+1 li2=3*(lin-1)+2 li3=3*lin if(kkk.eq.1) then X1=po(li1,jj) Y1=po(li2,jj) Z1=po(li3,jj) else X1=po(li1,jj) Y1=-po(li2,jj) Z1=-po(li3,jj) endif c call linep(x1,y1,z1,ii,lasl) c 40 continue c call lineen(lasl) c 42 continue 44 continue 100 CONTINUE RETURN END c --------------------- c subroutine symblv2 c --------------------- subroutine symblv2(r,g,b,x,y,z,h,chr,n) character chr*80,ica*80,ich(80)*1 equivalence (ica,ich) ica=chr write(10,10) c 10 format("Viewpoint { position 10 4 8 }"/"Transform {") 10 format("Transform {") write(10,12) x,y,z 12 format("translation ",3(1x,f5.1)) write(10,14) c 14 format("children ["/"Shape {"/"geometry Text {") 14 format("children ["/"Shape {"/"appearance Appearance {"/ & "material Material {") write(10,141) r,g,b 141 format("diffuseColor",3(1x,f3.1)/"}}"/"geometry Text {") c write(10,16) (ich(i),i=1,n) c 16 format("string ",'"',80a1,'"') write(10,*) "string ",'"',(ich(i),i=1,n),'"' write(10,18) 18 format("fontStyle FontStyle {") write(10,20) h 20 format("size ",f5.1) write(10,22) 22 format('family "SERIF"'/'style "ROMAN"') write(10,24) 24 format("}"/"}"/"}"/"]"/"}") return end c --------------------- c subroutine sphere1i c --------------------- subroutine sphere1(x,y,z,r,ir,ig,ib) write(10,10) 10 format("Transform {") write(10,12) x,y,z 12 format("translation ",3(1x,f5.1)) write(10,14) 14 format("children ["/"Shape{") write(10,16) r 16 format("geometry Sphere {radius",f5.1,"}") write(10,18) 18 format("appearance Appearance {") write(10,20) ir,ig,ib 20 format("material Material {diffuseColor",3i2,"} }") write(10,22) 22 format("}]}") return end c --------------------- c subroutine image3 c --------------------- subroutine image3(nx,ny,xb,yb,xl,yl,ico,icc,zcc,v,u,p) dimension p(1),u(1),v(1) ny0=ny-1 nx2=nx*2 nxy=nx*ny c do 10 i=1,nxy if(u(i).lt.0.0) u(i)=0.0 if(u(i).gt.1.0) u(i)=1.0 c u(i)=amax(u(i),0.0) c u(i)=amin(u(i),1.0) 10 continue c do 100 ii=1,ny0 c call trans2(nx,ny,ii,u,p) call facebe c call colorp(nx2,p,ico) c call coordbe c c call coordp(nx2,p) call coordp(ii,nx,ny,xb,yb,xl,yl,icc,zcc,p) c call coorden(nx) call faceen c 100 continue c return end c --------------------- c subroutine colorp5 c --------------------- subroutine colorp5(nx,p,ico) dimension p(1) r=1.0 g=1.0 b=1.0 do 100 i=1,nx x1=p(i) c c a0=1.20 if(ico.eq.2) go to 32 if(ico.eq.3) go to 33 if(ico.eq.4) go to 34 if(ico.eq.5) go to 35 if(ico.eq.6) go to 36 if(ico.eq.7) go to 37 c x1=2.0*x1 r=x1 g=2.0-x1 b=0.0 go to 40 32 continue r=x1 g=x1 b=x1 go to 40 33 continue x1=7.0*x1 b=x1 if(x1.ge.2.0) b=3.0-x1 if(x1.ge.5.0) b=0.5*(x1-5.0) g=x1-1.0 if(x1.ge.4.0) g=5.0-x1 if(x1.ge.5.0) g=0.5*(x1-5.0) r=x1-3.0 go to 40 34 continue x1=10.0*x1 b=x1 if(x1.ge.3.0) b=4.0-x1 if(x1.ge.6.0) b=0.25*(x1-6.0) g=x1-2.0 if(x1.ge.5.0) g=6.0-x1 if(x1.ge.6.0) g=0.25*(x1-6.0) r=x1 if(x1.ge.1.0) r=2.0-x1 if(x1.ge.4.0) r=x1-4.0 go to 40 35 continue a0=2.0 y1=x1*(a0+4.0) b=a0-y1 if(y1.ge.a0) b=y1-(a0+1.0) r=a0+3.0-y1 g=a0+1.0-y1 if(y1.ge.a0+2.0) g=y1-(a0+3.0) go to 40 36 continue a0=2.0 y1=x1*(a0+4.0) r=a0-y1 if(y1.ge.a0) r=y1-(a0+1.0) b=a0+3.0-y1 g=a0+1.0-y1 if(y1.ge.a0+2.0) g=y1-(a0+3.0) go to 40 37 continue x1=7.0*x1 b=x1 b=1.0 if(x1.ge.2.0) b=3.0-x1 if(x1.ge.5.0) b=0.5*(x1-5.0) g=1.0-x1 if(x1.ge.1.0) g=x1-1.0 if(x1.ge.4.0) g=5.0-x1 if(x1.ge.5.0) g=0.5*(x1-5.0) r=1.0-x1 if(x1.ge.1.0) r=x1-3.0 go to 40 c 40 continue c r=amax1(r,0.0) g=amax1(g,0.0) b=amax1(b,0.0) r=amin1(r,1.0) g=amin1(g,1.0) b=amin1(b,1.0) c if(i.ne.nx) write(10,20) r,g,b if(i.eq.nx) write(10,30) r,g,b 20 format(f4.2,1x,f4.2,1x,f4.2,",") 30 format(f4.2,1x,f4.2,1x,f4.2) 100 continue return end c --------------------- c subroutine coordp5 c --------------------- subroutine coordp5(ii,nx,ny,xb,yb,xl,yl,icc,zcc,p) dimension p(1) x=1.0 y=1.0 z=1.0 c xb=0.0 c yb=0.0 c dx=1.0 c dy=1.0 nxy=nx*ny dx=xl/float(nx-1) dy=yl/float(ny-1) c do 100 j1=1,2 do 100 i=1,nx do 100 j1=1,2 j=j1-1+ii ixy=i+nx*(j-1) c if (icc.eq.2) go to 52 if (icc.eq.3) go to 53 c x=xb+dx*float(i-1) y=yb+dy*float(j-1) z=zcc go to 60 52 continue c y=xb+dx*float(i-1) c z=yb+dy*float(j-1) z=xb+dx*float(i-1) y=yb+dy*float(j-1) x=zcc go to 60 53 continue x=xb+dx*float(i-1) z=yb+dy*float(j-1) y=zcc 60 continue c if(i.ne.nx) write(10,20) x,y,z if(i.eq.nx) write(10,30) x,y,z 20 format(f6.2,1x,f6.2,1x,f6.2,",") 30 format(f6.2,1x,f6.2,1x,f6.2) 100 continue return end c --------------------- c subroutine coordt1 c --------------------- subroutine coordt1(ia,aa,f,u) dimension ia(1),aa(1),f(1),u(1) c nx=ia(3) ny=ia(4) nz=ia(5) nxp=ia(6) nux=ia(9) nuy=ia(10) icc=ia(11) th0=aa(1) aru=aa(2) ar1=aa(3) arb=aa(4) xxl=aa(5) yyl=aa(6) zzl=aa(7) zcc=aa(25) c nx1=nx+1 ny1=ny+1 nz1=nz+1 nx2=nx+2 ny2=ny+2 nz2=nz+2 n1=nx2 n2=n1*ny2 n3=n2*nz2 hx=xxl/float(nx1) hy=yyl/float(ny1) hz=zzl/float(nz1) xb=aa(18) yb=0.0 zb=0.0 xl1=aa(19) yl1=aa(21) zl1=aa(21) xl=xl1-xb yl=yl1-yb zl=zl1-zb c if(icc.eq.2) go to 20 if(icc.eq.3) go to 30 c c xz plane y=zcc do 12 j=1,nuy do 12 i=1,nux z=zb+zl*float(j-1)/float(nuy-1) x=xb+xl*float(i-1)/float(nux-1) ak1=z/hz+1.5 aj1=y/hy+1.5 ai1=x/hx+0.5*float(1+nx2)-float(nxp) ka1=ak1 ja1=aj1 ia1=ai1 adz=ak1-float(ka1) ady=aj1-float(ja1) adx=ai1-float(ia1) adx1=1.0-adx ady1=1.0-ady adz1=1.0-adz i1=i+nux*(j-1) i2=ia1+n1*(ja1-1)+n2*(ka1-1)+n3*3 u(i1)=adx*ady*adz*f(i2+1+n1+n2)+adx1*ady*adz*f(i2+n1+n2) 1 +adx*ady1*adz*f(i2+1+n2)+adx1*ady1*adz*f(i2+n2) 2 +adx*ady*adz1*f(i2+1+n1)+adx1*ady*adz1*f(i2+n1) 3 +adx*ady1*adz1*f(i2+1)+adx1*ady1*adz1*f(i2) 12 continue go to 40 c c yz plane 20 continue x=zcc do 22 j=1,nuy do 22 i=1,nux z=zb+zl*float(j-1)/float(nuy-1) y=yb+yl*float(i-1)/float(nux-1) ak1=z/hz+1.5 aj1=y/hy+1.5 ai1=x/hx+0.5*float(1+nx2)-float(nxp) ka1=ak1 ja1=aj1 ia1=ai1 adz=ak1-float(ka1) ady=aj1-float(ja1) adx=ai1-float(ia1) adx1=1.0-adx ady1=1.0-ady adz1=1.0-adz i1=i+nux*(j-1) i2=ia1+n1*(ja1-1)+n2*(ka1-1)+n3*3 u(i1)=adx*ady*adz*f(i2+1+n1+n2)+adx1*ady*adz*f(i2+n1+n2) 1 +adx*ady1*adz*f(i2+1+n2)+adx1*ady1*adz*f(i2+n2) 2 +adx*ady*adz1*f(i2+1+n1)+adx1*ady*adz1*f(i2+n1) 3 +adx*ady1*adz1*f(i2+1)+adx1*ady1*adz1*f(i2) 22 continue go to 40 c c xy plane 30 continue z=zcc do 32 j=1,nuy do 32 i=1,nux y=yb+yl*float(j-1)/float(nuy-1) x=xb+xl*float(i-1)/float(nux-1) ak1=z/hz+1.5 aj1=y/hy+1.5 ai1=x/hx+0.5*float(1+nx2)-float(nxp) ka1=ak1 ja1=aj1 ia1=ai1 adz=ak1-float(ka1) ady=aj1-float(ja1) adx=ai1-float(ia1) adx1=1.0-adx ady1=1.0-ady adz1=1.0-adz i1=i+nux*(j-1) i2=ia1+n1*(ja1-1)+n2*(ka1-1)+n3*3 u(i1)=adx*ady*adz*f(i2+1+n1+n2)+adx1*ady*adz*f(i2+n1+n2) 1 +adx*ady1*adz*f(i2+1+n2)+adx1*ady1*adz*f(i2+n2) 2 +adx*ady*adz1*f(i2+1+n1)+adx1*ady*adz1*f(i2+n1) 3 +adx*ady1*adz1*f(i2+1)+adx1*ady1*adz1*f(i2) 32 continue 40 continue c write(6,*) icc,zcc,x,y,z return end c --------------------- c subroutine limit1 c --------------------- subroutine limit1(kk,fac,nux,nuy,vmax,vmin,umax,umin,u) dimension u(1) c if(kk.ge.2) go to 12 umax=u(1) umin=u(1) do 10 j=1,nuy do 10 i=1,nux i1=i+nux*(j-1) if(umax.lt.u(i1)) umax=u(i1) if(umin.gt.u(i1)) umin=u(i1) 10 continue c 12 continue c if(kk.eq.0) go to 30 do 20 j=1,nuy do 20 i=1,nux i1=i+nux*(j-1) if(kk.le.1) u(i1)=(u(i1)-umin)/(umax-umin) if(kk.ge.2) u(i1)=(u(i1)-vmin)/(vmax-vmin) 20 continue 30 continue return end c --------------------- SUBROUTINE AINTE1A(IA,AA,F,P) c --------------------- parameter(nnx=102,nny=102,nnz=2) PARAMETER(nline=9900,nlt=nline*3) DIMENSION IA(1),AA(1),F(1),P(1),Q(16),po(nlt) common app(nnx,nny,nnz),ie(nnx,nny,nnz) C IA 1 2 3 4 5 6 7 8 9 10 C NX NY MX MY MZ NXP MI MO C AA 1 2 3 4 5 6 7 8 9 10 C TH0 ARU AR1 ARB XXL YYL ZZL B0 DL SIGN C 11 12 13 14 15 16 17 18 19 20 C DIRE BOUN GX0 GY0 GXL GYL GTH GXMI GXMA EP1 C do 330 k=1,2 do 330 j=1,102 do 330 i=1,102 APP(I,J,K)=0.0 330 continue c P(1)=P(1) PI=3.1415926 NX=IA(1) NY=IA(2) MX=IA(3) MY=IA(4) MZ=IA(5) NXP=IA(6) MI=IA(7) MO=IA(8) TH0=AA(1)*PI/180.0 ARU=AA(2) AR1=AA(3) ARB=AA(4) XXL=AA(5) YYL=AA(6) ZZL=AA(7) B0=AA(8) GX0=AA(13) GY0=AA(14) GXL=AA(15) GYL=AA(16) GTH=AA(17)*PI/180.0 GXMI=AA(18) GXMA=AA(19) EP1=AA(20) ccc waku kaku tokoro C XX=999.0 XM=0.075 YM=0.075 X1=GX0 Y1=GY0-2.0*YM C CALL CHARA1(X1,Y1,0,0,-1,0,MI,GXMI,GXMA) C XC1=COS(GTH) XS1=SIN(GTH) Y=GY0+0.5*GYL c CALL PLOT(GX0,Y,3) X=GX0 Y=GY0+GYL c CALL PLOT(X,Y,2) X=GX0+GXL c CALL PLOT(X,Y,2) Y=GY0+0.5*GYL c CALL PLOT(X,Y,2) X=GX0 c CALL PLOT(X,Y,2) X=GX0+0.5*GYL*XC1 Y=GY0+0.5*GYL*XS1+0.5*GYL c CALL PLOT(X,Y,2) X=GX0+0.5*GYL*XC1+GXL c CALL PLOT(X,Y,2) X=GX0+GXL Y=GY0+0.5*GYL c CALL PLOT(X,Y,2) C y=gy0 c call plot(x,y,2) X=GX0 c CALL PLOT(X,Y,2) Y=GY0+0.5*GYL c CALL PLOT(X,Y,2) c XDD1=GXL*(0.0-GXMI)/(GXMA-GXMI) X=GX0+XDD1 Y=GY0+GYL c CALL PLOT(X,Y,3) Y=GY0+0.5*GYL c CALL PLOT(X,Y,2) X=X+0.5*GYL*XC1 Y=Y+0.5*GYL*XS1 c CALL PLOT(X,Y,2) c X=GX0+xdd1 Y=GY0+GYL*0.5 c CALL PLOT(X,Y,2) Y=GY0 c CALL PLOT(X,Y,2) c c call linee ccc waku kakiowatta! C ccc boundary?? condition ?? XL0=COS(TH0) XL01=XL0*1.0 XL=2.0*XL0 YL=XL NX1=NX+1 NX2=NX+2 NY1=NY+1 NY2=NY+2 N1=NX2 N2=N1*NY2 AR2=AR1*AR1 DX=XL/FLOAT(NX1) DY=YL/FLOAT(NY1) HX=XXL/FLOAT(MX+1) HY=YYL/FLOAT(MY+1) HZ=ZZL/FLOAT(MZ+1) XLMI=0.5*HX*FLOAT(2*NXP-MX-1) XLMA=0.5*HX*FLOAT(2*NXP+MX+1) DL0=HX*AA(9) DL=DL0 LAST=IFIX(XXL/DL) GYMI=-0.1*HY GYMA=AA(21) AMPY=0.5*(GXMA-GXMI)*GYL/(GYMA*GXL) c GZMI=-0.1*HZ gymi=-gyma gzmi=-gyma c GZMI GZMA ?? GZMA=GYMA AMP=GXL/(GXMA-GXMI) M1=MI M2=M1+1 M3=M2+1 M4=MI+8 M5=M4+1 M6=M5+1 c arcc0=hx*0.5 ardr0=hx*5.0 C cc polar cap kara start suru basyo wo kimeru JB1=NY2/2+1 do 100 jj=1,2 DO 100 J=JB1,NY2 cc DO 100 J=1,NY2 DO 100 I=1,NX2 X=0.5*DX*(2*I-NX2-1) Y=0.5*DY*(2*J-NY2-1) R=SQRT(X*X+Y*Y) IF(R.GT.XL01) GO TO 100 c z ga + sika naiyouna.... Z=SQRT(1.0-R*R) R=ARB*SQRT(ARB) X1=X*R Y1=Y*R Z1=ARB*ARB-X1*X1-Y1*Y1 c zga 0.0 yori tiisaitoki 100 he IF(Z1.LE.0.0) GO TO 100 c IF(Z1.LE.0.0) GO TO 1030 if(jj.eq.1) Z1=SQRT(Z1) if(jj.eq.2) Z1=-SQRT(Z1) 1040 CALL QUANT1(X1,Y1,Z1,HX,HY,HZ,1,MX,MY,MZ,NXP,AA,F,Q) B1=SQRT(Q(M1)*Q(M1)+Q(M2)*Q(M2)+Q(M3)*Q(M3)) IF(B1.LT.1.0E-20) GO TO 100 AA(11)=1.0 X2=Q(M1)*X1+Q(M2)*Y1+Q(M3)*Z1 IF(X2.LT.0.0) AA(11)=-1.0 C X2=0.0 Y2=0.0 Z2=0.0 X3=0.0 Y3=0.0 Z3=0.0 GX1=GX0+XDD1+AMP*(X1+AMPY*Y1*XC1) cc GY1=GY0+0.5*GYL+AMP*AMPY*(Z1+Y1*XS1) GY1=GY0+0.5*GYL+AMP*AMPY*(Z1+Y1*XS1) c CALL PLOT(GX1,GY1,3) lin=1 li1=3*(lin-1)+1 li2=3*(lin-1)+2 li3=3*lin po(li1)=X1 po(li2)=Y1 po(li3)=Z1 DX2=0.0 DY2=0.0 DZ2=0.0 DR2=SQRT(DX2*DX2+DY2*DY2+DZ2*DZ2) DR2=AMAX1(dr2,1.0e-8) arcc=0.0 drar=0.0 c**** DO 30 III=1,LAST R=SQRT((X3-X2)**2+(Y3-Y2)**2+(Z3-Z2)**2) DL=R/EP1 IF(DL.GT.0.5*HX) DL=0.5*HX IF(DL.LT.DL0) DL=DL0 DL1=DL*AA(11) X2=X1+DL1*Q(M1)/B1 Y2=Y1+DL1*Q(M2)/B1 Z2=Z1+DL1*Q(M3)/B1 R2=SQRT(X2*X2+Y2*Y2+Z2*Z2) IF(R2.LT.0.5*ARB) GO TO 210 IF(X2.LT.GXMI.OR.X2.GT.GXMA) GO TO 200 IF(Y2.LT.GYMI.OR.Y2.GT.GYMA) GO TO 200 IF(Z2.LT.GZMI.OR.Z2.GT.GZMA) GO TO 200 CALL QUANT1(X2,Y2,Z2,HX,HY,HZ,2,MX,MY,MZ,NXP,AA,F,Q) B2=SQRT(Q(M4)*Q(M4)+Q(M5)*Q(M5)+Q(M6)*Q(M6)) IF(B2.LT.1.0E-20) GO TO 100 DO 32 II=1,1 DX3=0.5*DL1*(Q(M1)/B1+Q(M4)/B2) DY3=0.5*DL1*(Q(M2)/B1+Q(M5)/B2) DZ3=0.5*DL1*(Q(M3)/B1+Q(M6)/B2) DR3=SQRT(DX3*DX3+DY3*DY3+DZ3*DZ3) X3=X1+DX3 Y3=Y1+DY3 Z3=Z1+DZ3 R2=SQRT(X3*X3+Y3*Y3+Z3*Z3) IF(R2.LT.0.5*ARB) GO TO 210 IF(X3.LT.GXMI.OR.X3.GT.GXMA) GO TO 200 IF(Y3.LT.GYMI.OR.Y3.GT.GYMA) GO TO 200 IF(Z3.LT.GZMI.OR.Z3.GT.GZMA) GO TO 200 CALL QUANT1(X3,Y3,Z3,HX,HY,HZ,2,MX,MY,MZ,NXP,AA,F,Q) B2=SQRT(Q(M4)*Q(M4)+Q(M5)*Q(M5)+Q(M6)*Q(M6)) IF(B2.LT.1.0E-20) GO TO 100 32 CONTINUE C GX1=GX0+XDD1+AMP*(X3+AMPY*Y3*XC1) cc GY1=GY0+0.5*GYL+AMP*AMPY*(Z3+Y3*XS1) GY1=GY0+0.5*GYL+AMP*AMPY*(Z3+Y3*XS1) c CALL PLOT(GX1,GY1,2) c arc1=sqrt((dx2/dr2-dx3/dr3)**2+(dy2/dr2-dy3/dr3)**2+ 1 (dz2/dr2-dz3/dr3)**2) arcc=arcc+arc1 ardr=ardr+dr3 DX2=DX3 DY2=DY3 DZ2=DZ3 DR2=DR3 if(arcc.gt.arcc0.or.ardr.gt.ardr0) go to 326 go to 327 326 continue c arcc=0.0 ardr=0.0 lin=lin+1 c lin=iii+1 li1=3*(lin-1)+1 li2=3*(lin-1)+2 li3=3*lin po(li1)=X3 po(li2)=Y3 po(li3)=Z3 327 continue X1=X3 Y1=Y3 Z1=Z3 B1=B2 DO 34 II=1,8 Q(II)=Q(II+8) 34 CONTINUE c if(mod(iii,100).eq.0) then c call linee c call plot(gx1,gy1,3) c end if 30 CONTINUE c**** 200 continue if(lin.eq.1) go to 100 c call newrgb(0.0,0.5,1.0) APP1=4.0 go to 300 210 continue if(lin.eq.1) go to 100 c call newrgb(0.0,1.0,0.0) APP1=1.0 300 continue lasl=lin lin=1 li1=3*(lin-1)+1 li2=3*(lin-1)+2 li3=3*lin X1=po(li1) Y1=po(li2) Z1=po(li3) c C c GX1=GX0+XDD1+AMP*(X1+AMPY*Y1*XC1) cc GY1=GY0+0.5*GYL+AMP*AMPY*(Z1+Y1*XS1) c GY1=GY0+0.5*GYL+AMP*AMPY*(Z1+Y1*XS1) c GX1=po(li1) c GY1=po(li2) c call plot(gx1,gy1,3) c do 40 ii=2,lasl c lin=ii c li1=3*(lin-1)+1 c li2=3*(lin-1)+2 c li3=3*lin c X1=po(li1) c Y1=po(li2) c Z1=po(li3) c C c GX1=GX0+XDD1+AMP*(X1+AMPY*Y1*XC1) cc GY1=GY0+0.5*GYL+AMP*AMPY*(Z1+Y1*XS1) c GY1=GY0+0.5*GYL+AMP*AMPY*(Z1+Y1*XS1) c GX1=po(li1) c GY1=po(li2) c call plot(gx1,gy1,2) c if(mod(ii,100).eq.0) then c call linee c call plot(gx1,gy1,3) c end if c 40 continue c call linee APP(I,J,JJ)=APP1 100 CONTINUE c write(20) app c call newrgb(0.0,0.0,0.0) RETURN END c --------------------- SUBROUTINE zsub33(ia,aa) c --------------------- parameter(nx=102,ny=102,nz=2) c parameter(last=1,nxa=102,nya=102) c parameter(last=1,nxa=102,nya=102) c dimension id(nx,ny,nz),ie(nx,ny,nz) c dimension id(nxa,nya,nz),ie(nxa,nya,nz) dimension ic(nx,ny,nz),id(nx,ny,nz) dimension ia(1),aa(1) common app(nx,ny,nz),ie(nx,ny,nz) c c c nxg=ia(1) nyg=ia(2) nxa=nxg+2 nya=nyg+2 c do 10 k=1,nz do 10 j=1,nya do 10 i=1,nxa ic(i,j,k)=ifix(app(i,j,k)) id(i,j,k)=ic(i,j,k) if(id(i,j,k).le.1) id(i,j,k)=0 if(id(i,j,k).gt.1) id(i,j,k)=8 if(id(i,j,k).le.1) ie(i,j,k)=id(i,j,k) if(id(i,j,k).gt.1) ie(i,j,k)=id(i,j,k) 10 continue c c nyah=nya/2 do 200 k=1,nz do 40 j=1,nya if(k.eq.1.and.j.le.nyah) j1=nya-j+1 if(k.eq.2.and.j.gt.nyah) j1=j do 40 i=1,nxa if(k.eq.1.and.j.le.nyah) id(i,j,1)=ie(i,j1,k) if(k.eq.2.and.j.gt.nyah) id(i,j,1)=ie(i,j1,k) 40 continue 200 continue c do 300 j=2,nya-1 do 300 i=2,nxa-1 id(i,j,2)=(id(i+1,j,1)+id(i-1,j,1)+id(i,j-1,1)+id(i,j+1,1) 1 +4*id(i,j,1))/8 300 continue c do 400 k=1,nz do 50 j=1,nya if(k.eq.1.and.j.le.nyah) j1=nya-j+1 if(k.eq.2.and.j.gt.nyah) j1=j do 50 i=1,nxa if(k.eq.1.and.j.le.nyah) ie(i,j1,k)=id(i,j,2) if(k.eq.2.and.j.gt.nyah) ie(i,j1,k)=id(i,j,2) 50 continue 400 continue c c c write(iotap) ie c c return end c --------------------- c subroutine backgr1 c --------------------- subroutine backgr1(r,g,b) write(10,10) 10 format("Background {", "groundColor [") write(10,12) r,g,b c 12 format(3(1x,f3.1),",") 12 format(f3.1,2(1x,f3.1),",") write(10,12) r,g,b write(10,14) r,g,b c 14 format(3(1x,f3.1)) 14 format(f3.1,2(1x,f3.1)) write(10,16) 16 format("]") write(10,18) 18 format("groundAngle [1.05, 1.57]") write(10,20) 20 format("skyColor [") write(10,12) r,g,b write(10,12) r,g,b write(10,14) r,g,b write(10,16) write(10,22) 22 format("skyAngle [1.05, 1.57]") write(10,24) 24 format("}") return end c c------------------------------ subroutine lat3d1(ic,aa) c------------------------------ dimension ic(1),aa(1) c nx=ic(1) ny=ic(2) nz=ic(3) xl=aa(1) yl=aa(2) zl=aa(3) r=aa(4) g=aa(5) b=aa(6) c lasl=max(nx,ny,nz)+1 las2=lasl*2 iid=las2*3 nxs=nx-1 nys=ny-1 nzs=nz-1 dx=xl/nx dy=yl/ny dz=zl/nz nz1=nz+1 c c ** y >= 0 do 10 k=1,nz z=dz*(k-1) call posetp1(dx,dy,nx,ny,r,g,b,z) 10 continue c ** y < 0 do 20 k=1,nz z=-1*dz*(k-1) call posetp1(dx,dy,nx,ny,r,g,b,z) 20 continue c h=1.0 x1=-1.1*yl y1=0.0 z1=0.0 call symblvc(x1,y1,z1,h,r,g,b,"X 40Re",6) x1=1.02*xl y1=0.0 z1=0.0 call symblvc(x1,y1,z1,h,r,g,b,"-160Re",6) x1=0.0 y1=1.1*yl z1=0.0 call symblvc(x1,y1,z1,h,r,g,b,"Z 40Re",6) x1=0.0 y1=0.0 z1=1.1*zl call symblvc(x1,y1,z1,h,r,g,b,"Y 40Re",6) c return end c c------------------------ c subroutine posetp1 c------------------------ subroutine posetp1(dx,dy,nx,ny,r,g,b,z) c dimension id(9000) nx1=nx+1 ny1=ny+1 c c ==== yokojiku ==== x2=dx*nx x1=-1*x2/4 call linebe(r,g,b) lasl=ny*2 do 10 j=1,ny1 y=dy*(j-1) call linepl(x1,y,z,1,lasl) call linepl(x2,y,z,j*2,lasl) 10 continue c do 11 j=1,ny1 ii=(j-1)*3+1 id(ii)=(j-1)*2 id(ii+1)=id(ii)+1 id(ii+2)=-1 11 continue idd=ny1*3 call coordenl(idd,id) call lineenl c c ==== tatejiku ==== y1=0.0 y2=dy*ny c ** x >= 0 call linebe(r,g,b) lasl=nx*2 do 20 i=1,nx1 x=dx*(i-1) call linepl(x,y1,z,1,lasl) call linepl(x,y2,z,i*2,lasl) 20 continue c do 21 i=1,nx1 ii=(i-1)*3+1 id(ii)=(i-1)*2 id(ii+1)=id(ii)+1 id(ii+2)=-1 21 continue idd=nx1*3 call coordenl(idd,id) call lineenl c ** x < 0 mx=ny call linebe(r,g,b) lasl=mx*2 do 30 i=1,mx x=-1*dx*i call linepl(x,y1,z,1,lasl) call linepl(x,y2,z,i*2,lasl) 30 continue c do 31 i=1,mx ii=(i-1)*3+1 id(ii)=(i-1)*2 id(ii+1)=id(ii)+1 id(ii+2)=-1 31 continue idd=mx*3 call coordenl(idd,id) call lineenl c 100 continue c return end c c c-------------------------- c subroutine coordenl c-------------------------- subroutine coordenl(nx,ib) c dimension ib(1) c write(10,21) 21 format("]") write(10,22) 22 format("}") write(10,23) 23 format("colorIndex [ ") write(10,24) (0,i=1,nx/3) 24 format(20i3) write(10,25) 25 format("]") write(10,26) 26 format("coordIndex [") c write(10,102) (ib(i),i=1,nx) 102 format(10i5) c return end c c -------------------- c subroutine lineenl c -------------------- subroutine lineenl c write(10,20) 20 format("]") write(10,21) 21 format("colorPerVertex FALSE") write(10,22) 22 format("}") write(10,23) 23 format("}") return end c c --------------------- c subroutine viewp c --------------------- subroutine viewp(iviewp) c for VRML plug-in diplay of personal computer c write(10,10) 10 format("Viewpoint {") write(10,20) 20 format("orientation 0 0 -1 0") write(10,30) iviewp 30 format("position 0 0 ",1x,I4) write(10,40) 40 format("}") return end c c --------------------- c subroutine linepl c --------------------- subroutine linepl(x,y,z,ii,lasl) c if(ii.lt.lasl) then write(10,10) x,y,z 10 format(f6.2,1x,f6.2,1x,f6.2,1h,) else write(10,20) x,y,z 20 format(f6.2,1x,f6.2,1x,f6.2) end if return end c------------------------------ subroutine lat3d2(ic,aa) c------------------------------ dimension ic(1),aa(1) c nx=ic(1) ny=ic(2) nz=ic(3) xl=aa(1) yl=aa(2) zl=aa(3) r=aa(4) g=aa(5) b=aa(6) c nxs=nx-1 nys=ny-1 nzs=nz-1 dx=xl/nx dy=yl/ny dz=zl/nz nz1=nz+1 c c ** y >= 0 do 10 k=1,nz z=dz*(k-1) call posetp2(dx,dy,nx,ny,r,g,b,z) 10 continue c ** y < 0 do 20 k=1,nz z=-1*dz*(k-1) call posetp2(dx,dy,nx,ny,r,g,b,z) 20 continue c c h=1.0 x1=-1.3*yl y1=0.0 z1=0.0 call symblvc(x1,y1,z1,h,r,g,b,"X 10Re",6) x1=1.1*xl y1=0.0 z1=0.0 call symblvc(x1,y1,z1,h,r,g,b,"-10Re",5) x1=0.0 y1=1.1*yl z1=0.0 call symblvc(x1,y1,z1,h,r,g,b,"Z 10Re",6) x1=0.0 y1=0.0 z1=1.1*zl call symblvc(x1,y1,z1,h,r,g,b,"Y 10Re",6) c return end c c------------------------ c subroutine posetp2 c------------------------ subroutine posetp2(dx,dy,nx,ny,r,g,b,z) c dimension id(9000) nx1=nx+1 ny1=ny+1 c c ==== yokojiku ==== x2=dx*nx x1=-1*x2 call linebe(r,g,b) lasl=ny*2 do 10 j=1,ny1 y=dy*(j-1) call linepl(x1,y,z,1,lasl) call linepl(x2,y,z,j*2,lasl) 10 continue c do 11 j=1,ny1 ii=(j-1)*3+1 id(ii)=(j-1)*2 id(ii+1)=id(ii)+1 id(ii+2)=-1 11 continue idd=ny1*3 call coordenl(idd,id) call lineenl c c ==== tatejiku ==== y1=0.0 y2=dy*ny c ** x >= 0 call linebe(r,g,b) lasl=nx*2 do 20 i=1,nx1 x=dx*(i-1) call linepl(x,y1,z,1,lasl) call linepl(x,y2,z,i*2,lasl) 20 continue c do 21 i=1,nx1 ii=(i-1)*3+1 id(ii)=(i-1)*2 id(ii+1)=id(ii)+1 id(ii+2)=-1 21 continue idd=nx1*3 call coordenl(idd,id) call lineenl c ** x < 0 mx=ny call linebe(r,g,b) lasl=mx*2 do 30 i=1,mx x=-1*dx*i call linepl(x,y1,z,1,lasl) call linepl(x,y2,z,i*2,lasl) 30 continue c do 31 i=1,mx ii=(i-1)*3+1 id(ii)=(i-1)*2 id(ii+1)=id(ii)+1 id(ii+2)=-1 31 continue idd=mx*3 call coordenl(idd,id) call lineenl c 100 continue c return end c------------------------------ subroutine lat3d3(ic,aa) c------------------------------ parameter (ang=45,lasl=90-ang+1,las2=lasl*2,iid=las2*3) dimension ic(1),aa(1),id(iid) c pi=3.141592 re1=aa(1) re2=aa(2) r=aa(3) g=aa(4) b=aa(5) c === tatejiku === do 40 j=45,90 th1=float(j)*pi/180.0 z1=re1*cos(th1) z2=re2*cos(th1) call linebe(r,g,b) c* z >= 0 do 10 i=45,90 th2=float(i)*pi/180.0 sin1=sin(th2)*sin(th1) cos1=cos(th2) call posetp3(re1,re2,lasl,sin1,cos1,z1,z2) 10 continue c do 12 i=1,las2 ii=(i-1)*3+1 id(ii)=(i-1)*2 id(ii+1)=id(ii)+1 id(ii+2)=-1 12 continue call coordenl(iid,id) call lineenl c c* z < 0 z1=-1*z1 z2=-1*z2 call linebe(r,g,b) do 30 i=45,90 th2=float(i)*pi/180.0 sin1=sin(th2)*sin(th1) cos1=cos(th2) call posetp3(re1,re2,lasl,sin1,cos1,z1,z2) 30 continue c do 32 i=1,las2 ii=(i-1)*3+1 id(ii)=(i-1)*2 id(ii+1)=id(ii)+1 id(ii+2)=-1 32 continue call coordenl(iid,id) call lineenl c 40 continue c c----------- yokojiku --------------- res=(re2-re1)/(lasl-1) c* z >= 0 do 70 k=45,90 th1=float(k)*pi/180.0 do 51 j=45,90 re=re1+res*(j-45) call linebe(r,g,b) do 50 i=45,90 th2=float(i)*pi/180.0 x1=re*cos(th2) y1=re*sin(th2)*sin(th1) z1=re*cos(th1) call linepl(x1,y1,z1,i,lasl) 50 continue call lineen(lasl) 51 continue c do 61 j=45,90 re=re1+res*(j-45) call linebe(r,g,b) do 60 i=45,90 th2=float(i)*pi/180.0 x1=-1*re*cos(th2) y1=re*sin(th2)*sin(th1) z1=re*cos(th1) call linepl(x1,y1,z1,i,lasl) 60 continue call lineen(lasl) 61 continue c 70 continue c c* z < 0 do 100 k=45,90 th1=float(k)*pi/180.0 do 81 j=45,90 re=re1+res*(j-45) call linebe(r,g,b) do 80 i=45,90 th2=float(i)*pi/180.0 x1=re*cos(th2) y1=re*sin(th2)*sin(th1) z1=-1*re*cos(th1) call linepl(x1,y1,z1,i,lasl) 80 continue call lineen(lasl) 81 continue c do 91 j=45,90 re=re1+res*(j-45) call linebe(r,g,b) do 90 i=45,90 th2=float(i)*pi/180.0 x1=-1*re*cos(th2) y1=re*sin(th2)*sin(th2) z1=-1*re*cos(th2) call linepl(x1,y1,z1,i,lasl) 90 continue call lineen(lasl) 91 continue c 100 continue c c return end c c ------------------ c subroutine posetp3 c ------------------ subroutine posetp3(re1,re2,lasl,sin1,cos1,z1,z2) c x1=re1*cos1 y1=re1*sin1 call linepl(x1,y1,z1,1,lasl) x2=re2*cos1 y2=re2*sin1 call linepl(x2,y2,z2,2,lasl) c x1=-1*x1 call linepl(x1,y1,z1,1,lasl) x2=-1*x2 call linepl(x2,y2,z2,2,lasl) c return end