C     PROGRAM ATTACH-2-EPIC-to-REG-2013
C
C                    Program ATTACH-2-EPIC-to-REG-2013		
C***********************************************************************
C	              	Author: GUENTER LEYDECKER
C			       Isernhagen, Germany
C                                April 04, 2013
C
C
C	ATTACHES EACH EPICENTER TO HIS SEISMOGEOGRAPHICAL REGION
C
C                          input           output
C		  ===========================================
C                   Epi.coordinates ---> Seismogeogr. Region
C                      -----------------------------
C                           File    --->    File
C                 or
C                          Screen   --->   Screen
C
C	In this program version the following limits are fixed:
C!!!!!	maximum number of regions STEZAHL = 200
C!!!!!	maximum number of polygonvertices STPOLY(I) for i'th region  = 30
C
C!!!!!	the computation of the cartesian coordinates in the 
C!!!!!	subroutines GEOKOO and GEOKO2 is related to
C!!!!!		50 degree N  and  10 degree E 
C!!!!!	as origin of the new coordinate system
C
C
      CHARACTER*2 STNAME(200),STNAM,CODE2,IDA23,BLANC2
      CHARACTER*35 KNAME(200),KNAM
      CHARACTER*80 FILE,OFILE,CODE1
      INTEGER STEZAHL,STPOLY(200)
      DIMENSION STLAT(200,30),STLONG(200,30),PX(30),PY(30)
      DIMENSION IDA(80)
      DATA LUI/2/,LUO/4/,BLANC2/'  '/
C
C
1108  WRITE (*,400) 
400   FORMAT(//,8X,' Attaches each EPICENTER to his SEISMOGEOGRAPHI',
     * 'CAL REGION ',/,9x,56(1H=),//,5X,
     *' How do you wish the input of the epicenter coordinates?'
     *//,5X,' 1. In- and Output on SCREEN  =  1  ',/,5X,
     *      ' 2. In- and Output on FILES   =  2',/,9X,
     *          '(adding only the abbreviation (2 letters) of the',/,
     *9X,       'Seismogeographical Region to output file)',///,
     *5X,'ATTENTION: in this program version, the COORDINATES of ',/,
     *5X,'Earthquake Data in IN- and OUTput FILE are in DEGREE AND ',
     *'MINUTES!!',/,5X,'If you  choose IN- and OUTput on Screen,',/, 
     *5X,'the input coordinates can be in Degree and Minutes or in ',
     *'Degree only',//,' Please give your decision: input 1 or 2') 
      READ (*,*) INOUT
      IF (INOUT .LT. 1 .OR. INOUT .GT. 2) GO TO 1108
      IF (INOUT .EQ. 2) THEN
C
      WRITE (*,'(1X,A,$)') 'INPUT FILE NAME of EARTHQUAKE DATA   '
      READ (*,'(Q,A)') N,FILE(:N)
      OPEN (UNIT=LUI,FILE=FILE,STATUS='OLD')
C
      WRITE (*,'(1X,A,$)') 'OUTPUT FILE NAME of EARTHQUAKE DATA   '
      READ (*,'(Q,A)') NP,OFILE(:NP)
      OPEN (UNIT=LUO,FILE=OFILE,STATUS='NEW',CARRIAGECONTROL='LIST')
C
      end if
C
C     INPUT: Data file of the SEISMOGEOGRAPHICAL REGIONS
C       and transformation from geographical to cartesian coordinates       
C
      CALL STEKIN (STLAT,STLONG,STNAME,STPOLY,STEZAHL,KNAME)
C
C     INPUT: list of earthquake epicenters
C
C	If input by file is selected, input data are in DEGREE and MINUTES
      IGRAD = 1
C
      IF (INOUT .EQ. 2) GO TO 1110
C
1107  WRITE (*,*) ' '
      WRITE(*,*)' INPUT of epic. coordinates in DEGREE and MINUTES',
     *  '  = 1'
      WRITE(*,*)' INPUT of epic. coordinates in DEGREE            ',
     *  '  = 2'
      READ *, IGRAD
      IF (IGRAD .LT. 1 .OR. IGRAD .GT. 2) GO TO 1107
C
1109  CONTINUE            !START new computation with input on screen
C
      IF (IGRAD .EQ. 1) THEN
      WRITE (*,*)(' LATITUDE in DEGREE and MINUTES (free format)')
      READ *, CLAT, YMIN
      WRITE (*,*)(' LONGITUDE in DEGREE and MINUTES (free format)')
      READ *, CLON, XMIN
      WRITE (*,333)
333   FORMAT (70(1H*),/)
      WRITE (*,334) CLAT,YMIN,CLON,XMIN
334     FORMAT (9X,' The Epicenter with the Coordinates',//,10X,F4.0,
     *  ' degr. ',F6.3,' min. N and ',F5.0,' degr. ',F6.3,' min. E',
     *  //,5X,' is situated in the SEISMOGEOGRAPHICAL REGION ')
      CALL GEOKO2 (CLAT,YMIN,CLON,XMIN,YLAT,XLONG)
C
      ELSE IF (IGRAD .EQ. 2) THEN
      WRITE (*,*)(' LATITUDE in DEGREE (free format)')
      READ *, CLAT
      WRITE (*,*)(' LONGITUDE in DEGREE (free format)')
      READ *, CLON
      WRITE (*,333)
      WRITE (*,335) CLAT,CLON
335     FORMAT (9X,' The Epicenter with the Coordinates',//,10X,F8.4,
     *  ' degree N and ',F9.4,' degree E',//,5X,
     *  ' is situated in the SEISMOGEOGRAPHICAL REGION ')
      CALL GEOKOO (CLAT,CLON,YLAT,XLONG)
C
      ENDIF
      GO TO 7777
C
C	INPUT: File with the earthquake data (catalogue)
C              
C       The following input file has the format of the German 
C       Earthquake Catalogue. 
C       CLAT = latitude N in degree
C       YMIN = latitude N in minutes 
C       CLON = longitude E or W in degree
C       XMIN = longitude E or W in minutes
C       LONS = "sign" of longitude: 0 or blanc = E
C                                   1          = W
C
1110  OPEN(UNIT=LUI,FILE=FILE,STATUS='OLD')
1111  READ(LUI,500,END=90) (IDA(K),K=1,80) 
500   FORMAT (7A2,A1,2(2A2,A1),8A1,3A2,A1,5A2,A1,A2,5(A1,A2),A2,    !44
     *7A1,A2,A1,27A1)                                               !36
C
       WRITE (CODE1,700,ERR=888)IDA(1),IDA(2),IDA(9),IDA(10),IDA(11),
     *IDA(12),IDA(13),IDA(14),IDA(15),IDA(27),IDA(31)
700    FORMAT(2A2,2(A2,A2,A1),A1,2A2)
      GO TO 892
888    WRITE(LUO,889)
       WRITE(*,889)
889    FORMAT(5X,12HENCODE ERROR)
C
892    READ (CODE1,701,ERR=890)JAHR,CLAT,YMIN,YZ,CLON,XMIN,XZ,LONS,
     *XML,XI0
701    FORMAT(I4,2(F2.0,F2.0,F1.0),I1,2F2.1)
C
         YMIN=YMIN+(YZ/10.)
	 XMIN=XMIN+(XZ/10.)
C	
      IF (LONS .EQ. 1) THEN		!Western longitude
          CLON = -CLON
          XMIN = -XMIN
      END IF
CC
       GO TO 893
890    WRITE(LUO,891)
       WRITE(*,891)
891    FORMAT(5X,12HDECODE ERROR)
C
C	TRANSFORMATION of the coordinates from DEGREE and MINUTES 
C       (geographic coordin.) in km (YLAT,XLONG) (cartesian coordin.)
C	
893   CALL GEOKO2 (CLAT,YMIN,CLON,XMIN,YLAT,XLONG)
C
C
C       Compute in which SEISMOGEOGRAPHICAL REGION  the 
C       EPICENTER is situated 
C
7777  DO 55 I=1,STEZAHL
      IPOLY=STPOLY(I)
C	store the coordinates of the i-th polygon in new variables
          DO 30 M=1,IPOLY
          PX(M)=STLONG(I,M)
          PY(M)=STLAT(I,M)
30    continue
C
      INSI=INSIDE (XLONG,YLAT,PX,PY,IPOLY)
C       return value:   0  if point outside
C		     +- 1  if point inside
C			2  if point is on an edge or vertex
C
      IF (INOUT .EQ. 1) GO TO 3333               !output screen
C
C	OUTPUT on FILE
C
      IF (INSI .EQ. -1 .OR. INSI .EQ. 1 .OR. INSI .EQ. 2) THEN
          WRITE (CODE2,710,ERR=790) STNAME(I)
710       FORMAT (A2)
        	GO TO 785
790       WRITE(LUO,889)
          WRITE(*,889)
785       READ (CODE2,711,ERR=791) NUM
711       FORMAT (A2)
          GO TO 786
791       WRITE(LUO,889)
          WRITE(*,889)
786       IDA(23)=NUM                 !seismogeographical region
          FLAG=1
          GO TO 60
      END IF
      GO TO 55
C
C	OUTPUT on SCREEN
3333  IF (INSI .EQ. -1 .OR. INSI .EQ. 1 .OR. INSI .EQ. 2) THEN
          STNAM=STNAME(I)     !abbreviation of seismogeographical region
          KNAM =KNAME(I)      !full name of seismogeographical region
      WRITE (*,336) STNAM,KNAM
336   FORMAT(/,15X,A,5X,A)
      WRITE (*,333)
          FLAG=1
          GO TO 60
      END IF
      GO TO 55
C
55    CONTINUE
C	
C       End of loop reached, NO seismogeographical region attachable
      FLAG=0
C
60    CONTINUE
C
C	OUTPUT: 
C       Earthquake data file with the attached seismogeographical regions
C	
      IF (INOUT .EQ. 1 .AND. FLAG .EQ. 0) THEN
          WRITE (*,*) '            ****   OUT OF RANGE   ****'
          WRITE (*,333)
      ENDIF
C
      IF (INOUT .EQ. 1) THEN
          WRITE (*,*) ' CONTINUE?  YES = 1;   NO = 0'
          READ *, ICONT
          IF (ICONT .EQ. 0) GO TO 9999
          GO TO 1109
      ENDIF
101   CONTINUE
      WRITE (LUO,500)(IDA(K),K=1,80)
C
      WRITE (IDA23,711,ERR=7901) IDA(23)
          GO TO 7855
7901  WRITE(LUO,889)
      WRITE(*,889)
C
7855  CONTINUE		!neu, verhindert das Schreiben "OUT OF RANGE"
C7855	IF (FLAG .EQ. 0 .AND. IDA23 .EQ. BLANC2) THEN
C		WRITE (LUO,*) '****   OUT OF RANGE   ****'
C	END IF
C
      GO TO 1111
90    CLOSE(UNIT=LUI,STATUS='KEEP')
      CLOSE (UNIT=LUO,STATUS='KEEP')
 9999 STOP
      END
C **********************************************************************
      SUBROUTINE STEKIN (STLAT,STLONG,STNAME,STPOLY,STEZAHL,KNAME)
C
C	to read the datafile with the SEISMOGEOGRAPHICAL REGIONS
C
C     input: polygon coordinates in degree ----> output: cartesian
C                                                coordin. in km
C               -------------------------------------------------
C	Format: Latitude      LAT          ----> STLAT(I,K)    
C		Longitude     LONG         ----> STLONG(I,K)    
C               plotcommand   ISTIFT                   
C		abbrev. name  NAM          ----> STNAME(I)     
C	         full name    KNAM         ----> KNAME(I)
C   last line of a polygon: NAM = "**"
C		
C           number of polygonvertices            STPOLY(I)
C	    number of seismogeographical regions STEZAHL
C
      CHARACTER STNAME(1)*2,NAM*2
      CHARACTER KNAME(1)*35, KNAM*35
      CHARACTER*80 FILE2
      REAL LAT,LONG
      INTEGER STEZAHL,STPOLY(200)
      DIMENSION STLAT(200,30),STLONG(200,30)
      DATA LUI2/12/
C
C	INPUT: coordinates of the corners of the polygons 
C
      WRITE (*,*) ' '
      WRITE (*,'(1X,A,$)') 'INPUT FILE with the SEISMOGEOGRAPHICAL',
     *  ' REGIONS    '
      READ (*,'(Q,A)') N,FILE2(:N)
      WRITE (*,*) ' '
      OPEN(UNIT=LUI2,FILE=FILE2,STATUS='OLD')
C
C       NN = number of the seismogeographical regions
      NN=0                            
      DO 1 N=1,200
      NN=NN+1
C
      K=0
      DO 5 I=1,30
      K = K+1
      IF (I .EQ. 1) THEN
      READ (LUI2,*,END=99) LAT,LONG,ISTIFT,NAM,KNAM
      ELSE 
      READ (LUI2,*,END=99) LAT,LONG,ISTIFT,NAM
        ENDIF
C
      IF (NAM .EQ. '**') GO TO 6          !last=first polygon corner
          IF (I .EQ. 1) THEN
              STNAME(N)=NAM       !abbreviation of name (2 letters)
              KNAME(N) =KNAM      !full name
          TYPE*,' STNAME(N)=',stname(N),'  NAM= ',nam
          END IF
C
      CALL GEOKOO (LAT,LONG,STY,STX)
C
      STLAT(N,I)=STY
      STLONG(N,I)=STX
C
5     CONTINUE
C
6     K = K-1
C
      STPOLY(N)=K                         !number of polygonvertices
1     CONTINUE
99    STEZAHL=NN-1                 !number of seismogeographical regions
C
      CLOSE(UNIT=LUI2,STATUS='KEEP')
C
      RETURN
      END
************************************************************************
      SUBROUTINE GEOKOO (LAT,LONG,YLAT,XLONG)
C	for input of LAT and LONG in degree
C
C	OUTPUT: YLAT = y-Coordinate in km
C	        XLONG = x-Coordinate in km
C
C      TRANSFORMATION of geographical coordinates in cartesian  
C      coordinates after C.F.RICHTER: ELEMENTARY SEISMOLOGY. pp 701-704.
C      W.H.FREEMAN and CO. SAN FRANCISCO and LONDON. 1958.
C
      REAL LAT,LONG
      DATA PPP/1.859061/,QQ/1.853842/
C
C!!!!!    the values for PPP and QQ are for 50 degree latitude  !!!!!
C
      DATA LORED/10.0/,LARED/50.0/
C
C!!!!!	the computation of the cartesian coordinates is related to
C!!!!!		50 degree N  and  10 degree E 
C!!!!!	as origin of the new coordinate system
C
        ZLAT=(LAT)/(57.29578)
C	PP and QQ : km/minute(angel)
      PP=PPP*COS(ZLAT)                    
      DX=(LONG-LORED)*60.*PP                       !cartesian X-value
      DY=(LAT-LARED)*60.*QQ                        !cartesian Y-value
      IF(DX*DX .LT. 0.000001)DX=0.001
      IF(DY*DY .LT. 0.000001)DY=0.001
      YLAT=DY
      XLONG=DX
C
      RETURN
      END
************************************************************************
      SUBROUTINE GEOKO2 (LAT,YMIN,LONG,XMIN,YLAT,XLONG)
C	for input of LAT and LONG in degree, XMIN and YMIN in minutes
C
C	OUTPUT: YLAT = y-coordinate in km
C	       XLONG = x-coordinate in km
C
C      TRANSFORMATION of geographical coordinates in cartesian  
C      coordinates after C.F.RICHTER: ELEMENTARY SEISMOLOGY. pp 701-704.
C      W.H.FREEMAN and CO. SAN FRANCISCO and LONDON. 1958.
C
      REAL LAT,LONG
      DATA PPP/1.859061/,QQ/1.853842/
C       the values of PPP and QQ are for 50 degree latitude
      DATA LORED/10.0/,LARED/50.0/
C
C!!!!!	the computation of the cartesian coordinates is related to
C!!!!!		50 degree N  and  10 degree E 
C!!!!!	as origin of the new coordinate system
C
        ZLAT=(LAT*60. + YMIN)/(60.*57.29578)
      PP=PPP*COS(ZLAT)
      DX=(LONG*60.+XMIN - 60.*LORED)*PP       !cartesian X-value
      DY=(LAT*60.+YMIN - 60.*LARED)*QQ        !cartesian Y-value
      IF(DX*DX .LT. 0.000001)DX=0.001
      IF(DY*DY .LT. 0.000001)DY=0.001
C	
      YLAT=DY
      XLONG=DX
C
      RETURN
      END
C***************************************************************************
C***************************************************************************
      FUNCTION INSIDE (X0,Y0,PX,PY,N)
C	GODKIN,C.B. AND J.J.PULLI: APPLICATION OF THE "WINDING-NUMBER
C	ALGORITHM" TO THE SPATIAL SORTING OF CATALOGED EARTHQUAKE DATA.
C	BSSA 75, 5, PP. 1845-1848, OCTOBER 1984
C
C	CHECK IF POINT X0,Y0 LIES INSIDE POLYGONE
C
C	PARAMETERS:
C	X0,Y0    POINT TO TEST
C	PX       ARRAY OF X-COORDINATES OF POLYGON
C	PY       ARRAY OF Y-COORDINATES OF POLYGON
C       N        NUMBER OF ELEMENTS OF PX AND PY OR OF SIDES OF POLYGON
C
C	RETURN VALUE:    0  IF POINT OUTSIDE
C                     +/-1  IF POINT INSIDE 
C                        2  IF POINT IS ON AN EDGE OR VERTEX
C
      REAL PX(30),PY(30)
C
C	ACCUMULATE SIGNED CROSSING NUMBERS WITH INSIDE
      INSIDE=0
      DO 200 I=1,N-1
C	PROCEED AROUND POLYGON CHECKING EACH SEGMENT TO SEE IF 
C	NEGATIVE X-AXIS WAS CROSSED
C       TRANSLATE COORDINATES OF POLYGON TO PUT TEST POINT AT ORIGIN
      ISICR=KSICR(PX(I)-X0,PY(I)-Y0,PX(I+1)-X0,PY(I+1)-Y0)
C	type *, I,'PX=',PX(I),'PY=',PY(I),' X0=',X0,' Y0=',Y0
C	STOP COUNTING IF POINT IS ON EDGE
      IF(ISICR .EQ. 4) GO TO 600
      INSIDE=INSIDE+ISICR
200   CONTINUE
C
C       CHECK SEGMENT FROM LAST VERTEX TO FIRST VERTEX
      ISICR=KSICR(PX(N)-X0,PY(N)-Y0,PX(1)-X0,PY(1)-Y0)
      IF(ISICR .EQ. 4) GO TO 600
      INSIDE=(INSIDE+ISICR)/2
C	type *, ' N= ',N,' INSIDE=',INSIDE,' ISICR=',ISICR
C
      RETURN
C
600     INSIDE=2
C
C	type *, ' N= ',N,' INSIDE=',INSIDE,' ISICR=',ISICR
      RETURN
      END
C**********************************************************************
      FUNCTION KSICR(X1,Y1,X2,Y2)
C	COMPUTE SIGNED CROSSING NUMBER
C
C	A "SIGNED CROSSING NUMBER" TELLS WETHER A SEGMENT 
C	(IN THIS CASE THE SEGMENT FROM (X1,Y1) TO (X2,Y2))
C	CROSSES THE NEGATIVE X-AXIS OR GOES THROUGH THE ORIGIN
C
C	THE RETURN VALUES ARE:
C             +2 IF SEGMENT CROSSES FROM BELOW
C	      +1 IF SEGMENT EITHER ENDS ON -X-AXIS FROM BELOW OR STARTS
C                           UPWARDS FROM -X-AXIS ("HALF CROSSING")
C	       0 IF THERE IS NO CROSSING
C	      -1 IF SEGMENT EITHER ENDS ON -X-AXIS FROM ABOVE OR STARTS
C                           DOWNWARDS FROM -X-AXIS ("HALF CROSSING")
C             -2 IF SEGMENT CROSSES FROM ABOVE
C             +4 IF SEGMENT CROSSES THROUGH THE ORIGIN
C
C
C	IF BOTH POINTS ARE ON THE SAME SIDE OF X-AXIS, RETURN 0
      IF (Y1*Y2 .GT. 0) GO TO 600
C
C       CHECK IF SEGMENT CROSSES THROUGH THE ORIGIN
        IF (X1*Y2 .NE. X2*Y1 .OR. X1*X2 .GT. 0) GO TO 100
      KSICR=4
      RETURN
C       
C	CHECK FOR COMPLETE CROSSING
100     IF(Y1*Y2 .LT. 0) GO TO 300
C
C       HALF CROSSING?
C       ONE END OF SEGMENT TOUCHES X-AXIS! WHICH END?       
        IF (Y2 .EQ. 0) GO TO 200
C       HERE Y1=0   CHECK IF SEGMENT TOUCHES +X-AXIS
        IF (X1 .GT. 0) GO TO 600
C       UPWARD OR DOWNWARD?
        IF (Y2 .GT. 0) GO TO 700
      GO TO 800
C
C       HERE Y2=0   CHECK IF SEGMENT TOUCHES +X-AXIS
200     IF (Y1 .EQ. 0 .OR. X2 .GT. 0) GO TO 600
C       UPWARD OR DOWNWARD?
        IF (Y1 .GT. 0) GO TO 800
      GO TO 700
C
C       COMPLETE CROSSING OF -X-AXIS?
C       BREAK INTO CASES ACCORDING TO CROSSING DIRECTION
300   IF (Y1 .GT. 0) GO TO 400
C
C       CASE   Y1 < 0 < Y2
      IF (X1*Y2 .GE. Y1*X2) GO TO 600
C	UPWARD CROSSING
      KSICR=2
      RETURN
C
C       CASE   Y1 > 0 > Y2
400   IF (Y1*X2 .GE. X1*Y2) GO TO 600
C	DOWNWARD CROSSING
      KSICR=-2
      RETURN
C
C	NO CROSSING
600   KSICR=0
      RETURN
C
C	UPWARD HALF CROSSING
700   KSICR=1
      RETURN
C
C	DOWNWARD HALF CROSSING
800   KSICR=-1
      RETURN
      END