PROGRAM SOD2COD
c last revision 13mar90 JSK, 18jul89 UFK,13dec94,17apr95,24feb96 UFK
C ***************************************************************
C                                                               *
C       SOD TO COD FORMAT OF ODF                                *
C       JOHN KALLEND SEPT 1987 & JUN 1988                       *
C                    (Glenwood)  (Los Alamos)                   *
C                    (  IL    )  (   NM     )                   *
C                                                               *
C       This version maintains the data in whatever format      *
c       (Kocks, Roe/Matthies or Bunge) they arrive in.          *
c                                                               *
C****************************************************************

      INTEGER*2 W(0:71,0:18,0:36),nsec,maxaz,isamp
      CHARACTER TITL*78,stuff*20,nomen*1,fname*40,secid*5

      PRINT *,'Convert SOD to VTK presentation of OD or vice versa'
      print*,'can read in the output .vtk file to Paraview'
      PRINT *,' Program by John Kallend, modified by ADR'
2000    FORMAT(/,A)

        CALL READSOD(NSEC,W,fname,TITL,stuff,nomen,MAXAZ,secid)
        PRINT 1000,NSEC
1000    FORMAT(1X,I3,' SECTIONS READ IN')
        IF(NSEC.GT.45)THEN
          ISAMP=3
        ELSE
          IF(NSEC.GT.25) THEN
            ISAMP=1
          ELSE
            ISAMP=0
          ENDIF
        ENDIF
        CALL vtkOUT(W,fname,TITL,ISAMP,stuff,nomen,MAXAZ,secid)
        END

C       ***************************************
        SUBROUTINE READSOD(I,W,fname,TITLE,stuff,nomen,MAXAZ,secid)
C               (may also read COD...  13dec94 UFK)
c       note; in this program, no conversion of incoming data to Roe
c       angles is performed.

        IMPLICIT INTEGER*2 (I-N)
        LOGICAL*1 TRIGNL,ORTHO
        INTEGER*2 W(0:71,0:18,0:36)
        integer*2 psi
        CHARACTER NAME*40,FNAME*40,TITLE*78,SECID*5,stuff*20
        character nomen*1,secin*5
        DATA TRIGNL/.FALSE./,ORTHO/.FALSE./
        logical sod2cod

C CODE::
23      PRINT '(A,$)','Input file (with .ext, default .SOD): '
        READ 1011,fname
1011    FORMAT(A)
c$$$      j=12
c$$$      do 100 i=1,12
c$$$      if(fname(i:i).eq.'.') j=i
c$$$ 100  continue
c$$$      if(j.eq.12) then
c$$$        name=fname(1:8)
c$$$        fname=name//'.SOD'
c$$$      else
c$$$        name=fname(1:j-1)
c$$$      endif

        sod2cod = .false.
        j = 0
        j = index(fname,'.sod')
        if(j.gt.0) then
           sod2cod = .true.
        else
           j = index(fname,'.SOD')
           if(j.gt.0) then
              sod2cod = .true.
           endif
        endif

        if(.not.sod2cod) then
           j = index(fname,'.cod')
           if(j.eq.0) then
              j = index(fname,'.COD')
              if(j.le.0) then
!                 stop 'neither SOD nor COD?!'
                 print*,'Careful: neither SOD nor COD'
              endif
           endif
        endif

        j = 0
        j = index(fname,'.')
        if(j.eq.0) then
           j = index(fname,' ')
        endif
        name = fname(1:j-1)
c$$$        if(sod2cod) then
c$$$           fname = name(1:j-1)//'.SOD'
c$$$        else
c$$$           fname = name(1:j-1)//'.COD'
c$$$        endif
        OPEN(FILE=FNAME,UNIT=3,STATUS='OLD',ERR=55)
c  now that we've opened the file, make the other name
c$$$        if(sod2cod) then
c$$$           fname = name(1:j-1)//'.COD'
c$$$        else
c$$$           fname = name(1:j-1)//'.SOD'
c$$$        endif
        fname = name(1:j-1)//'.vtk'

        I=0
12      IF(I.EQ.72)GOTO 35
        READ (3,'(A)',end=35)TITLE
        IF(I.EQ.0)PRINT '(/1X,A)',TITLE
        READ (3,1000)secin,PINC,PMAX,AINC,AZMAX,STUFF
1000    FORMAT(a5,4F5.1,A20)
1050    FORMAT(1X,19I4)
        if(secin(5:5).lt.'A') goto 35 !96
        secid=secin
        nomen=secid(5:)
        IF(I.EQ.0)THEN
           MAXAZ=AZMAX/5.+.01
           IF(MAXAZ.EQ.24)TRIGNL=.TRUE.
           IF(MAXAZ.EQ.36)ORTHO=.TRUE.
        ENDIF
        psi=i
        do 140 J=0,18
           IF(TRIGNL)THEN
              READ(3,'(1X,18I4)')(W(PSI,J,K),K=0,MAXAZ)
           ELSEIF(ORTHO) THEN
              READ(3,'(1X,18I4,/,1X,19I4)')(W(PSI,J,K),K=0,MAXAZ)
           ELSE
              read(3,1050) (W(psi,J,K), K=0,18)
           ENDIF
  140   continue
        READ(3,1011,END=35)
        I=I+1
        GOTO 12
35      CLOSE(UNIT=3)
        IF(I.LT.16)STOP'INCOMPLETE DATA FILE'
        if(secid(2:2).eq.'S') secid(2:2)='C'
        if(secin(2:2).eq.'C') secid(2:2)='S' !95
c      fname=name//'.'//secid(2:4)
c  taken care of above
        RETURN
55      PRINT*,' CANNOT OPEN FILE ',FNAME
        PRINT*
        GOTO 23
        END

C       ****************************************
        SUBROUTINE vtkOUT(W,fname,TITL,ISAMP,stuff,nomen,MAXAZ,secid)
C       OUTPUT COD (or SOD)
        IMPLICIT INTEGER*2(I-N)
        CHARACTER FNAME*40,TITL*78,SECID*5,YESNO
        character stuff*20,seclab*5,nomen*1,SECVAL*5
        INTEGER*2 W(0:71,0:18,0:36)
        INTEGER*4 PROJ(0:18,0:71)
        INTEGER*2 IQUAD(4)
        INTEGER*2 PHI,THET,WO(0:71),psi
        DATA IQUAD/4*18/

!  CODE::

        OPEN(UNIT=20,FILE=FNAME,STATUS='UNKNOWN')

        write(20,"('# vtk DataFile Version 2.0')")
!        write(20,"(' data set from ',a,a)") fname(1:dotpos+4),stamp
        write(20,"(' data set from ',a)") fname
        write(20,"('ASCII')") 
        write(20,"('DATASET STRUCTURED_POINTS')") 
c     write(20,"('DIMENSIONS ',3(2x,i8))") mx,my,num_z_values
        IF(ISAMP.EQ.0)THEN
           WRITE (20,"('DIMENSIONS ',3(2x,i8))") 19,19,19
        ENDIF
        IF(ISAMP.EQ.1)THEN
           WRITE (20,"('DIMENSIONS ',3(2x,i8))") 37,19,19
        ENDIF
        IF(ISAMP.EQ.3)THEN
           WRITE (20,"('DIMENSIONS ',3(2x,i8))") 72,19,19
        ENDIF
        write(20,"('ORIGIN 0.0 0.0 0.0')") 
        write(20,"('SPACING 1 1 1')") 
!        write(20,"('POINT_DATA ',i12)")xx*yy*zz
        IF(ISAMP.EQ.0)THEN
           WRITE (20,"('POINT_DATA ',i12)") 19*19*19
        ENDIF
        IF(ISAMP.EQ.1)THEN
           WRITE (20,"('POINT_DATA ',i12)") 37*19*19
        ENDIF
        IF(ISAMP.EQ.3)THEN
           WRITE (20,"('POINT_DATA ',i12)") 72*19*19
        ENDIF
        write(20,*)
        write(20,"('SCALARS Intensities float  1') ")
        write(20,"('LOOKUP_TABLE default')")

        print*,' Making file ',fname
        IF(ISAMP.EQ.0)THEN
           WRITE (20,"(10(1x,g10.3))")
     $          (((float(w(i,j,k))/100.,i=0,18),j=0,18),k=0,18)
        ENDIF
        IF(ISAMP.EQ.1)THEN
           WRITE (20,"(10(1x,g10.3))")
     $          (((float(w(i,j,k))/100.,i=0,36),j=0,18),k=0,18)
        ENDIF
        IF(ISAMP.EQ.3)THEN
           WRITE (20,"(10(1x,g10.3))")
     $          (((float(w(i,j,k))/100.,i=0,72),j=0,18),k=0,18)
        ENDIF
        write(20,*)

        RETURN
        END