C
C Program SRFSECT to calculate the crossection of a triangulated surface C with a parallelogram. Each of the three sides of each triangle of the C triangulated surface is tested whether it has a crossection with the C specified parallelogram. The crossection points, if found, are written C to the output file. C C Version: 7.50 C Date: 2018, December 11 C C Coded by: Petr Bulant C Department of Geophysics, Charles University Prague, C Ke Karlovu 3, 121 16 Praha 2, Czech Republic, C http://sw3d.cz/staff/bulant.htm C C....................................................................... C C Description of data files: C C Input data read from the standard input device (*): C The data are read by the list directed input (free format) and C consist of a single string 'SEP': C 'SEP'...String in apostrophes containing the name of the input C SEP parameter or history file with the input data. C No default, 'SEP' must be specified and cannot be blank. C C C Input data file 'SEP': C File 'SEP' has the form of the SEP C parameter file. The parameters, which do not differ from their C defaults, need not be specified in file 'SEP'. C Data specifying input files with the triangulated surface: C VRTX='string'... Name of the file with vertices of the polygons. C Description of file VRTX C Default: VRTX='vrtx.out' C TRGL='string'... Name of the file describing the triangles or C polygons. C Description of file TRGL C Default: TRGL='trgl.out' C Data specifying the intersection parallelogram: C A1=real,A2=real,A3=real,B1=real,B2=real,B3=real,C1=real,C2=real, C C3=real... Cartesian coordinates of three points A,B and C C defining the intersection plane. The intersection of the C given triangulated surface with the parallelogram given C by points A,B and C is searched for. C Defaults: A1=0.,A2=0.,A3=0.,B1=0.,B2=0.,B3=0., C C1=0.,C2=0.,C3=0. C Note that the points A,B and C must be different and thus C at least two of the values A1,A2,A3,B1,B2,B3,C1,C2,C3 must C be different from the default values. C Output file: C SRFSECT='string'... Name of the file to store the calculated C crossection of the surface with the parallelogram. C Each of the three sides of each triangle from file TRGL C is tested whether it has a crossection with the specified C parallelogram. The crossection points, if found, are C written to this output file. C Description of file SRFSECT. C Default: SRFSECT='srfsect.out' C Optional parameters specifying the form of the real quantities C written in the output formatted files: C MINDIG,MAXDIG=positive integers ... See the description in file C forms.for. C C C Input file VRTX with the vertices: C (1) None to several strings terminated by / (a slash) C (2) For each vertex data (2.1): C (2.1) 'NAME',X1,X2,X3,/ C 'NAME'... Name of the vertex. Not considered. May be blank. C X1,X2,X3... Coordinates of the vertex. C /... None to several values terminated by a slash. C (3) / or end of file. C C C Input file TRGL with the triangles: C (1) For each triangle data (1.1): C (1.1) I1,I2,I3,/ C I1,I2,I3... Indices of 3 vertices of the triangle, right-handed C with respect to the given surface normals. C The vertices in file VRTX are indexed by positive integers C according to their order. C /... List of vertices is terminated by a slash. C (2) / or end of file. C C C Output formatted file SRFSECT (format LIN): C (1) / (a slash). C (2) For each triangle from the input file TRGL whose side(s) have a C crossection with the parallelogram given by points A,B and C C (2.1), (2.2), and (2.3): C (2.1) 'wf',/ (text 'wf' and a slash) C (2.2) Coordinates of the endpoints of the crossection segment: C (2.2.1) X1,X2,X3,/ C X1,X2,X3... Coordinates of the first endpoint. C /... A slash. C (2.2.2) Y1,Y2,Y3,/ C Y1,Y2,Y3... Coordinates of the second endpoint. C /... A slash. C (2.3) / (a slash). C (3) / (a slash). C General description of format LIN C C======================================================================= C C Common block /RAMC/: INCLUDE 'ram.inc' C ram.inc C INTEGER IRAM(MRAM) EQUIVALENCE (IRAM,RAM) C C....................................................................... C C External functions and subroutines: EXTERNAL LENGTH,RSEP1,RSEP3T,RSEP3R,ERROR,FORM1 INTEGER LENGTH C C Filenames and parameters: CHARACTER*80 FSEP,FVRTX,FTRGL,FOUT INTEGER LU1,LU2,LU3,IUNDEF PARAMETER (LU1=1,LU2=2,LU3=3,IUNDEF=-999999) C Other variables: CHARACTER*255 TEXT REAL U1,U2,U3,V1,V2,V3,A1,A2,A3,B1,B2,B3,C1,C2,C3, * BA1,BA2,BA3,CA1,CA2,CA3,AA,BB,CC,DD LOGICAL LCROS INTEGER NVRTX,NPLGN,OPLGN,NQ,I1,I C LCROS.. Says whether a crossection of the current triangle with C the parallelogram was already found. C C....................................................................... C C Reading main input data: WRITE(*,'(A)') '+SRFSECT: Enter input filename: ' FSEP=' ' READ (*,*) FSEP IF(FSEP.EQ.' ') THEN C SRFSECT-01 CALL ERROR('SRFSECT-01: No input file specified') C Input file in the form of the SEP (Stanford Exploration Project) C parameter or history file must be specified. C There is no default filename. END IF CALL RSEP1(LU1,FSEP) WRITE(*,'(A)') '+SRFSECT: Working... ' C C Reading input and output filenames: CALL RSEP3T('VRTX' ,FVRTX,'vrtx.out') CALL RSEP3T('TRGL' ,FTRGL,'trgl.out') CALL RSEP3T('SRFSECT',FOUT,'srfsect.out' ) C C Reading input parameters for the intersection parallelogram: CALL RSEP3R('A1',A1,0.) CALL RSEP3R('A2',A2,0.) CALL RSEP3R('A3',A3,0.) CALL RSEP3R('B1',B1,0.) CALL RSEP3R('B2',B2,0.) CALL RSEP3R('B3',B3,0.) CALL RSEP3R('C1',C1,0.) CALL RSEP3R('C2',C2,0.) CALL RSEP3R('C3',C3,0.) C The three points must be mutually different: IF (((A1.EQ.B1).AND.(A2.EQ.B2).AND.(A3.EQ.B3)).OR. * ((A1.EQ.C1).AND.(A2.EQ.C2).AND.(A3.EQ.C3)).OR. * ((B1.EQ.C1).AND.(B2.EQ.C2).AND.(B3.EQ.C3))) THEN C SRFSECT-02 CALL ERROR('SRFSECT-02: Wrong input parallelogram') C The points A,B and C must be mutually different. ENDIF C C Reading vertices: NQ=3 OPEN(LU1,FILE=FVRTX,STATUS='OLD') READ(LU1,*) (TEXT,I=1,20) NVRTX=0 20 CONTINUE IF(NVRTX+NQ.GT.MRAM) THEN C SRFSECT-03 CALL ERROR('SRFSECT-03: Too small array RAM') END IF TEXT='$' DO 21 I=NVRTX+1,NVRTX+NQ RAM(I)=0. 21 CONTINUE READ(LU1,*,END=29) TEXT,(RAM(I),I=NVRTX+1,NVRTX+NQ) IF(TEXT.EQ.'$') THEN GO TO 29 END IF C Number of storage locations in RAM used for the vertices NVRTX=NVRTX+NQ GO TO 20 29 CONTINUE CLOSE(LU1) C C Reading the triangles: DO 81 I=NVRTX+1,MRAM IRAM(I)=0 81 CONTINUE OPEN(LU1,FILE=FTRGL) NPLGN=0 OPLGN=NVRTX 82 CONTINUE IF(OPLGN+NPLGN+3.GT.MRAM) THEN C SRFSECT-04 CALL ERROR('SRFSECT-04: Too small array RAM') END IF IRAM(OPLGN+NPLGN+1)=IUNDEF IRAM(OPLGN+NPLGN+2)=IUNDEF IRAM(OPLGN+NPLGN+3)=IUNDEF READ(LU1,*,END=89) (IRAM(I),I=OPLGN+NPLGN+1,OPLGN+NPLGN+3) IF(IRAM(OPLGN+NPLGN+1).EQ.IUNDEF) THEN C End of triangles: GO TO 89 END IF IF ((IRAM(OPLGN+NPLGN+2).EQ.IUNDEF).OR. * (IRAM(OPLGN+NPLGN+3).EQ.IUNDEF)) THEN C SRFSECT-05 WRITE(TEXT,'(A,I6)') * 'SRFSECT-05: Missing vertices in triangle',NPLGN/3+1 CALL ERROR(TEXT(1:LENGTH(TEXT))) END IF IF ((IRAM(OPLGN+NPLGN+1).EQ.IRAM(OPLGN+NPLGN+2)).OR. * (IRAM(OPLGN+NPLGN+1).EQ.IRAM(OPLGN+NPLGN+3)).OR. * (IRAM(OPLGN+NPLGN+2).EQ.IRAM(OPLGN+NPLGN+3))) THEN C SRFSECT-06 WRITE(TEXT,'(A,I6)') * 'SRFSECT-06: Equal vertices in triangle',NPLGN/3+1 CALL ERROR(TEXT(1:LENGTH(TEXT))) END IF NPLGN=NPLGN+3 GO TO 82 89 CONTINUE CLOSE(LU1) C C Calculating the crossections: OPEN(LU3,FILE=FOUT) WRITE(LU3,'(A)') '/' C Parametric equation of the parallelogram: C X = A + s*[B-A] + t*[C-A] BA1=B1-A1 BA2=B2-A2 BA3=B3-A3 CA1=C1-A1 CA2=C2-A2 CA3=C3-A3 C General equation of the plane with the parallelogram: C AA*X1 + BB*X2 + CC*X3 + DD = 0 AA=BA2*CA3-BA3*CA2 BB=BA3*CA1-BA1*CA3 CC=BA1*CA2-BA2*CA1 DD=-1.*(AA*A1+BB*A2+CC*A3) C Loop over the triangles: DO 200, I1=OPLGN,OPLGN+NPLGN-3,3 LCROS=.FALSE. C First side of the triangle: U1=RAM((IRAM(I1+1)-1)*NQ+1) U2=RAM((IRAM(I1+1)-1)*NQ+2) U3=RAM((IRAM(I1+1)-1)*NQ+3) V1=RAM((IRAM(I1+2)-1)*NQ+1) V2=RAM((IRAM(I1+2)-1)*NQ+2) V3=RAM((IRAM(I1+2)-1)*NQ+3) CALL SRFCRO(U1,U2,U3,V1,V2,V3, * A1,A2,A3,BA1,BA2,BA3,CA1,CA2,CA3,AA,BB,CC,DD,LCROS,LU3) C Second side of the triangle: U1=RAM((IRAM(I1+1)-1)*NQ+1) U2=RAM((IRAM(I1+1)-1)*NQ+2) U3=RAM((IRAM(I1+1)-1)*NQ+3) V1=RAM((IRAM(I1+3)-1)*NQ+1) V2=RAM((IRAM(I1+3)-1)*NQ+2) V3=RAM((IRAM(I1+3)-1)*NQ+3) CALL SRFCRO(U1,U2,U3,V1,V2,V3, * A1,A2,A3,BA1,BA2,BA3,CA1,CA2,CA3,AA,BB,CC,DD,LCROS,LU3) C Third side of the triangle: U1=RAM((IRAM(I1+2)-1)*NQ+1) U2=RAM((IRAM(I1+2)-1)*NQ+2) U3=RAM((IRAM(I1+2)-1)*NQ+3) V1=RAM((IRAM(I1+3)-1)*NQ+1) V2=RAM((IRAM(I1+3)-1)*NQ+2) V3=RAM((IRAM(I1+3)-1)*NQ+3) CALL SRFCRO(U1,U2,U3,V1,V2,V3, * A1,A2,A3,BA1,BA2,BA3,CA1,CA2,CA3,AA,BB,CC,DD,LCROS,LU3) IF (LCROS) THEN WRITE(LU3,'(A)') '/' ENDIF 200 CONTINUE WRITE(LU3,'(A)') '/' C WRITE(*,'(A)') '+SRFSECT: Done. ' STOP END C C======================================================================= C SUBROUTINE SRFCRO(U1,U2,U3,V1,V2,V3, * A1,A2,A3,BA1,BA2,BA3,CA1,CA2,CA3,AA,BB,CC,DD,LCROS,LU3) INTEGER LU3 REAL U1,U2,U3,V1,V2,V3, * A1,A2,A3,BA1,BA2,BA3,CA1,CA2,CA3,AA,BB,CC,DD LOGICAL LCROS C C Subroutine designed to determine whether the abscissa formed by points C U and V intersects the parallelogram defined by points A, B and C. C C Input: C LU3 ... Number of logical unit connected to the output file. C U1,U2,U3,V1,V2,V3... Coordinates of the points U and V. C A1,A2,A3... Coordinates of the point A. C BA1,BA2,BA3,CA1,CA2,CA3... Coordinates of the vectors B-A and C-A. C AA,BB,CC,DD... Coefficients of the general equation of the plane C formed by points A, B and C. C LCROS... Identifies whether any crossection was already found. C C----------------------------------------------------------------------- C C Common block /RAMC/: INCLUDE 'ram.inc' C ram.inc C INTEGER IRAM(MRAM) EQUIVALENCE (IRAM,RAM) C C....................................................................... C Local storage locations: INTEGER J1,J2 CHARACTER*(25) FORMAT REAL VU1,VU2,VU3,W,SS,TT,X1,X2,X3,X(3),A(3),BA(3),CA(3),AUX1,AUX2 C C....................................................................... C C Parametric equation of the side: C X = U + w*[V-U] VU1=V1-U1 VU2=V2-U2 VU3=V3-U3 C Value of w for the crossection point: AUX1=(VU1*AA+VU2*BB+VU3*CC) IF (AUX1.EQ.0.) THEN C Vector VU is parallel to the parallelogram: AUX2=AA*X1+BB*X2+CC*X3+DD IF (AUX2.EQ.0.) THEN C Abscissa VU lies in the plane of the parallelogram. C Such abscissa is excluded in the current version of the code. W=-999. ELSE C Abscissa VU lies away of the plane of the parallelogram. W=-999. ENDIF ELSE W=-1.*(U1*AA+U2*BB+U3*CC+DD)/AUX1 ENDIF IF ((W.LE.1.).AND.(W.GE.0.)) THEN C Crossection of the triangle side with the parallelogram plane: X1=U1+W*VU1 X2=U2+W*VU2 X3=U3+W*VU3 C Computing values of parameters s and t for the parallelogram: C Selecting the two of the three parametric equations of the C parallelogram with at least one nonzero coefficient for s or t: J1=0 J2=0 IF ((BA1.NE.0.).OR.(CA1.NE.0.)) THEN J1=1 ENDIF IF ((BA2.NE.0.).OR.(CA2.NE.0.)) THEN IF (J1.EQ.0) THEN J1=2 IF ((BA3.NE.0.).OR.(CA3.NE.0.)) THEN J2=3 ELSE C SRFSECT-07 CALL ERROR('SRFSECT-07: Error in parametric equation') ENDIF ELSE J2=2 ENDIF ELSE IF (((BA3.NE.0.).OR.(CA3.NE.0.)).AND.(J1.NE.0)) THEN J2=3 ELSE C SRFSECT-08 CALL ERROR('SRFSECT-08: Error in parametric equation') ENDIF ENDIF A(1)=A1 A(2)=A2 A(3)=A3 BA(1)=BA1 BA(2)=BA2 BA(3)=BA3 CA(1)=CA1 CA(2)=CA2 CA(3)=CA3 X(1)=X1 X(2)=X2 X(3)=X3 C Calculating parameters s and t: SS=(X(J1)*CA(J2)-X(J2)*CA(J1)-A(J1)*CA(J2)+A(J2)*CA(J1))/ * (BA(J1)*CA(J2)-BA(J2)*CA(J1)) TT=(X(J1)*BA(J2)-X(J2)*BA(J1)-A(J1)*BA(J2)+A(J2)*BA(J1))/ * (CA(J1)*BA(J2)-CA(J2)*BA(J1)) C Values of s and t for the parallelogram must be within [0,1]: IF ((SS.GE.0.).AND.(SS.LT.1.)) THEN IF ((TT.GE.0.).AND.(TT.LT.1.)) THEN IF (.NOT.LCROS) THEN C First crossection of this triangle: WRITE(LU3,'(2A)') '''wf''',' /' LCROS=.TRUE. ENDIF FORMAT='(3(' CALL FORM1(AMIN1(X1,X2,X3),AMAX1(X1,X2,X3),FORMAT(4:11)) FORMAT(11:12)='))' WRITE(LU3,FORMAT) X1,' ',X2,' ',X3,' /' ENDIF ENDIF ENDIF RETURN END C C======================================================================= C INCLUDE 'error.for' C error.for INCLUDE 'sep.for' C sep.for INCLUDE 'length.for' C length.for INCLUDE 'forms.for' C forms.for C C======================================================================= C