From ishisaki@mpuax5.phys.metro-u.ac.jp Fri Nov  5 18:58:48 1999
Received: from mpuax5.phys.metro-u.ac.jp (mpuax5.phys.metro-u.ac.jp [133.86.36.53])
	by olegacy.gsfc.nasa.gov (8.9.1a/8.9.1) with ESMTP id SAA01460
	for <ebisawa@olegacy.gsfc.nasa.gov>; Fri, 5 Nov 1999 18:58:42 -0500 (EST)
Received: (from ishisaki@localhost)
	by mpuax5.phys.metro-u.ac.jp (8.8.8/3.6W-98050602) id IAA26486;
	Sat, 6 Nov 1999 08:58:09 +0900 (JST)
From: Yoshitaka ISHISAKI <ishisaki@mpuax5.phys.metro-u.ac.jp>
Message-Id: <199911052358.IAA26486@mpuax5.phys.metro-u.ac.jp>
To: ebisawa@olegacy.gsfc.nasa.gov
CC: gis@amalthea.phys.s.u-tokyo.ac.jp
In-reply-to: <199911052241.RAA30161@olegacy.gsfc.nasa.gov> (message from Ken
 Ebisawa on Fri, 5 Nov 1999 17:41:08 -0500 (EST))
Date: 6 Nov 1999 08:55:20 +0900
Subject: RE: new XRT files
Reply-to: ishisaki@phys.metro-u.ac.jp
Status: RO

$B3$O7Bt$5$^(B

>> On Fri, 5 Nov 1999 17:41:08 -0500 (EST),
>> Ken Ebisawa <ebisawa@olegacy.gsfc.nasa.gov> said:

> http://www-x.phys.metro-u.ac.jp/~ishisaki/crab990210/Report.html
> $B$K$h$j$^$9$H!"(BGIS$B%A!<%`$,(Btune$B$7$?(Bray-tracing$B$O!"(Bv6.2t1$B$H$$$&$b$N$G$9$M!#(B
> $B$^$@$3$l$G(B xrt_ea_v2_0.fits,xrt_psf_v2_0.fits$B$KBP1~$9$k(Bcalibration file$B$O(B
> $B:n@.$7$F$$$J$$!"$HM}2r$7$F$$$^$9!#(B

> $B@P:j$/$s!"$3$l$r:n$C$FD:$/$+!"(Bv6.2t1$B$N(Bray-tracing$B$r8x3+$7$FD:$1$^$9$+!)(B
> $B8x3+$7$FD:$1$l$P!"$3$A$i$G(BXRT cal file$B$O:n@.$G$-$k$H;W$$$^$9!#(B

$BD9$$4V$[$C$?$i$+$7$K$7$F$$$F$9$_$^$;$s!#(B
$B$H$j$"$($:!"(Bv6.2t1 $B$N(B ray-tracing $B$r$*$o$?$7$7$^$9!#(B
XRT cal file $B$,$=$A$i$G:n$l$k$H$$$&$3$H$G$"$l$P!"$*4j$$$7$^$9!#(B

---
$B@P:j(B $B6U>0(B (ISHISAKI Yoshitaka)
$BEl5~ETN)Bg3X(B $BM}3X8&5f2J(B    TEL: 0426-77-2485  FAX: 0426-77-2483
E-Mail: ishisaki@phys.metro-u.ac.jp


raycore_v6.2_t1.f:
c--------1---------2---------3---------4---------5---------6---------7--
c
c ASCA XRT raytracing
c
c 94/10/27 Y.Ishisaki version 1.0
c     modified from NAGOYA version for SimASCA
c 94/03/27 A.Furuzawa version 2.0
c     introduced so called horn (tsuno)
c 95/10/20 M.Watanabe version 2.1
c     Floating error at gold edge is fixed
c 95/07/03 version 3.0
c     Crab tune sono2
c 95/09/30 version 4.0
c     Image tuning by GRO***
c       nvfluct, etc
c 95/10/27 Y.Ishisaki version 4.0a
c     for Infinit Loop
c       write(*,*) variables in Raytrace_raycha() after theeta  = - dasin(st)
c       theeta = -theeta in Raytrace_scat() if ( theeta .lt. 0 )
c 95/10/27 M.Watanabe version 5.0
c     Crab tune sono3
c 96/11/21 A. Furuzawa version 6.0
c        made a common module(This code) used in the Nagoya Univ. version
c        and SimASCA
c 96/01/21 A. Furuzawa version 6.1
c     Infinite loop bug fixed
c 97/10/7 Ishisaki & Watanabe version 6.1_t1
c     M IV -36eV, M V -34eV shift
c 97/10/10 Ishisaki & Watanabe versiton 6.1_t1_1
c     Energy shift of optical const. in V6.1_t1 is limited at Au.
c 97/10/10 Ishisaki & Watanabe version 6.1_t2
c     1.5-2.2keV optical constant value approach the new henke value.
c 99/02/03 Y.ISHISAKI version 6.2_t1
c     M IV -25eV, M V -25eV shift
c
c=======================================================================
      Subroutine Raytrace_datred(iunt, iflct)
c     read astro-d xrt data
c  input   nothing
c  output  : ** in common /struct/**
c     real*8 r(120,1)    :   radius of primary top
c     real*8 r(120,2)    :   radius of primary bottom
c     real*8 r(120,3)    :   radius of secondary top
c     real*8 r(120,4)    :   radius of secondary bottom
c     real*8 beta(120,1) :   angle of primary mirror
c     real*8 beta(120,2) :   angle of secondary mirror
c     real*8 z(4)        :   distance from focal plane
c=======================================================================
      Implicit None
c input
      Integer iunt, iflct
c common
      Real*8 r, beta, z
      Common /Raytrace_struct/ r(120,4), beta(120,2), z(4)
c local
      Integer i, j, k
c function
      Real*8 Raytrace_ran
c
      Do i=1, 120
         Read(iunt,*,end=999) (r(i,j),j=1,4),(beta(i,k),k=1,2)
         If (iflct.eq.1) then
            Do j=1, 4
               r(i,j) = r(i,j)+((Raytrace_ran()-0.5d0)*0.035d0)
            End do
         End if
         beta(i,1) = (r(i,1)-r(i,2))/101.596
         beta(i,2) = (r(i,3)-r(i,4))/101.596
      End do
 999  Continue
c
      z(1) = 3605.3423525d0
      z(2) = 3503.7463525d0
      z(3) = 3496.2536475d0
      z(4) = 3394.6576475d0
c
      Return
      End

c=====================================================================
      Subroutine Raytrace_raytra(eng,angin,angphi,xp,yp,fai,rad,
     &                           wei,xf,yf,icond,ihitp,ihits)
c
c     << ray tracing routine in astro-d mirror >>
c     input  : eng      / energy of incident photon (keV)
c              angin,angphi / direction of incident photon (degree)
c              xp,yp    / position of incident photon at mirror top (mm)
c              fai,rad  / position of incident photon at mirror top (radian)
c     input/output:
c              wei      / photon weight, 0 for discard mode /
c     output : xf,yf    / position of reflected photon at focal plane
c              icond    / condition code
c              icond =   0 : normal end
c                    = 910 : hidden by primary housing
c                    = 911 : hidden by primary housing side wall
c                    = 915 : hidden by secondary housing
c                    = 916 : hidden by secondary housing side wall
c                    = 920 : hidden by 1st alignment bar
c                    = 921 : hidden by 2nd alignment bar
c                    = 923 : hidden by 3rd alignment bar
c                    = 924 : hidden by 4th alignment bar
c                    = 925 : hidden by 1&14 sector
c                    = 930 : hidden by primary foil
c                    = 935 : hidden by secondary foil
c                    = 940 : absorbed by primary inner mirror
c                    = 945 : absorbed by secondarr inner mirror
c                    = 950 : absorbed on primary foil
c                    = 955 : absorbed on secondary foil
c                    = 960 : only primary reflect
c                    = 970 : only secondary reflect
c                    = 980 : no reflect
c                    = 990 : abnormal reflected
c
c=====================================================================
      Implicit None
c input
      Real*8 eng, angin, angphi, xp, yp, fai, rad
c input/output
      Real*8 wei
c output
      Real*8 xf, yf
      Integer icond, ihitp, ihits
c common
      Real*8 r, beta, z
      Common /Raytrace_struct/ r(120,4),beta(120,2),z(4)
c constant
      Integer mskhs
      Parameter (mskhs=1)
      Real*8 M_PI, M_90DEG, DEG2RAD, WALBAR, SECTOR
      Parameter (M_PI = 3.1415926535897932385d0)
      Parameter (M_90DEG = M_PI/2d0)
      Parameter (DEG2RAD = 0.017453292519943295769d0)
      Parameter (WALBAR = 1.574738d0)
      Parameter (SECTOR = 6.428d0*DEG2RAD)
c local
      Integer ihit, ihitm, km, knest, inhc, irc
      Real*8 theeta, phi
      Real*8 ray0(3), rayn(3), ray0s(2)
      Real*8 t, sint, widfai, diffai, secfai, rs
      Real*8 acref, aref, zcref, zref
c function
      Real*8 Raytrace_ran, Raytrace_crefrd, Raytrace_refrd

c-------------------------- initializasiot -----------------------------
      icond    = 0
      ray0(1)  = xp
      ray0(2)  = yp
      ray0(3)  = z(1)
      ihitp = 0
      ihits = 0

      theeta   = angin  * DEG2RAD
      phi      = angphi * DEG2RAD

      sint = dsin(theeta)
      rayn(1)  = sint * dcos(phi)
      rayn(2)  = sint * dsin(phi)
      rayn(3)  = - dcos(theeta)

      Do km=1, 2
        If (km.eq.2) then
          t = -( ray0(3) - z(3) ) / rayn(3)
          ray0(1) = ray0(1) + t * rayn(1)
          ray0(2) = ray0(2) + t * rayn(2)
          ray0(3) = z(3)
          rad     = ray0(1)*ray0(1)+ray0(2)*ray0(2)
          rad     = dsqrt(rad)
          fai     = datan(ray0(2)/ray0(1))
          If (ray0(1).lt.0.0d0) then
            fai   = fai + M_PI
          Else if (ray0(2).lt.0.0d0) then
            fai   = fai + 2.0d0 * M_PI
          End if
        End if
        ray0s(1) = ray0(1)
        ray0s(2) = ray0(2)
c------------- finding the hiddeness by housing ------------------------
        If ( km.eq.1 ) then
           If ((dabs(ray0(1)).lt.5.0036d0)) goto 910
           If ((dabs(ray0(2)).lt.5.0036d0)) goto 910
        Else
           If ((dabs(ray0(1)).lt.5.0036d0)) goto 915
           If ((dabs(ray0(2)).lt.5.0036d0)) goto 915
        End if
c------------- finding the hiddeness by alignment bar ------------------
        widfai = WALBAR / rad
        diffai = dmod(fai+widfai/2.0d0, SECTOR)
        If ( diffai.lt.widfai ) then
           If ( km.eq.1 ) goto 920
           If ( km.eq.2 ) goto 923
        End if
c------------- finding the photon in 1&14 sector -----------------------
        If ((mskhs.eq.1).and.(km.eq.1)) then
          secfai = dmod(fai, M_90DEG)
          if (secfai.lt. SECTOR) goto 925
          if (secfai.gt. M_90DEG - SECTOR) goto 925
        End if
c------------------ finding the hiddeness by foils ---------------------
        If (km.eq.1) then
           rs = 58.9903745d0
        Else
           rs = 57.6391477d0
        End if
        Do knest=1, 120
           If ((rad.gt.rs).and.(rad.lt.r(knest,km*2-1))) goto 22
           rs = r(knest,km*2-1) + 0.15d0
        End do
        if (km.eq.1) goto 930
        if (km.eq.2) goto 935
 22     Continue
c------------------- ray tracing in knest'th space ---------------------
        ihitm = 1
        Do ihit=1, 100
          Call Raytrace_hinm(theeta,ray0,rayn,knest,km,inhc)
          If (inhc.eq.1) then
            zcref = Raytrace_crefrd(dabs(theeta), eng)
            If ( wei.gt.0d0 ) then
               wei = wei * zcref
               If ( km.eq.1 ) then
                  ihitp = ihitp + ihitm
               Else
                  ihits = ihits + ihitm
               End if
            Else
               acref = Raytrace_ran()
               If ( km.eq.1 ) then
                  If ( zcref.lt.acref ) goto 940
                  ihitp = ihitp + ihitm
               Else
                  If ( zcref.lt.acref ) goto 945
                  ihits = ihits + ihitm
               End if
            End if
          Endif
          If ((inhc.eq.0).and.(ihit.ne.1)) goto 24
c          write(*,*) 'O.K.3'
c          Call Raytrace_raycha(eng,theeta,ray0,rayn,
c     &                         rad,fai,knest,km,irc)
          Call Raytrace_raycha(eng,theeta,ray0,rayn,
     &                         rad,fai,knest,km,irc,ihitp,ihits,ihit)
c          write(*,*) 'O.K.6'
          ihitm = ihitm * 2
          If (irc.eq.1) then
            zref = Raytrace_refrd(dabs(theeta), eng)
            If ( wei.gt.0d0 ) then
               wei = wei * zref
               If ( km.eq.1 ) then
                  ihitp = ihitp + ihitm
               Else
                  ihits = ihits + ihitm
               End if
            Else
               aref = Raytrace_ran()
               If ( km.eq.1 ) then
                  If ( zref.lt.aref ) goto 950
                  ihitp = ihitp + ihitm
               Else
                  if ( zref.lt.aref ) goto 955
                  ihits = ihits + ihitm
               End if
            End if
          End if
          if (irc.eq.0) goto 24
          ihitm = ihitm * 2
        End do
 24     Continue
        t = -( ray0(3) - z(2*km) ) / rayn(3)
        ray0(1) = ray0(1) + t * rayn(1)
        ray0(2) = ray0(2) + t * rayn(2)
        ray0(3) = z(2*km)
        rad     = ray0(1)*ray0(1)+ray0(2)*ray0(2)
        rad     = dsqrt(rad)
        fai     = datan(ray0(2)/ray0(1))
        If (ray0(1).lt.0.0d0) then
          fai   = fai + M_PI
        Else if (ray0(2).lt.0.0d0) then
          fai   = fai + 2.0d0 * M_PI
        End if
c-------- finding the hiddeness by housing side wall -------------------
        If (km.eq.1) then
          if((dabs(ray0(1)).lt.5.0036d0).or.
     &       (dabs(ray0(2)).lt.5.0036d0).or.
     &       ((ray0(1)*ray0s(1)).lt.0.0d0).or.
     &       ((ray0(2)*ray0s(2)).lt.0.0d0)) goto 911
        Else
          if((dabs(ray0(1)).lt.5.0036d0).or.
     &       (dabs(ray0(2)).lt.5.0036d0).or.
     &       ((ray0(1)*ray0s(1)).lt.0.0d0).or.
     &       ((ray0(2)*ray0s(2)).lt.0.0d0)) goto 916
        Endif
c-------- finding the hiddeness by alignment bar -----------------------
        widfai = WALBAR / rad
        diffai = dmod(fai+widfai/2.0d0, SECTOR)
        If ( diffai.lt.widfai ) then
           If ( km.eq.1 ) goto 921
           if ( km.eq.2 ) goto 924
        End if
      End do
c---------------- calcuration of the focal plane position --------------
      t      = - ray0(3) / rayn(3)
      xf     = ray0(1) + t * rayn(1)
      yf     = ray0(2) + t * rayn(2)

      If (ihitp.eq.2) then
         If (ihits.eq.2) then
            icond = 0
            Return
         Else if (ihits.eq.0) then
            goto 960
         End if
      Else if (ihitp.eq.0) then
         If (ihits.eq.2) then
            goto 970
         Else if (ihits.eq.0) then
            goto 980
         End if
      End if

      goto 990
c----------------------- extraordinary end treatment -------------------
 910  icond   = 910
      wei     = 0d0
      xf      = 999d0
      yf      = 999d0
      return

 911  icond   = 911
      wei     = 0d0
      xf      = 999d0
      yf      = 999d0
      return

 915  icond   = 915
      wei     = 0d0
      xf      = 999d0
      yf      = 999d0
      return

 916  icond   = 916
      wei     = 0d0
      xf      = 999d0
      yf      = 999d0
      return

 920  icond   = 920
      wei     = 0d0
      xf      = 999d0
      yf      = 999d0
      return

 921  icond   = 921
      wei     = 0d0
      xf      = 999d0
      yf      = 999d0
      return

 923  icond   = 923
      wei     = 0d0
      xf      = 999d0
      yf      = 999d0
      return

 924  icond   = 924
      wei     = 0d0
      xf      = 999d0
      yf      = 999d0
      return

 925  icond   = 925
      wei     = 0d0
      xf      = 999d0
      yf      = 999d0
      return

 930  icond   = 930
      wei     = 0d0
      xf      = 999d0
      yf      = 999d0
      return

 935  icond   = 935
      wei     = 0d0
      xf      = 999d0
      yf      = 999d0
      return

 940  icond   = 940
      wei     = 0d0
      xf      = 999d0
      yf      = 999d0
      return

 945  icond   = 945
      wei     = 0d0
      xf      = 999d0
      yf      = 999d0
      return

 950  icond   = 950
      wei     = 0d0
      xf      = 999d0
      yf      = 999d0
      return

 955  icond   = 955
      wei     = 0d0
      xf      = 999d0
      yf      = 999d0
      return

 960  icond   = 960
      return

 970  icond   = 970
      return

 980  icond   = 980
      return

 990  icond   = 990
      return

      End

c=====================================================================
      Subroutine Raytrace_hinm(theeta,ray0,rayn,knest,km,icond)
c
c     << inner mirror hitting condition >>
c     input  : ray0   / start point of ray
c              rayn   / direction vector of ray
c              knest  / mirror nesting number
c              km     / mirror number
c     output : icond  / inner mirror hitting condition
c                       = 0 .not hit
c                       = 1 .hit
c
c====================================================================
      Implicit None
c input
      Integer knest, km
c input/output
      Real*8 theeta, ray0(3), rayn(3)
c output
      Integer icond
c common
      Real*8 r, beta, z
      Common /Raytrace_struct/ r(120,4),beta(120,2),z(4)
c local
      Real*8 tb, tzb, a, b, c, d
      Real*8 w1, w2, ww1, ww2, rh, vn1, vn2, vn3, st, cosb
      Real*8 wwt, wwc, xx, yy
      Real*8 xh1,yh1,zh1

      Real*8 xh2,yh2,zh2
c

      xh2 = ray0(1)
      yh2 = ray0(2)
      zh2 = ray0(3)

      icond   = 0

      if( knest.ne.1 ) then
         tb     = dtan(beta(knest-1,km))
ccc        brod   = gus(sigma,amean)
ccc        tb     = dtan(beta(knest-1,km)+brod)

         tzb = ( z(km*2-1) - ray0(3) )*tb - r(knest-1,km*2-1)-0.15d0
         a   = rayn(1)**2 + rayn(2)**2 - rayn(3)**2 * tb ** 2
         b   = ray0(1)*rayn(1) + ray0(2)*rayn(2) + tb*tzb*rayn(3)
         c   = ray0(1)**2 + ray0(2)**2 - tzb**2
         d   = b**2 - a*c
         If ( d.ge.1.d-3 ) then
            w1     = -b/a
            w2     = dsqrt(d)/a
            ww1    = w1 - w2
            ww2    = w1 + w2
            wwc    = (ww1 - 1.d-3)*(ww2 - 1.d-3)
            if (ww1.gt.1.d-3.and.ww2.gt.1.d-3) then
               wwt = min(ww1,ww2)
            elseif (wwc.le.0.d0) then
               wwt = max(ww1,ww2)
            else
               return
            endif

            xh1 = ray0(1) + wwt * rayn(1)
            yh1 = ray0(2) + wwt * rayn(2)
            zh1 = ray0(3) + wwt * rayn(3)

            xx  = xh1 * ray0(1)
            yy  = yh1 * ray0(2)

            If (zh1.gt.z(km*2-1) .or. zh1.lt.z(km*2)) then
               return
            else
               if (xx.lt.0.d0 .or. yy.lt.0.d0) then
                  return
               else
                  icond=1
                  ray0(1)=xh1
                  ray0(2)=yh1
                  ray0(3)=zh1
               endif
            endif

            If   ( icond.eq.1 ) then
               rh     = ray0(1)**2 + ray0(2)**2
               rh     = dsqrt(rh)
               cosb   = dcos(beta(knest,km)) / rh
               vn1    = cosb * ray0(1)
               vn2    = cosb * ray0(2)
               vn3    = - dsin(beta(knest,km))
ccc            vn1    = dcos(beta(knest,km)+brod) * ray0(1) / rh
ccc            vn2    = dcos(beta(knest,km)+brod) * ray0(2) / rh
ccc            vn3    = - dsin(beta(knest,km)+brod)

               st        = rayn(1)*vn1 + rayn(2)*vn2 + rayn(3)*vn3
               theeta    = - dasin(st)

c               if (theeta.le.0.d0) then
c                  write(6,*) 'theta less than 0 in hinm'
c                  write(6,*) 'habetsu-shiki= ',d
c                  write(6,*) ' b**2        = ',b**2
c                  write(6,*) ' a*c         = ',a*c
c                  write(6,*) 'XRT stage    = ',km
c                  write(6,*) 'radius       = ',rh
c                  write(6,*) 'initial position1= ',xh2
c                  write(6,*) 'initial position2= ',yh2
c                  write(6,*) 'initial position3= ',zh2
c                  write(6,*) 'reflect position1= ',ray0(1)
c                  write(6,*) 'reflect position2= ',ray0(2)
c                  write(6,*) 'reflect position3= ',ray0(3)
c                  write(6,*) 'ray vector 1 = ',rayn(1)
c                  write(6,*) 'ray vector 2 = ',rayn(2)
c                  write(6,*) 'ray vector 3 = ',rayn(3)
c                  write(6,*) 'norm vector1 = ',vn1
c                  write(6,*) 'norm vector2 = ',vn2
c                  write(6,*) 'norm vector3 = ',vn3
c                  write(6,*) 'nest number  = ',knest
c                  write(6,*) 'theta        = ', theeta
c               endif
               rayn(1)   = rayn(1) - 2d0*st*vn1
               rayn(2)   = rayn(2) - 2d0*st*vn2
               rayn(3)   = rayn(3) - 2d0*st*vn3
            End if
         Else
            icond   = 0
         End if
      Else
         icond   = 0
      End if

      Return
      End

c=====================================================================
c      Subroutine Raytrace_raycha(eng,theeta,ray0,rayn,
c     &                           rad,fai,knest,km,icond)
      Subroutine Raytrace_raycha(eng,theeta,ray0,rayn,
     &                           rad,fai,knest,km,icond,ihitp,ihits,
     &                           ihit)
c
c     << ray line changing for reflection >>
c     input  : theeta / incident angle for the mirror (radiun)
c              ray0   / start point of ray
c              rayn   / direction vector of ray
c              knest  / mirror nesting number
c              km     / mirror number
c     output : ray0   / start point of ray
c              rayn   / direction vector of ray
c              rad    / start point of ray at secondary top
c              fai    / start point of ray at secondary top
c              icond  / outer mirror hitting condition
c                       = 0 .not hit
c                       = 1 .hit
c
c====================================================================
      Implicit None
c input
      Real*8 eng
      Integer knest, km
c input/output
      Real*8 theeta, ray0(3), rayn(3)
c output
      Real*8 rad, fai
      Integer icond
c common
      Real*8 r, beta, z
      Common /Raytrace_struct/ r(120,4),beta(120,2),z(4)
      Real*8 sigma1, sigma2, amean, consct
      Common /Raytrace_broad/ sigma1, sigma2, amean, consct
c local
      Real*8 tb, tzb, a, b, c, d, brod, raythe
      Real*8 w1, w2, ww1, ww2, rh, vn1, vn2, vn3, st, cosb
      Real*8 wwt, wwc,xx,yy
      Real*8 xh1,yh1,zh1
      Real*8 rpph
c constants
      Real*8 DEG2RAD, SCT2EDGE1, SCT2EDGE2, SCT13EDGE1, SCT13EDGE2
      Parameter (DEG2RAD = 0.017453292519943295769d0)
      parameter (SCT2EDGE1 = 6.428d0)
      parameter (SCT2EDGE2 = 12.857d0)
      parameter (SCT13EDGE1 = 90.d0 - SCT2EDGE2)
      parameter (SCT13EDGE2 = 90.d0 - SCT2EDGE1)
c function
      Real*8 Raytrace_scat, Raytrace_flor, Raytrace_nvfluct
      Real*8 Raytrace_ran

c For test
      integer ihitp,ihits,ihit

      icond = 0

c
      tb     = dtan(beta(knest,km))
ccc      tb     = dtan(beta(knest,km)+brod)
      tzb    = ( z(km*2-1) - ray0(3) )*tb - r(knest,km*2-1)
      a      = rayn(1)**2 + rayn(2)**2 - rayn(3)**2 * tb ** 2
      b      = ray0(1)*rayn(1) + ray0(2)*rayn(2) + tb*tzb*rayn(3)
      c      = ray0(1)**2 + ray0(2)**2 - tzb**2
      d      = b**2 - a*c

      If ( d.ge.0.0 ) then
         w1     = -b/a
         w2     = dsqrt(d)/a
         ww1    = w1-w2
         ww2    = w1+w2
         wwc    = (ww1 - 1.d-3)*(ww2 - 1.d-3)
         if (ww1.gt.1.d-3.and.ww2.gt.1.d-3) then
            wwt = min(ww1,ww2)
c            write(6,*) 'wwt = ',wwt
         elseif (wwc.le.0.d0) then
            wwt = max(ww1,ww2)
c            write(6,*) 'wwt = ',wwt
         else
            return
         endif
         xh1 = ray0(1) + wwt * rayn(1)
         yh1 = ray0(2) + wwt * rayn(2)
         zh1 = ray0(3) + wwt * rayn(3)

         xx  = xh1 * ray0(1)
         yy  = yh1 * ray0(2)
         If (zh1.gt.z(km*2-1) .or. zh1.lt.z(km*2)) then
            return
         else
            if (xx.lt.0.d0 .or. yy.lt.0.d0) then
               return
            else
               icond=1
               ray0(1)=xh1
               ray0(2)=yh1
               ray0(3)=zh1
            endif
         endif

         brod = 0.0
         If ( icond.ne.0 ) then
c            write(*,*) 'O.K.6'
            rh      = ray0(1)**2 + ray0(2)**2
            rh      = dsqrt(rh)
            cosb    = - dcos(beta(knest,km)) / rh
            vn1     = cosb * ray0(1)
            vn2     = cosb * ray0(2)
            vn3     = dsin(beta(knest,km))
            st      = rayn(1)*vn1 + rayn(2)*vn2 + rayn(3)*vn3
            theeta  = - dasin(st)

c            if (theeta.le.0.d0) then
c               write(6,*) 'theta less than 0 in raycha'
c               write(6,*) 'habetsu-shiki= ',d
c               write(6,*) 'XRT stage    = ',km
c               write(6,*) 'History (Hit)= ',ihit
c               write(6,*) 'Histiry (Pri)= ',ihitp
c               write(6,*) '        (Sce)= ',ihits
c               write(6,*) 'radius       = ',rh
c               write(6,*) 'ray position1= ',ray0(1)
c               write(6,*) 'ray position2= ',ray0(2)
c               write(6,*) 'ray position3= ',ray0(3)
c               write(6,*) 'ray vector 1 = ',rayn(1)
c               write(6,*) 'ray vector 2 = ',rayn(2)
c               write(6,*) 'ray vector 3 = ',rayn(3)
c               write(6,*) 'norm vector1 = ',vn1
c               write(6,*) 'norm vector2 = ',vn2
c               write(6,*) 'norm vector3 = ',vn3
c               write(6,*) 'nest number  = ',knest
c            endif

            rayn(1) = rayn(1) - 2d0*st*vn1
            rayn(2) = rayn(2) - 2d0*st*vn2
            rayn(3) = rayn(3) - 2d0*st*vn3

            If ((consct.ge.1.0).or.(sigma1.gt.0.0)) then
 10            If (consct.ge.1d0) then
                  brod = Raytrace_scat(theeta,eng)
               End if
               If (sigma1.gt.0d0) then
                  If (sigma2.gt.0.0) then
                     rpph=datan(dabs(ray0(2)/ray0(1)))/DEG2RAD
                     If ((rpph.gt.SCT2EDGE1.and.rpph.lt.SCT2EDGE2)
     & .or.(rpph.gt.SCT13EDGE1.and.rpph.lt.SCT13EDGE2)) then
                        brod = brod + Raytrace_flor(sigma2)
                     Else
                        brod = brod + Raytrace_nvfluct()
                     End if
                  Else
                     brod = brod + Raytrace_nvfluct()
                  End if
               End if
               If ((theeta+brod).le.0d0) then
                  brod = 0d0
                  goto 10
               End if
               raythe  = dacos(-rayn(3))
               cosb = dsin(raythe+brod) / dsin(raythe)
               rayn(1) = rayn(1) * cosb
               rayn(2) = rayn(2) * cosb
               rayn(3) = -dcos(raythe+brod)
c               if (rayn(3).gt.0.d0) then
c                  write(6,*) 'Reverse ray'
c               endif
            End if
         End if
      Else
         icond = 0
      End if

      Return
      End

c=====================================================================
      Real*8 function Raytrace_refrd (ang, eng)
c
c     << reflectivity table read >>
c     input  :  ang     / incident angle (radian)
c     output :  refrd   / reflectivity
c
c=====================================================================
      Implicit None
c input
      Real*8 ang, eng
c constant
      Real*8 M_PI
      Parameter (M_PI = 3.1415926535897932385d0)
c local
      Real*8 energy /-999d0/, delta, beta, rough /0d0/
      Save energy, delta, beta, rough
      Real*8 fw1, fw2, b, w, w1, w2, a2, b2, w3, w4, sa
c function
      Real*8 Raytrace_clcdel, Raytrace_clcbet
c
      If ( eng.ne.energy ) then
         energy = eng
         delta = Raytrace_clcdel(79,eng)
         beta = Raytrace_clcbet(79,eng)
      End if
c
      w = ( 1d0 - delta )**2
      b = beta**2
      fw1  = w - b
      fw2  = 4d0 * w * b
      w2    = fw1 - dcos(ang)**2
      w1    = w2*w2 + fw2
      w = dsqrt(w1)
      a2    = 0.5d0 * ( w + w2 )
      b2    = 0.5d0 * ( w - w2 )
      sa = dsin(ang)
      w = a2 + b2 + sa**2
      b = 2.0*dsqrt(a2)*sa
      w3    = w - b
      w4    = w + b

      Raytrace_refrd = w3/w4
      If ( rough .gt. 0d0 ) then
        Raytrace_refrd = Raytrace_refrd *
     &                   exp(-(4d0*M_PI*sa*rough*eng/12.4d0)**2)
      End if


      return
      end

c=====================================================================
      Real*8 function Raytrace_crefrd (ang, eng)
c
c     << reflectivity table read >>
c     input  :  ang     / incident angle (radian)
c     output :  crefrd  / reflectivity
c
c=====================================================================
      Implicit None
c input
      Real*8 ang, eng
c local
      Real*8 energy /-999d0/, delta, beta
      Save energy, delta, beta
      Real*8 fw1, fw2, b, w, w1, w2, a2, b2, w3, w4, sa
c function
      Real*8 Raytrace_clcdel, Raytrace_clcbet
c
      If ( eng.ne.energy ) then
         energy = eng
         delta = Raytrace_clcdel(6,eng)
         beta = Raytrace_clcbet(6,eng)
      End if
c
      w = ( 1d0 - delta )**2
      b = beta**2
      fw1  = w - b
      fw2  = 4d0 * w * b
      w2    = fw1 - dcos(ang)**2
      w1    = w2*w2 + fw2
      w = dsqrt(w1)
      a2    = 0.5d0 * ( w + w2 )
      b2    = 0.5d0 * ( w - w2 )
      sa = dsin(ang)
      w = a2 + b2 + sa**2
      b = 2d0*dsqrt(a2)*sa
      w3    = w - b
      w4    = w + b

      Raytrace_crefrd = w3/w4

c      if (Raytrace_crefrd.ge.1.d0) then
c         write(6,*) 'Abnormal reflectivity', Raytrace_crefrd
c      endif

      return
      end

c==============================================================
       Real*8 function Raytrace_flor (gamma)
c                    input gamma(arcmin.)
c
c                    output flor(rad.)
c
c
c==============================================================
      Implicit None
c input
      Real*8 gamma
c constant
      Real*8 M_PI
      Parameter (M_PI = 3.1415926535897932385d0)
c function
      Real*8 Raytrace_ran
c local
      Real*8 gamma1, a, b, flor
c
      gamma1 = gamma / 60d0
      a = -2d0 / gamma1
      b =  2d0 / gamma1

      flor = dtan((datan(b) - datan(a)) * Raytrace_ran() + datan(a))
      flor = flor * gamma / 2d0

      flor = flor * M_PI / 1.8d2

      Raytrace_flor = flor

      Return
      End

c=====================================================================
      Real*8 function Raytrace_scat (theeta,eng)
c
c     << make composite distribution >>
c     input  :  theeta(rad.)
c               eng(kev)
c     output :  scat(rad.)
c
c=====================================================================
      Implicit None
c input
      Real*8 theeta, eng
c constant
      Real*8 DEG2RAD, MIN2RAD, RAD2MIN
      Parameter (DEG2RAD = 0.017453292519943295769d0)
      Parameter (MIN2RAD = DEG2RAD/60d0)
      Parameter (RAD2MIN = 1.d0/MIN2RAD)
c common
      Real*8 sigma1 ,sigma2, amean, consct
      Common /Raytrace_broad/  sigma1, sigma2, amean, consct
      Real*8 s3, c3
      Common /Raytrace_calsct/ s3, c3
c local
      Integer i
      Real*8 core, r1, r2, ar1, rthr, thm
      Real*8 pwave, apw, fscat, f, fc
c function
      Real*8 Raytrace_ran
c
      If ( theeta .lt. 0 ) then
         write(*,*) 'theeta=', theeta
         theeta = -theeta
      End if

      thm = theeta*RAD2MIN
c      athm = dabs(thm)
c---------- definition of core width ---
      core = 0.8840d0
c---------------------------------------
 20   continue
c---------------------------------------------------------------------
c     generation of random number
c--------------------------------------------------------------------
c      do 10 i=1,41
c          r1=i-21.
c Y.ISHISAKI 95/09/19
cYI      r1 = (SimASCA_DRNDTS()-0.5d0)*100d0
c      r1 = SimASCA_DRNDTS()*(50.d0+thm)-thm
c      r2 = SimASCA_DRNDTS()
      r1=(Raytrace_ran()-0.5)*100.
c      r1 = Raytrace_ran()*(50.d0+thm)-thm
      r2=Raytrace_ran()
      ar1 = dabs(r1)
      rthr = r1*MIN2RAD
c--------------------------------------------------------------
c      calculation of core distribution
c--------------------------------------------------------------
      If ((-rthr).gt.theeta) goto 20
      If (ar1.ge.2d0) then
         If (ar1.lt.10d0) then
            If (r2.gt.0.2d0) goto 20
         Else
            If (r2.gt.0.015d0) goto 20
         Endif
         f = dlog10(ar1)
         fc = 10**(-0.5456d0-1.71887d0*f-0.70664d0*f*f)
      Else
         f = r1/core
         fc=dexp(-f*f/2d0)
      End if
c--------------------------------------------------------------
c      calculation of scattering component
c--------------------------------------------------------------
      pwave = 0.50671d0*eng*dsin(theeta+rthr/2d0)*dsin(rthr)
      apw = dabs(pwave)
      fscat = 68.8797d0*(1.d0+s3)*eng**3*dsin(theeta)*dsin(theeta+rthr)
     &     * dsin(theeta+rthr)* dexp(-9329d0*(1.d0+c3)*apw)
      If ((consct.ge.2d0).and.(consct.lt.3d0)) then
         fscat = 0d0
      Else if (consct.ge.3d0) then
         fc    = 0d0
      End if
      f=(fc+fscat)*0.8d0
c--------------------------------------------------------------
c       test print of fscat
c--------------------------------------------------------------
c      write(*,1003)eng,thm,r1,r2,f
c 1003  format(1x,5f12.5)
c--------------------------------------------------------------
      if (r2.gt.f) goto 20
 10   continue

      Raytrace_scat = rthr

      Return
      End

c=======================================================================
      Real*8 function Raytrace_clcdel(jz,er0)
c
c
c***********************************************************************
      Implicit None
c input
      Integer jz
      Real*8 er0
c constant
      Real*8 M_PI
      Parameter (M_PI = 3.1415926535897932385d0)
c local
      Integer i, j, k, l, n, kz, izh
      Real*8 sum, x, zh, wl, clcdel, er
c function
      Real*8 Raytrace_energy_change
c
      Integer iz(7)
      Real*8 a(7), rho(7), e(22,7)
      Real*8 c, h, rcl, avn

      data c, h /2.99793d+10, 4.1354d-18/
      data rcl, avn /2.81785d-13, 6.02486d+23/

c   carbon's constant
      data iz(7),a(7),rho(7) /6, 12.011d0, 1.54d0/

c   nickel's constant
      data iz(4),a(4),rho(4) /28, 58.71d0, 6.92d0/

c   gold's constant
      data iz(6),a(6),rho(6) /79, 197.0d0, 18.88d0/

c   energy level of k,l(i),l(ii),.....
c   carbon
      data e(1,7),e(2,7),e(3,7),e(4,7)
     &     /0.2842d0, 0.0d0, 0.0d0, 0.0d0/

c   nickel
      data e(1,4),e(2,4),e(3,4),e(4,4)
     &     /8.333d0, 1.008d0, 0.872d0, 0.855d0/
      data e(5,4),e(6,4),e(7,4),e(8,4),e(9,4)
     &     /0.112d0, 0.068d0, 0.068d0, 0.004d0, 0.004d0/

c   gold
      data e(1,6),e(2,6),e(3,6),e(4,6)
     &   /80.72d0, 14.35d0, 13.73d0, 11.919d0/
      data e(5,6),e(6,6),e(7,6),e(8,6),e(9,6)
     &   /3.425d0, 3.148d0, 2.743d0, 2.291d0, 2.206d0/
      data e(10,6),e(11,6),e(12,6),e(13,6),e(14,6),e(15,6),e(16,6)
     &   /0.759d0, 0.644d0, 0.545d0, 0.352d0, 0.334d0, 0.086d0, 0.083d0/
      data e(17,6),e(18,6),e(19,6),e(20,6),e(21,6),e(22,6)
     &   /0.108d0, 0.072d0, 0.054d0, 0.003d0, 0.003d0, 0.001d0/
c
      er = er0
      if(jz.eq.79) er = Raytrace_energy_change(er0)
      do 10 i=1,7
        if(jz.eq.iz(i)) goto 20
 10   continue
      Raytrace_clcdel=0.0d0
      return
 20   continue
      sum=iz(i)
      kz=0
      do 40 n=1,7
        do 35 l=1,n
          do 30 j=l,l+1
            if(j.eq.1)goto 30
            k=(n-1)**2+l+j-2
            if(k.gt.22)goto 50
            if(e(k,i)) 40,40,25
 25         x=er/e(k,i)
            if(x.eq.1.d0) then
              sum = 0.d0
              goto 50
            else
              zh= dble(jmin0(2*(j-1),iz(i)-kz))
              sum=sum+zh*dlog(dabs(x**2-1.0d0))/(x**2)
              izh=int(zh)
              kz=kz+izh
            endif
            if(kz.ge.iz(i)) goto 50
 30       continue
 35     continue
 40   continue
 50   continue

      wl=h*c/er
      clcdel=(rcl*wl*wl*avn*rho(i)/(2.0*M_PI*a(i)))*sum
      if(clcdel.lt.0.0) clcdel=0.0
      if(jz.eq.79) then
         if(er.lt.0.152) then
           clcdel = 1.065040  * dexp(-24.29330 * er)
         else if((er.ge.0.152).and.(er.lt.0.199)) then
           clcdel = 0.418700  * dexp(-18.15499 * er)
         else if((er.ge.0.199).and.(er.lt.0.286)) then
           clcdel = 0.031777  * dexp( -5.19690 * er)
         else if((er.ge.0.286).and.(er.lt.0.980)) then
           clcdel = 0.0119859 * dexp( -1.78744 * er)
         else if((er.ge.0.980).and.(er.lt.2.145)) then
           clcdel = 0.008700  * dexp( -1.46040 * er)
         else if((er.ge.2.145).and.(er.lt.2.206)) then
           clcdel = -0.0133595  *(er-2.09461)*(er-2.09461) + 0.000412683
         else if((er.ge.2.206).and.(er.lt.2.291)) then
           clcdel = -0.0174020  *(er-2.23930)*(er-2.23930) + 0.000309502
         else if((er.ge.2.291).and.(er.lt.2.397)) then
           clcdel = -0.00349211 *(er-2.39680)*(er-2.39680) + 0.000327174
         else if((er.ge.2.397).and.(er.lt.2.743)) then
           clcdel = -0.000520486*(er-2.39680)*(er-2.39680) + 0.000327174
         else if((er.ge.2.743).and.(er.lt.2.800)) then
           clcdel = -0.00908291 *(er-2.78728)*(er-2.78728) + 0.000282599
         else if((er.ge.2.800).and.(er.lt.3.450)) then
           clcdel = -0.00012313 * er + 0.0006241
         endif
         if((er.ge.2.3).and.(er.lt.3.68)) then
            clcdel = clcdel*(1.+0.092*(er-1.9465)/1.7335)
         elseif((er.ge.3.68).and.(er.lt.5.32)) then
            clcdel = clcdel*(1.+0.092*(5.32-er)/1.64)
         else if(er.ge.5.32) then
            clcdel = clcdel*(1.+0.044422*(er-5.32)
     &                         -0.0016108*(er*er-5.32*5.32))
         endif
cccccccccccccccccccccc  added at 1995.7.3 (crab tune sono2)
                 if((er.ge.2.7).and.(er.lt.4.5)) then
            clcdel = clcdel*(-0.01638*er+1.0437)
         elseif((er.ge.4.5).and.(er.lt.5.8)) then
            clcdel = clcdel*(0.02361*er+0.8656)
         elseif((er.ge.5.8).and.(er.lt.7.2)) then
            clcdel = clcdel*(-0.02071*er+1.1221)
         elseif((er.ge.7.2).and.(er.le.10.0)) then
            clcdel = clcdel*(0.005406*er+0.9345)
         endif
ccccccccccccccccccccc  added by  1995.10.27 (crab tune sono3)
                 if((er.ge.5.58).and.(er.lt.6.40)) then
            clcdel = clcdel*(-0.02955*er+1.167)
         elseif((er.ge.6.40).and.(er.lt.6.89)) then
            clcdel = clcdel*(0.02655*er+0.8080)
         elseif((er.ge.6.89).and.(er.lt.7.94)) then
            clcdel = clcdel*(0.002545*er+0.9740)
         elseif((er.ge.7.94).and.(er.le.8.53)) then
            clcdel = clcdel*(-0.0225*er+1.173)
         elseif((er.ge.8.53).and.(er.le.10.0)) then
            clcdel = clcdel*(0.001786*er+0.9659)
         endif
ccccccccccccccccccccc  added at 1997.10.10 (edge tune)
         if((er.ge.1.425).and.(er.lt.1.897)) then
            clcdel = clcdel*(0.1108*er+0.8421)
         elseif((er.ge.1.897).and.(er.lt.2.163)) then
            clcdel = clcdel*1.066*exp(-(er-2.014)*(er-2.014)/0.5308/2.0)
         endif
cccccccccccccccccccccc
c         if((er.ge.2.7).and.(er.lt.3.6)) then
c            clcdel = clcdel * (1 + (er-2.7)*0.04)
c         elseif((er.ge.3.6).and.(er.lt.4.5)) then
c            clcdel = clcdel * (1 + (4.5-er)*0.04)
c         else if(er.ge.4.500) then
c           clcdel = clcdel * ( 1.+ 6.*(er-4.5)/ 350.)
c         endif
      endif

      Raytrace_clcdel = clcdel
      return
      end

c*******************************************************************
      Real*8 function Raytrace_clcbet(jz,e0)
c    clcbet:imaginary part of reflaction
c    jz:atomic number
c    e:energy (kev)
c*******************************************************************
      Implicit None
c input
      Integer jz
      Real*8 e0
c local
      Integer i, j
      Real*8 xmu, xml, el2, el, xlogmu, e
c
      Real*8 d(3,60,13), factor(13)
      integer z(13)/1,2,4,5,6,8,9,13,22,24,28,47,79/
c function
      Real*8 Raytrace_energy_change
c     carbon
      data factor(5)/0.1513d-7/
      data d(1,1,5),d(2,1,5),d(3,1,5)/0.1085d0,2*1.92d4/
      data d(1,2,5),d(2,2,5),d(3,2,5)/0.2122d0,2*4.50d3/
      data d(1,3,5),d(2,3,5),d(3,3,5)/0.2833d0,2303.0d0,58392.0d0/
      data d(1,4,5),d(2,4,5),d(3,4,5)/0.5728d0,2*9.9033d0/
      data d(1,5,5),d(2,5,5),d(3,5,5)/1.000d0,2*2.18d3/
      data d(1,6,5),d(2,6,5),d(3,6,5)/1.012d0,2*2.15d3/
      data d(1,7,5),d(2,7,5),d(3,7,5)/5.000,2*18.6d0/
      data d(1,8,5),d(2,8,5),d(3,8,5)/10.00d0,2*1.98d0/
      data d(1,9,5),d(2,9,5),d(3,9,5)/20.00d0,2*0.2141d0/
      data d(1,10,5),d(2,10,5),d(3,10,5)/40.00d0,2*0.03284d0/
      data d(1,11,5),d(2,11,5),d(3,11,5)/80.00d0,2*0.02031d0/
      data d(1,12,5),d(2,12,5),d(3,12,5)/100.00d0,2*0.02146d0/
c     nickel
      data factor(11)/0.8694d-7/
      data d(1,1,11),d(2,1,11),d(3,1,11)/0.100d0,2*8.00d4/
      data d(1,2,11),d(2,2,11),d(3,2,11)/0.150d0,2*5.70d4/
      data d(1,3,11),d(2,3,11),d(3,3,11)/0.200d0,2*3.60d4/
      data d(1,4,11),d(2,4,11),d(3,4,11)/0.853d0,2160.0d0,7800.0d0/
      data d(1,5,11),d(2,5,11),d(3,5,11)/0.871d0,7400.0d0,9000.0d0/
      data d(1,6,11),d(2,6,11),d(3,6,11)/1.015d0,6400.0d0,9500.0d0/
      data d(1,7,11),d(2,7,11),d(3,7,11)/1.500d0,2*3.75d3/
      data d(1,8,11),d(2,8,11),d(3,8,11)/8.331d0,38.7d0,345.0d0/
      data d(1,9,11),d(2,9,11),d(3,9,11)/10.000d0,2*2.24d2/
      data d(1,10,11),d(2,10,11),d(3,10,11)/100.00d0,2*0.224d0/
      data d(1,11,11),d(2,11,11),d(3,11,11)/100.00d0,2*0.224d0/
      data d(1,12,11),d(2,12,11),d(3,12,11)/100.00d0,2*0.224d0/
c   gold
      data factor(13)/1.896d-7/
      data d(1,1,13),d(2,1,13),d(3,1,13)/0.0915d0,2*33400/
      data d(1,2,13),d(2,2,13),d(3,2,13)/0.1085d0,2*17100/
      data d(1,3,13),d(2,3,13),d(3,3,13)/0.1140d0,2*14000/
      data d(1,4,13),d(2,4,13),d(3,4,13)/0.1328d0,2*8670/
      data d(1,5,13),d(2,5,13),d(3,5,13)/0.1487d0,2*7980/
      data d(1,6,13),d(2,6,13),d(3,6,13)/0.1511d0,2*8020/
      data d(1,7,13),d(2,7,13),d(3,7,13)/0.1717d0,2*9370/
      data d(1,8,13),d(2,8,13),d(3,8,13)/0.1833d0,2*10400/
      data d(1,9,13),d(2,9,13),d(3,9,13)/0.1926d0,2*11300/
      data d(1,10,13),d(2,10,13),d(3,10,13)/0.2122d0,2*13000/
      data d(1,11,13),d(2,11,13),d(3,11,13)/0.2770d0,2*15200/
      data d(1,12,13),d(2,12,13),d(3,12,13)/0.3177d0,2*15300/
      data d(1,13,13),d(2,13,13),d(3,13,13)/0.3924d0,2*15400/
      data d(1,14,13),d(2,14,13),d(3,14,13)/0.3953d0,2*15400/
      data d(1,15,13),d(2,15,13),d(3,15,13)/0.4522d0,2*13700/
      data d(1,16,13),d(2,16,13),d(3,16,13)/0.5113d0,2*12100/
      data d(1,17,13),d(2,17,13),d(3,17,13)/0.5249d0,2*11800/
      data d(1,18,13),d(2,18,13),d(3,18,13)/0.5563d0,2*10900/
      data d(1,19,13),d(2,19,13),d(3,19,13)/0.5728d0,2*10500/
      data d(1,20,13),d(2,20,13),d(3,20,13)/0.6374d0,2*10200/
      data d(1,21,13),d(2,21,13),d(3,21,13)/0.6768d0,2*9290/
      data d(1,22,13),d(2,22,13),d(3,22,13)/0.7050d0,2*8790/
      data d(1,23,13),d(2,23,13),d(3,23,13)/0.7762d0,2*7950/
      data d(1,24,13),d(2,24,13),d(3,24,13)/0.8515d0,2*6760/
      data d(1,25,13),d(2,25,13),d(3,25,13)/0.9297d0,2*5790/
      data d(1,26,13),d(2,26,13),d(3,26,13)/1.0117d0,2*4970/
      data d(1,27,13),d(2,27,13),d(3,27,13)/1.0410d0,2*4700/
      data d(1,28,13),d(2,28,13),d(3,28,13)/1.1880d0,2*3620/
      data d(1,29,13),d(2,29,13),d(3,29,13)/1.2536d0,2*3250/
      data d(1,30,13),d(2,30,13),d(3,30,13)/1.4867d0,2*2310/
      data d(1,31,13),d(2,31,13),d(3,31,13)/1.7400d0,2*1660/
      data d(1,32,13),d(2,32,13),d(3,32,13)/2.0424d0,2*1180/
      data d(1,33,13),d(2,33,13),d(3,33,13)/2.1659d0,2*1040/
      data d(1,34,13),d(2,34,13),d(3,34,13)/2.206d0,914.34,2489.2/
      data d(1,35,13),d(2,35,13),d(3,35,13)/2.291d0,2284.3,3333.2/
      data d(1,36,13),d(2,36,13),d(3,36,13)/2.743d0,2189.5,2522.9/
      data d(1,37,13),d(2,37,13),d(3,37,13)/3.150d0,1816.5,1929.6/
      data d(1,38,13),d(2,38,13),d(3,38,13)/3.425d0,1577.9,1648.3/
      data d(1,39,13),d(2,39,13),d(3,39,13)/3.6917d0,2*1510/
      data d(1,40,13),d(2,40,13),d(3,40,13)/4.4663d0,2*953/
      data d(1,41,13),d(2,41,13),d(3,41,13)/4.5108d0,2*930/
      data d(1,42,13),d(2,42,13),d(3,42,13)/4.9522d0,2*737/
      data d(1,43,13),d(2,43,13),d(3,43,13)/5.4147d0,2*589/
      data d(1,44,13),d(2,44,13),d(3,44,13)/5.8988d0,2*474/
      data d(1,45,13),d(2,45,13),d(3,45,13)/6.9303d0,2*314/
      data d(1,46,13),d(2,46,13),d(3,46,13)/7.4782d0,2*258/
      data d(1,47,13),d(2,47,13),d(3,47,13)/8.0478d0,2*214/
      data d(1,48,13),d(2,48,13),d(3,48,13)/8.6389d0,2*178/
      data d(1,49,13),d(2,49,13),d(3,49,13)/9.8864d0,2*126/
      data d(1,50,13),d(2,50,13),d(3,50,13)/11.92d0,75.2d0,187.1d0/
      data d(1,51,13),d(2,51,13),d(3,51,13)/13.73d0,128.1d0,177.36d0/

      e = e0
      if(jz.eq.79) e = Raytrace_energy_change(e0)
      Raytrace_clcbet=0.0
      do 10 i=1,13
        if(jz.eq.z(i)) goto 20
 10   continue

      xmu=310415926d+31
      return

 20   continue
      if(jz.eq.79) goto26
      do 25 j=2,12
        if(e.lt.d(1,j,i)) goto 30
 25   continue
      j=12
      goto 30
 26   continue
      do 27 j=2,51
        if(e.lt.d(1,j,i)) goto 30
 27   continue


 30   xml=dlog(d(2,j,i)/d(3,j-1,i))
      el2=dlog(d(1,j,i)/d(1,j-1,i))
      el=dlog(e/d(1,j-1,i))

      xlogmu=dlog(d(3,j-1,i))+el*xml/el2
      xmu=exp(xlogmu)
c     Raytrace_clcbet=1.896d-7*xmu/e
      Raytrace_clcbet=factor(i)*xmu/e
ccccccccccccccccccccc  added at 1997.10.10 version 6.1_t2 (edge tune)
      if((e.ge.1.222).and.(e.lt.2.206)) then
            Raytrace_clcbet = Raytrace_clcbet*(-0.1652*e+1.202)
      endif
      return
      end

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C--------1---------2---------3---------4---------5---------6---------7--
C FUNCTION : Raytrace_energy_change
C     M IV -36eV, M V -34eV shift 97/10/7 Ishisaki & Watanabe version 6.1_t1
C     M IV -25eV, M V -25eV shift 99/02/3 Y.ISHISAKI version 6.2_t1
C
      Real*8 function Raytrace_energy_change(e0)
      Implicit None
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

c input
      Real*8 e0
c local
      Real*8 e
c begin
      if (e0.le.2.145d0 .or. e0.ge.2.743d0) then
         e = 0.d0
      else if (e0.gt.2.145d0 .and. e0.le.2.206d0) then
         e = (e0-2.145d0)*0.025d0/(2.206d0-2.145d0)
      else if (e0.gt.2.206d0 .and. e0.le.2.291d0) then
         e = 0.025d0
      else
         e = 0.025d0 - (e0-2.291d0)*0.025d0/(2.743d0-2.291d0)
      endif
c
      Raytrace_energy_change = e0 - e
      return
      end

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C--------1---------2---------3---------4---------5---------6---------7--
C SUBROUTINE : STPRNDINIT
C     step function initialization for stprnd     94.09.22 A. Furuzawa
C
      subroutine Raytrace_stprndinit()
      Implicit None
C
C     Output
C       stpx1 :R8  :left edge position of mid step
C       stpx2 :R8  :left edge position of low step
C       stpy0 :R8  :hight of top step
C       stpy1 :R8  :hight of mid step
C       stpy2 :R8  :hight of low step
C       stps0 :R8  :area of top step
C       stps1 :R8  :area of mid step
C       stps2 :R8  :area of low step
C
C       * : not available now
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      common /Raytrace_v_stprnd/ stpx1,stpx2,stpy0,stpy1,stpy2,stps0,
     &     stps1,stps2,stpc
      common /Raytrace_v_fluct/ ard,brd,crd,xmaxrd

c constant
      Real*8 M_PI
      Parameter (M_PI = 3.1415926535897932385d0)
C     For output
      real*8 stpx1,stpx2,stpy0,stpy1,stpy2,stps0,stps1,stps2,stpc
C     LOCAL
      real*8 y0,y1,y2,s0,s1,s2,x1,x2,c,ard,brd,crd,xmaxrd
C     EXTERNAL
      real*8 Raytrace_f_fluct

      data x1/5.d-1/,x2/5.d0/

      write(6,*) 'stprnd init'
      stpx1=x1/6.d1/1.8d2*M_PI
      stpx2=x2/6.d1/1.8d2*M_PI
c      write(6,*) stpx1,stpx2
      y0=1.1d0*Raytrace_f_fluct(0.d0)
      s0=y0*stpx1
c      write(6,*) y0,s0
      y1=1.1d0*Raytrace_f_fluct(stpx1)
      s1=y1*(stpx2-stpx1)
c      write(6,*) y1,s1
      y2=1.1d0*Raytrace_f_fluct(stpx2)
      s2=y2*(xmaxrd-stpx2)
c      write(6,*) y2,s2
      stpc=2.d0*(s0+s1+s2)

      stpy0=y0/stpc
      stpy1=y1/stpc
      stpy2=y2/stpc
      stps0=s0/stpc
      stps1=s1/stpc
      stps2=s2/stpc

c      write(6,*) stps0+stps1+stps2

      return
      end

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C--------1---------2---------3---------4---------5---------6---------7--
C SUBROUTINE : STPRND
C     generate random number of step-like probability function
C                                                 94.09.22 A. Furuzawa
C
      subroutine Raytrace_stprnd (rx,y)
      Implicit None
C
C     Input
C       stpx1 :R8  :left edge position of mid step
C       stpx2 :R8  :left edge position of low step
C       stpy0 :R8  :hight of top step
C       stpy1 :R8  :hight of mid step
C       stpy2 :R8  :hight of low step
C       stps0 :R8  :area of top step
C       stps1 :R8  :area of mid step
C       stps2 :R8  :area of low step
C       stpc  :R8  :factor
C     Output
C       rx    :R8  :generated random number(rad.)
C       y     :R8  :height of step function * stpc
C       * : not available now
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      common /Raytrace_v_stprnd/ stpx1,stpx2,stpy0,stpy1,stpy2,stps0,
     &     stps1,stps2,stpc

C     For INPUT
      real*8 stpx1,stpx2,stpy0,stpy1,stpy2,stps0,stps1,stps2,stpc
C     For OUTPUT
      real*8 rx,y
C     For LOCAL
      real*8 ars, rs
c function
      Real*8 Raytrace_ran
c
      rs=Raytrace_ran()-5.d-1
      ars=dabs(rs)

      if (ars.le.stps0) then
         rx=ars/stpy0
         y=stpy0*stpc
      else if (ars.gt.stps0.and.ars.le.(stps0+stps1)) then
         rx=(ars-stps0)/stpy1+stpx1
         y=stpy1*stpc
      else if (ars.gt.(stps0+stps1).and.ars.le.5.d-1) then
         rx=(ars-(stps0+stps1))/stpy2+stpx2
         y=stpy2*stpc
      else
         write(6,*) 'ERROR : RX outside of region'
         stop
      endif

      if (rs.le.0.0d0) rx = rx * -1.d0

      return
      end

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C--------1---------2---------3---------4---------5---------6---------7--
C SUBROUTINE : NVFLUCT
C     generate random number of user-defined probability function
C                                                 94.09.22 A. Furuzawa
C
      Real*8 function Raytrace_nvfluct()
      Implicit None
C
C     Input
C       gamma :R8  :lorenzian's gamma(arcmin.)
C       a     :R8  :inclination of linear function(/arcmin.)
C
C     Output
C       nvfluct:R8 :user-defined distribution random number(rad.)
C       * : not available now
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

C     For OUTPUT
c      real*8 nvfluct

C     LOCAL
      real*8 ry,fx

C     For STPRND
      real*8 rx,y

C     EXTERNAL
      real*8 Raytrace_f_fluct, Raytrace_ran

 10   call Raytrace_stprnd(rx,y)
      ry=Raytrace_ran()*y
      fx=Raytrace_f_fluct(rx)

      if (ry.gt.fx) goto 10

      Raytrace_nvfluct=rx

      return
      end

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C--------1---------2---------3---------4---------5---------6---------7--
C SUBROUTINE : F_FLUCTINIT
C     initialize user-defined probability function
C                                               94.09.22 A. Furuzawa
C
      subroutine Raytrace_f_fluctinit(gamma,a)
      Implicit None
C
C     Input
C       gamma :R8  :lorenzian's gamma(arcmin.)
C       a     :R8  :inclination of linear function(/arcmin.)
C
C       * : not available now
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      common /Raytrace_v_fluct/ ard,brd,crd,xmaxrd

c constant
      Real*8 M_PI
      Parameter (M_PI = 3.1415926535897932385d0)
C     For INPUT
      real*8 gamma,a
C     For OUTPUT
      real*8 ard,brd,crd,xmaxrd
C     For LOCAL
      real*8 b,xmax

      data xmax/1.2d2/

      write(6,*) 'f_func init'
      ard = a*6.d1*1.8d2/M_PI
      b=gamma/2.d0
      brd = b/6.d1/1.8d2*M_PI
      xmaxrd=xmax/6.d1/1.8d2*M_PI
c      write(6,*) a,b,xmax
c      write(6,*) ard,brd,crd,xmaxrd
      crd = 1.d0/(ard*dlog(xmax*xmax/b/b+1.d0)+2.d0/brd*datan(xmax/b))

      return
      end

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C--------1---------2---------3---------4---------5---------6---------7--
C FUNCTION : F_FLUCT
C     user-defined probability function
C                                               94.09.22 A. Furuzawa
C
      Real*8 function Raytrace_f_fluct(x)
      Implicit None
C
C     Input
C       x     :R8 :x(rad.)
C     Output
C       nvfluct:R8 :user-defined distribution random number(rad.)
C       * : not available now
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      common /Raytrace_v_fluct/ ard,brd,crd,xmaxrd

C     For INPUT
      real*8 x,ard,brd,crd,xmaxrd
C     For OUTPUT
c      real*8 f_fluct

      Raytrace_f_fluct=crd*(ard*dabs(x)+1.d0)/(x*x+brd*brd)

      return
      end

