*
* wc.kumac - Geant++ macros to accept drawing origin arguments u,v
*            in world coordinates (cm) instead of the viewport coordinates
*            expected by the Geant3 drawing package commands.  Otherwise
*            the arguments are the same as those of the Geant++ commands.
*
* Usage :  exec wc#dcut site x .1 10 10 .01 .01
*
* The following drawing commands are supported:
*  draw - plot a view of a named volume in its own reference frame
*  dcut - plot a planar intersection with a named volume in its frame
*  dxcut - same as dcut, but plane is general instead of only xy, xz or yz
*  dvolume - same as draw, but can use master reference frame and select
*            particular instances of a given volume in the geometry tree
* The following graphical queries are supported
*  pick - clicks of the mouse are converted into 3d points using the
*         current drawing parameters and cut plane (see below).
* Multiple clicks of the mouse with the left button select the points,
* while the pick input is terminated with a click of the right button.
* Special arguments to the pick command produce special results:
*  pick m - pick two points followed by a right-click, and the
*           distance between the two points is computed and printed.
*  pick o - pick a single point, and that will be the new origin for
*           subsequence wc drawing commands.
* The following calls modify the behavior of the above macros
*  setorigin x y z - wc origin (drawing focus) is set to (x,y,z)
*  setview u0 v0 su sv theta phi psio
*      - sets location of drawing focus in graphics window to (u0,v0)
*      - sets scale factors for drawing in graphics window to (su,sv)
*      - sets viewing angles (degrees) to Geant angles theta,phi,psio
*        where psio differs from the Euler psi of the camera in that
*        psio=0 is defined as the psi for which the y axis projects to
*        the up direction on the graphics window, unless the projection
*        of yhat is zero, in which case xhat is chosen as up
*  setangles theta phi psio
*      - attempts to calculate the actual Euler psi corresponding to the
*        Geant convention for psio (degrees) and save them all as globals
*  setcutplane ctheta cphi xint
*      - defines the virtual cut plane in the graphics image as the
*        plane perpendicular to the ctheta,cphi (degrees) direction,
*        whose intercept with the axis ctheta,cphi is xint (cm).
*
* author: Richard Jones
* date: May 10, 2004
*
MACRO wc
  mess 'wc.kumac - Geant++ macros to accept drawing origin arguments u,v'
  mess '           in world coordinates (cm) instead of the viewport coordinates'
  mess '           expected by the Geant3 drawing package commands.  Otherwise'
  mess '           the arguments are the same as those of the Geant++ commands.'
  mess
  mess 'Usage :  exec wc#dcut site x .1 10 10 .01 .01'
  mess
  mess 'The following drawing commands are supported:'
  mess ' draw - plot a view of a named volume in its own reference frame'
  mess ' dcut - plot a planar intersection with a named volume in its frame'
  mess ' dxcut - same as dcut, but plane is general instead of only xy, xz or yz'
  mess ' dvolume - same as draw, but can use master reference frame and select'
  mess '           particular instances of a given volume in the geometry tree'
  mess 'The following graphical queries are supported'
  mess ' pick - clicks of the mouse are converted into 3d points using the'
  mess '        current drawing parameters and cut plane (see below).'
  mess ' Multiple clicks of the mouse with the left button select the points,'
  mess ' while the pick input is terminated with a click of the right button.'
  mess ' Special arguments to the pick command produce special results:'
  mess '   pick m - pick two points followed by a right-click, and the'
  mess '            distance between the two points is computed and printed.'
  mess '   pick o - pick a single point, and that will be the new origin for'
  mess '            subsequence wc drawing commands.'
  mess ' The following calls modify the behavior of the above macros'
  mess ' setorigin x y z - wc origin (drawing focus) is set to (x,y,z)'
  mess ' setview u0 v0 su sv theta phi psio'
  mess '     - sets location of drawing focus in graphics window to (u0,v0)'
  mess '     - sets scale factors for drawing in graphics window to (su,sv)'
  mess '     - sets viewing angles (degrees) to Geant angles theta,phi,psio'
  mess '       where psio differs from the Euler psi of the camera in that'
  mess '       psio=0 is defined as the psi for which the y axis projects to'
  mess '       the up direction on the graphics window, unless the projection'
  mess '       of yhat is zero, in which case xhat is chosen as up'
  mess ' setangles theta phi psio'
  mess '     - attempts to calculate the actual Euler psi corresponding to the'
  mess '       Geant convention for psio (degrees) and save them all as globals'
  mess ' setcutplane ctheta cphi xint'
  mess '     - defines the virtual cut plane in the graphics image as the'
  mess '       plane perpendicular to the ctheta,cphi (degrees) direction,'
  mess '       whose intercept with the axis ctheta,cphi is xint (cm).'
  mess
  mess 'author: Richard Jones'
  mess 'date: May 10, 2004'
  mess
RETURN

MACRO draw
  if ([1].eq.' ') then
    message "Error: first argument (volume name) is mandatory"
    exitm
  endif
  global/import wc*
  if ([2].ne.' ') then
    theta=[2]
  elseif ($defined(wctheta).eq.wctheta) then
    theta=[wctheta]
  else
    theta=30
  endif
  if ([3].ne.' ') then
    phi=[3]
  elseif ($defined(wcphi).eq.wcphi) then
    phi=[wcphi]
  else
    phi=30
  endif
  if ([4].ne.' ') then
    psio=[4]
  elseif ($defined(wcpsio).eq.wcpsio) then
    psio=[wcpsio]
  else
    psio=0
  endif
  if ([5].ne.' ') then
    u0=[5]
  elseif ($defined(wcu0).eq.wcu0) then
    u0=[wcu0]
  else
    u0=0
  endif
  if ([6].ne.' ') then
    v0=[6]
  elseif ($defined(wcv0).eq.wcv0) then
    v0=[wcv0]
  else
    v0=0
  endif
  if ([7].ne.' ') then
    su=[7]
  elseif ($defined(wcsu).eq.wcsu) then
    su=[wcsu]
  else
    su=1
  endif
  if ([8].ne.' ') then
    sv=[8]
  elseif ($defined(wcsv).eq.wcsv) then
    sv=[wcsv]
  else
    sv=1
  endif
  if ($defined(wcoriginx).ne.wcoriginx) then
    exec setorigin 0 0 0
  endif
  if (([su].eq.0).or.([sv].eq.0)) then
    message "Error: horizontal and vertical scale factors must be non-zero"
    exitm
  endif
  exec setangles [theta] [phi] [psio]
  global/import wc*
  sigma wcx0=([wcoriginx])*cos([wcphi]*pi/180)+([wcoriginy])*sin([wcphi]*pi/180)
  sigma wcy0=-([wcoriginx])*sin([wcphi]*pi/180)+([wcoriginy])*cos([wcphi]*pi/180)
  sigma wcz0=[wcoriginz]
  sigma wcz1=wcz0*cos([wctheta]*pi/180)+wcx0*sin([wctheta]*pi/180)
  sigma wcx1=-wcz0*sin([wctheta]*pi/180)+wcx0*cos([wctheta]*pi/180)
  sigma wcy1=wcy0
  sigma wcx2=wcx1*cos([wcpsi]*pi/180)+wcy1*sin([wcpsi]*pi/180)
  sigma wcy2=-wcx1*sin([wcpsi]*pi/180)+wcy1*cos([wcpsi]*pi/180)
  sigma wcz2=wcz1
  u=$sigma([u0]-(wcx2*[su]))
  v=$sigma([v0]-(wcy2*[sv]))
  draw [1] [theta] [phi] [psio] [u] [v] [su] [sv]
  exec setview [theta] [phi] [psio] [u0] [v0] [su] [sv]
  exec setcutplane [theta] [phi] 0
RETURN

MACRO dcut
  if ([1].eq.' ') then
    message "Error: first argument (volume name) is mandatory"
    exitm
  endif
  global/import wc*
  if ([2].ne.' ') then
    caxis=$lower([2])
  elseif ($defined(wcaxis).eq.wcaxis) then
    caxis=[wcaxis]
  else
    caxis=z
  endif
  if ([3].ne.' ') then
    cxing=[3]
  elseif ($defined(wcxing).eq.wcxing) then
    cxing=[wcxing]
  else
    cxing=0
  endif
  if ([4].ne.' ') then
    u0=[4]
  elseif ($defined(wcu0).eq.wcu0) then
    u0=[wcu0]
  else
    u0=10
  endif
  if ([5].ne.' ') then
    v0=[5]
  elseif ($defined(wcv0).eq.wcv0) then
    v0=[wcv0]
  else
    v0=10
  endif
  if ([6].ne.' ') then
    su=[6]
  elseif ($defined(wcsu).eq.wcsu) then
    su=[wcsu]
  else
    su=1
  endif
  if ([7].ne.' ') then
    sv=[7]
  elseif ($defined(wcsv).eq.wcsv) then
    sv=[wcsv]
  else
    sv=1
  endif
  if ($defined(wcoriginx).ne.wcoriginx) then
    exec setorigin 0 0 0
  endif
  if (([su].eq.0).or.([sv].eq.0)) then
    message "Error: horizontal and vertical scale factors must be non-zero"
    exitm
  endif
  global/import wc*
  if ([caxis].eq.x) then
    u=$sigma([u0]-([wcoriginz]*[su]))
    v=$sigma([v0]-([wcoriginy]*[sv]))
    h=$sigma([cxing]+([wcoriginx]))
    phi=0; theta=-90; psi=0
    cutthe=90; cutphi=0
  elseif ([caxis].eq.y) then
    u=$sigma([u0]-([wcoriginz]*[su]))
    v=$sigma([v0]-([wcoriginx]*[sv]))
    h=$sigma([cxing]+([wcoriginy]))
    phi=-90; theta=-90; psi=0
    cutthe=90; cutphi=90
  elseif ([caxis].eq.z) then
    u=$sigma([u0]-([wcoriginx]*[su]))
    v=$sigma([v0]-([wcoriginy]*[sv]))
    h=$sigma([cxing]+([wcoriginz]))
    phi=0; theta=0; psi=0
    cutthe=0; cutphi=0
  else
    message Error: argument 2=[caxis] must be either x, y or z
    exitm
  endif
  dcut [1] [caxis] [h] [u] [v] [su] [sv]
  global/create wcaxis [caxis] 'cut axis for wc'
  global/create wcxing [cxing] 'cut axis crossing value for wc'
  exec setview [theta] [phi] 0 [u0] [v0] [su] [sv]
  exec setcutplane [cutthe] [cutphi] [cxing]
RETURN

MACRO dxcut
  if ([1].eq.' ') then
    message "Error: first argument (volume name) is mandatory"
    exitm
  endif
  global/import wc*
  if ([2].ne.' ') then
    cutthe=[2]
  elseif ($defined(wcthecut).eq.wcthecut) then
    cutthe=[wcthecut]
  else
    message "Error: no default for argument CUTTHE"
    exitm
  endif
  if ([3].ne.' ') then
    cutphi=[3]
  elseif ($defined(wcphicut).eq.wcphicut) then
    cutphi=[wcphicut]
  else
    message "Error: no default for argument CUTPHI"
    exitm
  endif
  if ([4].ne.' ') then
    cutval=[4]
  elseif ($defined(wcvalcut).eq.wcvalcut) then
    cutval=[wcvalcut]
  else
    cutval=0
  endif
  if ([5].ne.' ') then
    theta=[5]
  elseif ($defined(wctheta).eq.wctheta) then
    theta=[wctheta]
  else
    theta=30
  endif
  if ([6].ne.' ') then
    phi=[6]
  elseif ($defined(wcphi).eq.wcphi) then
    phi=[wcphi]
  else
    phi=30
  endif
  if ([7].ne.' ') then
    u0=[7]
  elseif ($defined(wcu0).eq.wcu0) then
    u0=[wcu0]
  else
    u0=10
  endif
  if ([8].ne.' ') then
    v0=[8]
  elseif ($defined(wcv0).eq.wcv0) then
    v0=[wcv0]
  else
    v0=10
  endif
  if ([9].ne.' ') then
    su=[9]
  elseif ($defined(wcsu).eq.wcsu) then
    su=[wcsu]
  else
    su=1
  endif
  if ([10].ne.' ') then
    sv=[10]
  elseif ($defined(wcsv).eq.wcsv) then
    sv=[wcsv]
  else
    sv=1
  endif
  if ($defined(wcoriginx).ne.wcoriginx) then
    exec setorigin 0 0 0
  endif
  if (([su].eq.0).or.([sv].eq.0)) then
    message "Error: horizontal and vertical scale factors must be non-zero"
    exitm
  endif
  exec setangles [theta] [phi] 0
  global/import wc*
  sigma wcnx=sin([cutthe]*pi/180)*cos([cutphi]*pi/180)
  sigma wcny=sin([cutthe]*pi/180)*sin([cutphi]*pi/180)
  sigma wcnz=cos([cutthe]*pi/180)
  sigma wcx0=([wcoriginx])*cos([wcphi]*pi/180)+([wcoriginy])*sin([wcphi]*pi/180)
  sigma wcy0=-([wcoriginx])*sin([wcphi]*pi/180)+([wcoriginy])*cos([wcphi]*pi/180)
  sigma wcz0=[wcoriginz]
  sigma wcz1=wcz0*cos([wctheta]*pi/180)+wcx0*sin([wctheta]*pi/180)
  sigma wcx1=-wcz0*sin([wctheta]*pi/180)+wcx0*cos([wctheta]*pi/180)
  sigma wcy1=wcy0
  sigma wcx2=wcx1*cos([wcpsi]*pi/180)+wcy1*sin([wcpsi]*pi/180)
  sigma wcy2=-wcx1*sin([wcpsi]*pi/180)+wcy1*cos([wcpsi]*pi/180)
  sigma wcz2=wcz1
  u=$sigma([u0]-(wcx2*[su]))
  v=$sigma([v0]-(wcy2*[sv]))
  cut=$sigma([cutval]+([wcoriginx]*wcnx)+([wcoriginy]*wcny)+([wcoriginz]*wcnz))
  dxcut [1] [cutthe] [cutphi] [cut] [theta] [phi] [u] [v] [su] [sv]
  exec setview [theta] [phi] 0 [u0] [v0] [su] [sv]
  exec setcutplane [cutthe] [cutphi] [cutval]
RETURN

MACRO dvolume
  if ([1].eq.' ') then
    message "Error: first argument (n) is mandatory"
    exitm
  endif
  if ([2].eq.' ') then
    message "Error: second argument (volume list) is mandatory"
    exitm
  endif
  if ([3].eq.' ') then
    chnrs=[3]
  else
    chnrs=MARS
  endif
  global/import wc*
  if ([4].ne.' ') then
    theta=[4]
  elseif ($defined(wctheta).eq.wctheta) then
    theta=[wctheta]
  else
    theta=30
  endif
  if ([5].ne.' ') then
    phi=[5]
  elseif ($defined(wcphi).eq.wcphi) then
    phi=[wcphi]
  else
    phi=30
  endif
  if ([6].ne.' ') then
    psi=[6]
  elseif ($defined(wcpsi).eq.wcpsi) then
    psi=[wcpsi]
  else
    psi=30
  endif
  if ([7].ne.' ') then
    u0=[7]
  elseif ($defined(wcu0).eq.wcu0) then
    u0=[wcu0]
  else
    u0=10
  endif
  if ([8].ne.' ') then
    v0=[8]
  elseif ($defined(wcv0).eq.wcv0) then
    v0=[wcv0]
  else
    v0=10
  endif
  if ([9].ne.' ') then
    su=[9]
  elseif ($defined(wcsu).eq.wcsu) then
    su=[wcsu]
  else
    su=1
  endif
  if ([10].ne.' ') then
    sv=[10]
  elseif ($defined(wcsv).eq.wcsv) then
    sv=[wcsv]
  else
    sv=1
  endif
  if ($defined(wcoriginx).ne.wcoriginx) then
    exec setorigin 0 0 0
  endif
  if (([su].eq.0).or.([sv].eq.0)) then
    message "Error: horizontal and vertical scale factors must be non-zero"
    exitm
  endif
  exec setangles [theta] [phi] [psi]
  global/import wc*
  sigma wcx0=([wcoriginx])*cos([wcphi]*pi/180)+([wcoriginy])*sin([wcphi]*pi/180)
  sigma wcy0=-([wcoriginx])*sin([wcphi]*pi/180)+([wcoriginy])*cos([wcphi]*pi/180)
  sigma wcz0=[wcoriginz]
  sigma wcz1=wcz0*cos([wctheta]*pi/180)+wcx0*sin([wctheta]*pi/180)
  sigma wcx1=-wcz0*sin([wctheta]*pi/180)+wcx0*cos([wctheta]*pi/180)
  sigma wcy1=wcy0
  sigma wcx2=wcx1*cos([wcpsi]*pi/180)+wcy1*sin([wcpsi]*pi/180)
  sigma wcy2=-wcx1*sin([wcpsi]*pi/180)+wcy1*cos([wcpsi]*pi/180)
  sigma wcz2=wcz1
  u=$sigma([u0]-(wcx2*[su]))
  v=$sigma([v0]-(wcy2*[sv]))
  dvolume [1] [2] [chnrs] [theta] [phi] [psi] [u] [v] [su] [sv]
  exec setview [theta] [phi] [psi] [u0] [v0] [su] [sv]
  exec setcutplane [theta] [phi] 0
RETURN

MACRO pick opt
  global/import wc*
  if (($defined(wcu0).ne.wcu0).or. _
      ($defined(wcv0).ne.wcv0).or. _
      ($defined(wcsu).ne.wcsu).or. _
      ($defined(wcsv).ne.wcsv).or. _
      ($defined(wcthecut).ne.wcthecut).or. _
      ($defined(wcphicut).ne.wcphicut).or. _
      ($defined(wcvalcut).ne.wcvalcut)) then
    message "Error: you must use wc#xxx to create a drawing first"
    exitm
  endif
  sigma wcnx=sin([wcthecut]*pi/180)*cos([wcphicut]*pi/180)
  sigma wcny=sin([wcthecut]*pi/180)*sin([wcphicut]*pi/180)
  sigma wcnz=cos([wcthecut]*pi/180)
  sigma wcn1=wcnx*cos([wcphi]*pi/180)+wcny*sin([wcphi]*pi/180)
  sigma wcn2=-wcnx*sin([wcphi]*pi/180)+wcny*cos([wcphi]*pi/180)
  sigma wcn3=wcnz
  sigma wcn6=wcn3*cos([wctheta]*pi/180)+wcn1*sin([wctheta]*pi/180)
  sigma wcn4=-wcn3*sin([wctheta]*pi/180)+wcn1*cos([wctheta]*pi/180)
  sigma wcn5=wcn2
  sigma wcnu=wcn4*cos([wcpsi]*pi/180)+wcn5*sin([wcpsi]*pi/180)
  sigma wcnv=-wcn4*sin([wcpsi]*pi/180)+wcn5*cos([wcpsi]*pi/180)
  sigma wcnw=wcn6
  vec/del wcvec*
  vlocate wcvecu wcvecv
  if ($vexist(wcvecu).eq.0) then
    message '  no points entered'
    exitm
  else
    nhit=$vdim(wcvecu,1)
  endif
  sigma wcx0=(wcvecu-[wcu0])/[wcsu]
  sigma wcy0=(wcvecv-[wcv0])/[wcsv]
  sigma wcz0=([wcvalcut]-(wcx0*wcnu+wcy0*wcnv))/wcnw
  sigma wcx1=wcx0*cos([wcpsi]*pi/180)-wcy0*sin([wcpsi]*pi/180)
  sigma wcy1=wcx0*sin([wcpsi]*pi/180)+wcy0*cos([wcpsi]*pi/180)
  sigma wcz1=wcz0
  sigma wcz2=wcz1*cos([wctheta]*pi/180)-wcx1*sin([wctheta]*pi/180)
  sigma wcx2=wcz1*sin([wctheta]*pi/180)+wcx1*cos([wctheta]*pi/180)
  sigma wcy2=wcy1
  sigma wcx3=wcx2*cos([wcphi]*pi/180)-wcy2*sin([wcphi]*pi/180)
  sigma wcy3=wcx2*sin([wcphi]*pi/180)+wcy2*cos([wcphi]*pi/180)
  sigma wcz3=wcz2
  sigma wcx=[wcoriginx]+wcx3
  sigma wcy=[wcoriginy]+wcy3
  sigma wcz=[wcoriginz]+wcz3
  if ([opt].eq.o) then
    exec setorigin $sigma(wcx(1)) $sigma(wcy(1)) $sigma(wcz(1))
    message wc origin reset to ([wcoriginx],[wcoriginy],[wcoriginz])
  elseif ([opt].eq.m) then
    sigma d=sqrt((wcx([nhit])-wcx(1))**2_
                +(wcy([nhit])-wcy(1))**2_
                +(wcz([nhit])-wcz(1))**2)
    message measured distance $sigma(d) cm
  elseif ([opt].eq.l) then
    line $sigma(wcvecu(1)) $sigma(wcvecv(1)) $sigma(wcvecu(2)) $sigma(wcvecv(2))
  elseif ($fexist(wc.f).gt.0) then
    call wc.f77
    call wcprint([nhit],wcx,wcy,wcz)
  else
    do i=1,[nhit]
      message '   point' [i]: $sigma(wcx([i])) $sigma(wcy([i])) $sigma(wcz([i]))
    enddo
    message 'To print volume and medium info you need to copy wc.f to your working directory'
  endif
RETURN

MACRO setorigin 1=0 2=0 3=0
  global/create wcoriginx [1] 'x of wc origin in current DRS (cm)'
  global/create wcoriginy [2] 'y of wc origin in current DRS (cm)'
  global/create wcoriginz [3] 'z of wc origin in current DRS (cm)'
RETURN

MACRO setview theta phi psio u0 v0 su sv
  global/create wcu0 [u0] 'u offset of wc origin in drawing coordinates'
  global/create wcv0 [v0] 'v offset of wc origin in drawing coordinates'
  global/create wcsu [su] 'x scale factor of drawing in wc'
  global/create wcsv [sv] 'y scale factor of drawing in wc'
  exec setangles [theta] [phi] [psio]
RETURN

MACRO setcutplane theta phi w
  global/create wcthecut [theta] 'euler theta of normal to cut plane'
  global/create wcphicut [phi] 'euler phi of normal to cut plane'
  global/create wcvalcut [w] 'w-intercept of cut plane (for picking)'
RETURN

MACRO setangles theta phi psio
  global/import wc*
  sigma tanphi=tan([phi]*pi/180)
  sigma costhe=cos([theta]*pi/180)
  sigma tanpsi=-costhe*tanphi
  sigma cospsi=1/sqrt(1+tanpsi**2)
  sigma cosphi=cos([phi]*pi/180)
  sigma sinphi=sin([phi]*pi/180)
  if ($sigma(abs(costhe)).lt.0.005) then
    if ($sigma(abs(cosphi)).lt.0.005) then
      sigma psi=pi
      theta=90
      phi=90
    elseif ($sigma(cosphi).ge.0) then
      sigma psi=0
    else
      sigma psi=pi
    endif
  elseif ($sigma(cosphi*cospsi).ge.0) then
    sigma psi=atan(tanpsi)
  else
    sigma psi=atan(tanpsi)+pi
  endif
  psi=$sigma((psi*180/pi)-([psio]))
  global/create wcpsi [psi] 'euler psi of view in wc frame'
  global/create wctheta [theta] 'euler theta of view in wc frame'
  global/create wcphi [phi] 'euler phi of view in wc frame'
  global/create wcpsio [psio] 'euler psi offset of view in wc frame'
RETURN