#define WITH_PAINT 1
#define OUTER_SURFACE_QUALITY 0
#define MIDDLE_SURFACE_QUALITY 0.99
#define INNER_SURFACE_QUALITY 0.99

      subroutine define_geometry()
      real elevation
      real arc_radius(2)
      real straight_length(2)
      real fiber_width(3)
      real arc_spread
      real arc_length
      real pi
      parameter (pi=3.141596)

      character*20 natmed
      integer nmat
      integer isvol
      integer ifield
      real fieldm
      real tmaxfd
      real stemax
      real deemax
      real epsil
      real stmin
      real ubuf(9)
      integer nwbuf
       
      integer npckov
      real ppckov(2)
      real absco(2)
      real effic(2)
      real rindex(2)

      character*4 chname
      character*4 chshape
      integer nmed
      real par(9)
      integer npar
      integer ivolu
      
      integer nr
      character*4 chmoth
      real x,y,z
      integer irot
      character*4 chonly

c define the geometry dimensions
      elevation = 7.0
      arc_radius(1) = 7.0
      arc_radius(2) = 7.0
      arc_length = acos(1-elevation/(arc_radius(1)+arc_radius(2)))
      arc_length = arc_length*180/pi
      fiber_width(1) = 0.176
      fiber_width(2) = 0.192
      fiber_width(3) = 0.200
      straight_length(1) = 10.
      arc_spread = (arc_radius(1)+arc_radius(2))*sin(arc_length*pi/180)
      straight_length(2) = 180 - arc_spread - straight_length(1)
      print *, 'arc length of each curved region is',arc_length,' deg'
      print *, 'arc spread of the complete S-bend is',arc_spread,' cm'

c define the tracking media
      isvol = 0
      ifield = 0
      fieldm = 0
      tmaxfd = 0
      stemax = 0
      deemax = 0
      epsil = 1e-5
      stmin = 1e-5
      nwbuf = 0
      nmed = 1
      natmed = 'atmosphere medium'
      nmat = 15 ! air
      call gstmed(nmed,natmed,nmat,isvol,ifield,fieldm,tmaxfd,stemax,
     +            deemax,epsil,stmin,ubuf,nwbuf)
      nmed = 100
      natmed = 'fiber core medium'
      nmat = 100
      call gstmed(nmed,natmed,nmat,isvol,ifield,fieldm,tmaxfd,stemax,
     +            deemax,epsil,stmin,ubuf,nwbuf)
      nmed = 101
      natmed = 'inner clad medium'
      nmat = 101
      call gstmed(nmed,natmed,nmat,isvol,ifield,fieldm,tmaxfd,stemax,
     +            deemax,epsil,stmin,ubuf,nwbuf)
      nmed = 102
      natmed = 'outer clad medium'
      nmat = 102
      call gstmed(nmed,natmed,nmat,isvol,ifield,fieldm,tmaxfd,stemax,
     +            deemax,epsil,stmin,ubuf,nwbuf)

c assign the refractive properties of these media
      nmed = 1
      npckov = 2
      ppckov(1) = 1e-9
      ppckov(2) = 4e-9
      absco(1) = 1e6
      absco(2) = 1e6
      effic(1) = 0
      effic(2) = 0
#if WITH_PAINT
      rindex(1) = 2.5
      rindex(2) = 2.5
#else
      rindex(1) = 1
      rindex(2) = 1
#endif
      call gsckov(nmed,npckov,ppckov,absco,effic,rindex)
      nmed = 100
      rindex(1) = 1.60
      rindex(2) = 1.60
      absco(1) = 1000
      absco(2) = 1000
      call gsckov(nmed,npckov,ppckov,absco,effic,rindex)
      nmed = 101
      rindex(1) = 1.49
      rindex(2) = 1.49
      absco(1) = 1000
      absco(2) = 1000
      call gsckov(nmed,npckov,ppckov,absco,effic,rindex)
      nmed = 102
      rindex(1) = 1.42
      rindex(2) = 1.42
      absco(1) = 1000
      absco(2) = 1000
      call gsckov(nmed,npckov,ppckov,absco,effic,rindex)
      
c define the world box
      nmed = 1
      chname = 'SITE'
      chshape = 'BOX'
      par(1) = 300.
      par(2) = 300.
      par(3) = 300.
      npar = 3
      call gsvolu(chname,chshape,nmed,par,npar,ivolu)

c define the first straight section
      nmed = 100
      chname = 'COS1'
      chshape = 'BOX'
      par(1) = straight_length(1)/2
      par(2) = fiber_width(1)/2
      par(3) = fiber_width(1)/2
      npar = 3
      call gsvolu(chname,chshape,nmed,par,npar,ivolu)
      nmed = 101
      chname = 'ICS1'
      chshape = 'BOX'
      par(1) = straight_length(1)/2
      par(2) = fiber_width(2)/2
      par(3) = fiber_width(2)/2
      npar = 3
      call gsvolu(chname,chshape,nmed,par,npar,ivolu)
      nmed = 102
      chname = 'OCS1'
      chshape = 'BOX'
      par(1) = straight_length(1)/2
      par(2) = fiber_width(3)/2
      par(3) = fiber_width(3)/2
      npar = 3
      call gsvolu(chname,chshape,nmed,par,npar,ivolu)

c place the components of the first straight section
      chname = 'COS1'
      chmoth = 'ICS1'
      nr = 1
      x = 0
      y = 0
      z = 0
      irot = 0
      chonly = 'ONLY'
      call gspos(chname,nr,chmoth,x,y,z,irot,chonly)
      chname = 'ICS1'
      chmoth = 'OCS1'
      nr = 1
      x = 0
      y = 0
      z = 0
      irot = 0
      chonly = 'ONLY'
      call gspos(chname,nr,chmoth,x,y,z,irot,chonly)
      chname = 'OCS1'
      chmoth = 'SITE'
      nr = 1
      x = straight_length(1)/2
      y = 0
      z = 0
      irot = 0
      chonly = 'ONLY'
      call gspos(chname,nr,chmoth,x,y,z,irot,chonly)

c define the first arc section
      nmed = 100
      chname = 'COA1'
      chshape = 'TUBS'
      par(1) = arc_radius(1) - fiber_width(1)/2
      par(2) = arc_radius(1) + fiber_width(1)/2
      par(3) = fiber_width(1)/2
      par(4) = 90 - arc_length
      par(5) = 90
      npar = 5
      call gsvolu(chname,chshape,nmed,par,npar,ivolu)
      nmed = 101
      chname = 'ICA1'
      chshape = 'TUBS'
      par(1) = arc_radius(1) - fiber_width(2)/2
      par(2) = arc_radius(1) + fiber_width(2)/2
      par(3) = fiber_width(2)/2
      par(4) = 90 - arc_length
      par(5) = 90
      npar = 5
      call gsvolu(chname,chshape,nmed,par,npar,ivolu)
      nmed = 102
      chname = 'OCA1'
      chshape = 'TUBS'
      par(1) = arc_radius(1) - fiber_width(3)/2
      par(2) = arc_radius(1) + fiber_width(3)/2
      par(3) = fiber_width(3)/2
      par(4) = 90 - arc_length
      par(5) = 90
      npar = 5
      call gsvolu(chname,chshape,nmed,par,npar,ivolu)

c place the components of the first arc section
      chname = 'COA1'
      chmoth = 'ICA1'
      nr = 1
      x = 0
      y = 0
      z = 0
      irot = 0
      chonly = 'ONLY'
      call gspos(chname,nr,chmoth,x,y,z,irot,chonly)
      chname = 'ICA1'
      chmoth = 'OCA1'
      x = 0
      y = 0
      z = 0
      irot = 0
      chonly = 'ONLY'
      call gspos(chname,nr,chmoth,x,y,z,irot,chonly)
      chname = 'OCA1'
      chmoth = 'SITE'
      nr = 1
      x = straight_length(1)
      y = -arc_radius(1)
      z = 0
      irot = 0
      chonly = 'ONLY'
      call gspos(chname,nr,chmoth,x,y,z,irot,chonly)

c define the second arc section
      nmed = 100
      chname = 'COA2'
      chshape = 'TUBS'
      par(1) = arc_radius(2) - fiber_width(1)/2
      par(2) = arc_radius(2) + fiber_width(1)/2
      par(3) = fiber_width(1)/2
      par(4) = 270 - arc_length
      par(5) = 270
      npar = 5
      call gsvolu(chname,chshape,nmed,par,npar,ivolu)
      nmed = 101
      chname = 'ICA2'
      chshape = 'TUBS'
      par(1) = arc_radius(2) - fiber_width(2)/2
      par(2) = arc_radius(2) + fiber_width(2)/2
      par(3) = fiber_width(2)/2
      par(4) = 270 - arc_length
      par(5) = 270
      npar = 5
      call gsvolu(chname,chshape,nmed,par,npar,ivolu)
      nmed = 102
      chname = 'OCA2'
      chshape = 'TUBS'
      par(1) = arc_radius(2) - fiber_width(3)/2
      par(2) = arc_radius(2) + fiber_width(3)/2
      par(3) = fiber_width(3)/2
      par(4) = 270 - arc_length
      par(5) = 270
      npar = 5
      call gsvolu(chname,chshape,nmed,par,npar,ivolu)

c place the components of the second arc section
      chname = 'COA2'
      chmoth = 'ICA2'
      nr = 1
      x = 0
      y = 0
      z = 0
      irot = 0
      chonly = 'ONLY'
      call gspos(chname,nr,chmoth,x,y,z,irot,chonly)
      chname = 'ICA2'
      chmoth = 'OCA2'
      nr = 1
      x = 0
      y = 0
      z = 0
      irot = 0
      chonly = 'ONLY'
      call gspos(chname,nr,chmoth,x,y,z,irot,chonly)
      chname = 'OCA2'
      chmoth = 'SITE'
      x = straight_length(1) + arc_spread
      y = arc_radius(2) - elevation
      z = 0
      irot = 0
      chonly = 'ONLY'
      call gspos(chname,nr,chmoth,x,y,z,irot,chonly)

c define the second straight section
      nmed = 100
      chname = 'COS2'
      chshape = 'BOX'
      par(1) = straight_length(2)/2
      par(2) = fiber_width(1)/2
      par(3) = fiber_width(1)/2
      npar = 3
      call gsvolu(chname,chshape,nmed,par,npar,ivolu)
      nmed = 101
      chname = 'ICS2'
      chshape = 'BOX'
      par(1) = straight_length(2)/2
      par(2) = fiber_width(2)/2
      par(3) = fiber_width(2)/2
      npar = 3
      call gsvolu(chname,chshape,nmed,par,npar,ivolu)
      nmed = 102
      chname = 'OCS2'
      chshape = 'BOX'
      par(1) = straight_length(2)/2
      par(2) = fiber_width(3)/2
      par(3) = fiber_width(3)/2
      npar = 3
      call gsvolu(chname,chshape,nmed,par,npar,ivolu)

c place the components of the second straight section
      chname = 'COS2'
      chmoth = 'ICS2'
      nr = 1
      x = 0
      y = 0
      z = 0
      irot = 0
      chonly = 'ONLY'
      call gspos(chname,nr,chmoth,x,y,z,irot,chonly)
      chname = 'ICS2'
      chmoth = 'OCS2'
      nr = 1
      x = 0
      y = 0
      z = 0
      irot = 0
      chonly = 'ONLY'
      call gspos(chname,nr,chmoth,x,y,z,irot,chonly)
      chname = 'OCS2'
      chmoth = 'SITE'
      nr = 1
      x = straight_length(1) + arc_spread + straight_length(2)/2
      y = -elevation
      z = 0
      irot = 0
      chonly = 'ONLY'
      call gspos(chname,nr,chmoth,x,y,z,irot,chonly)
      end

      real function guplsh(medi0,medi1)
      integer medi0,medi1
      if (medi0 .eq. 102 .and. medi1 .eq. 1) then
        guplsh = OUTER_SURFACE_QUALITY
      else if (medi0 .eq. 101 .and. medi1 .eq. 102) then
        guplsh = MIDDLE_SURFACE_QUALITY
      else if (medi0 .eq. 100 .and. medi1 .eq. 101) then
        guplsh = INNER_SURFACE_QUALITY
      else
        guplsh = 1.
      endif
      end