C 
C Program CRTCART converting the rays calculated by program CRT from
C curvilinear to Cartesian coordinates.
C
C Version: 5.20
C Date: 1998, February 28
C
C Coded by: Ludek Klimes
C     Department of Geophysics, Charles University Prague,
C     Ke Karlovu 3, 121 16 Praha 2, Czech Republic,
C     E-mail: klimes@seis.karlov.mff.cuni.cz
C
C.......................................................................
C
C                                                    
C Description of data files:
C
C Main input data read from external interactive device (*):
C     The data consist of character strings, read by list directed (free
C     format) input.  The strings have thus to be enclosed in
C     apostrophes.  The interactive * external unit may be redirected to
C     the file containing the data.
C (1) 'MODEL','RAYALL','INIALL','RAY2P','INI2P',/
C     'MODEL'... Formatted file with the input data for the model.
C             Only integer KOORS specifying the type of the coordinate
C             system is read from data file MODEL.
C     'RAYALL'... Input file CRT-R with the quantities stored along rays
C             (see C.R.T.5.5.1), or input file CRT-S with the quantities
C             stored at the specified surfaces (see C.R.T.5.5.2).
C     'INIALL'... Input file CRT-I with the quantities at the initial
C             points of rays, corresponding to file RAYALL (see
C             C.R.T.6.1).
C     'RAY2P'... Output file CRT-R or CRT-S, with the quantities
C             corresponding to the 2-point rays only.
C     'INI2P'... Output file CRT-I with the quantities at the initial
C             points of 2-point rays, corresponding to file RAY2P (see
C             C.R.T.6.1).
C Default: 'MODEL'='model.dat', 'RAYALL'='r01.out', 'INIALL'='r01i.out',
C     'RAY2P'='ray2p.out', 'INI2P'='ini2p.out'.
C
C Formatted file MODEL:
C     See the description within source code file 'model.for'.
C     Description of data MODEL 
C
C Unformatted files RAYALL and RAY2P:
C     See the description within source code file 'writ.for'.
C     Description of files CRT-R
C     Description of files CRT-S
C
C Unformatted files INIALL and INI2P:
C     See the description within source code file 'writ.for'.
C     Description of files CRT-I
C
C-----------------------------------------------------------------------
C
C Common block /POINTC/ to store the results of complete ray tracing:
      INCLUDE 'pointc.inc'
C     pointc.inc
C None of the storage locations of the common block are altered.
C
C Subroutines and external functions required:
      EXTERNAL AP00
C     AP00... File 'ap.for'.
C
C-----------------------------------------------------------------------
C
C     Auxiliary storage locations:
      INTEGER LU1,LU2,LU3,LU4,I
      PARAMETER (LU1=1,LU2=2,LU3=3,LU4=4)
      CHARACTER*80 FILE0,FILE1,FILE2,FILE3,FILE4
      CHARACTER*1 TEXTM
C
C.......................................................................
C
C     Opening input and output files:
      FILE0='model.dat'
      FILE1='r01.out'
      FILE2='r01i.out'
      FILE3='ray2p.out'
      FILE4='ini2p.out'
      WRITE(*,'(2A)')
     *  ' Enter 4 filenames',
     *  ' (model.dat, r01.out, r01i.out, ray2p.out, ini2p.out): '
      READ(*,*) FILE0,FILE1,FILE2,FILE3,FILE4
      OPEN(LU1,FILE=FILE0,STATUS='OLD')
      READ(LUN,*) TEXTM
      I=0
      READ(LUN,*) I
      CALL METR1(I)
      CLOSE(LU1)
      OPEN(LU1,FILE=FILE1,FORM='UNFORMATTED',STATUS='OLD')
      OPEN(LU2,FILE=FILE2,FORM='UNFORMATTED',STATUS='OLD')
      OPEN(LU3,FILE=FILE3,FORM='UNFORMATTED')
      OPEN(LU4,FILE=FILE4,FORM='UNFORMATTED')
C
C     Loop for the points of rays
   10 CONTINUE
C       Reading the results of the complete ray tracing
        CALL AP00(0,LU1,LU2)
        IF(IWAVE.LT.1)THEN
C         End of rays
          GO TO 80
        END IF
C       Conversion to Cartesian coordinates
        CALL TOCART(Y,YL)
        WRITE(LU3) IWAVE,IRAY,NY,ICB1,ISRF,X,YL,(Y(I),I=1,NY)
        IF(IPT.LE.1)THEN
C         New ray - recording the initial point
          CALL TOCART(YI,YLI)
          WRITE(LU4) -IWAVE,IRAY,ICB1I,IEND,ISHEET,IREC,YLI,YI
        END IF
      GO TO 10
C
   80 CONTINUE
      CLOSE(LU1)
      CLOSE(LU2)
      CLOSE(LU3)
      CLOSE(LU4)
      STOP
      END
C
C=======================================================================
C
C     
C
      SUBROUTINE TOCART(Y,YL)
      REAL Y(11),YL(6)
C
C This subroutine transforms points, slowness vectors, R.C.C.S. basis
C vector and velocity gradient from the model coordinates to Cartesian
C coordinates.
C
C Input and output:
C     Y...    Array containing basic quantities computed along the ray
C             (C.R.T.5.2.1).
C             Description of Y
C     YL...   Array containing local quantities at the point of the ray
C             (C.R.T.5.5.4).
C             Description of YL
C
C Subroutines and external functions required:
      EXTERNAL CARTES
C
C Date: 1998, February 28
C Coded by Ludek Klimes
C
C-----------------------------------------------------------------------
C
C     Auxiliary storage locations:
      REAL CART(3),CART1,CART2,CART3
      REAL PDER(9),PDER1,PDER2,PDER3,PDER4,PDER5,PDER6,PDER7,PDER8,PDER9
      EQUIVALENCE (CART(1),CART1),(CART(2),CART2),(CART(3),CART3)
      EQUIVALENCE (PDER(1),PDER1),(PDER(2),PDER2),(PDER(3),PDER3)
      EQUIVALENCE (PDER(4),PDER4),(PDER(5),PDER5),(PDER(6),PDER6)
      EQUIVALENCE (PDER(7),PDER7),(PDER(8),PDER8),(PDER(9),PDER9)
C
C.......................................................................
C
C     Conversion of model coordinates to Cartesian coordinates CART:
      CALL CARTES(Y(3),.TRUE. ,CART,PDER)
C     Transposed transformation matrix PDER of covariant vectors:
      CALL CARTES(Y(3),.FALSE.,CART,PDER)
C
C     Transformations:
      Y( 3)=CART1
      Y( 4)=CART2
      Y( 5)=CART3
      CART1=PDER1*Y( 6)+PDER2*Y( 7)+PDER3*Y( 8)
      CART2=PDER4*Y( 6)+PDER5*Y( 7)+PDER6*Y( 8)
      CART3=PDER7*Y( 6)+PDER8*Y( 7)+PDER9*Y( 8)
      Y( 6)=CART1
      Y( 7)=CART2
      Y( 8)=CART3
      CART1=PDER1*Y( 9)+PDER2*Y(10)+PDER3*Y(11)
      CART2=PDER4*Y( 9)+PDER5*Y(10)+PDER6*Y(11)
      CART3=PDER7*Y( 9)+PDER8*Y(10)+PDER9*Y(11)
      Y( 9)=CART1
      Y(10)=CART2
      Y(11)=CART3
      CART1=PDER1*YL(4)+PDER2*YL(5)+PDER3*YL(6)
      CART2=PDER4*YL(4)+PDER5*YL(5)+PDER6*YL(6)
      CART3=PDER7*YL(4)+PDER8*YL(5)+PDER9*YL(6)
      YL(4)=CART1
      YL(5)=CART2
      YL(6)=CART3
      END
C
C=======================================================================
C
      INCLUDE 'error.for'
C     error.for
      INCLUDE 'metric.for'
C     metric.for
      INCLUDE 'ap.for'
C     ap.for
C
C=======================================================================
C