PROGRAM absor

c      f77 -o absor absorption.f
C
C   Absorption correction for powder diffraction for annular cylinder
C****************************************************************************
C     By D.Schmitt and B. Ouladdiaf,  J.Appl.Crys.31(1998)620
C****************************************************************************

c
       IMPLICIT REAL*8 (A-H,O-Z)
       CHARACTER filout*16,rep*1
C
       common /rayon/ RP, precis, PI
c
     9 format (8x,'teta ',3x,'sin2(teta)',4x,'mu*R ',6x,
      1 'transm (Rp=',F7.4,')')
    13 format (1x,3(2x,f10.4),3x,e13.6)
c
       precis = 1.d-14
       pi = 4.D0*datan(1.D0)

       TYPE *,'   '
       write(*,*) '           Program absor to calculate Absorption '
       write(*,*) '       Ceofficients for annular cylinder By Simpson
      &method '

       TYPE *,'   '
C
  200      write (*,*)' Give the output filename ?'
       read(*,'(a)') filout
       open(32,file=filout,status='unknown')

        write (*,*)' Give sin^2(theta_ini), sin^2(theta_fin) and Step in
      &sin^2(theta): ?'


        read *, stetai,stetaf,sdteta
        write(*,*),'Give the Linear Absorption Coefficient muR(muRext):?'

        read *,xmui


    15 print *, 'Give the Radius ratio = Rint/Rext ( <1 ) ? '

       read *, RP

       if (RP.GT.1. .OR. RP.LT.0.) go to 15


        write (*,*)'Please wait ...'
       nr = 100
       nf = 100


c
       dr = (1.D0-RP)/dfloat(nr)
       df = 2.D0*pi/dfloat(nf)
       if (sdteta.EQ.0.) sdteta = 1.
       if (dmu.EQ.0.) dmu = 1.

       nteta = dabs(stetaf - stetai)/dabs(sdteta) + 1

       write(32,9) RP

       do 99 iteta = 1,nteta
       steta = stetai + (iteta - 1) * sdteta
       teta = dasin(dsqrt(steta))
       tetaD = teta * 180. / pi
c

        xmu = xmui
c
       xxl = 0.
c
       r = RP
       xl2 = 0.
       fi=0.
       call chemin (r,fi,teta,y2)
c      call chemin (r,0.,teta,y2)
c
       do 30 ifi = 1, nf
       y0 = y2
       fi = ifi * df
       call chemin (r,fi,teta,y2)
       fi = fi - 0.5 * df
       call chemin (r,fi,teta,y1)
       xl2 = xl2 + r * (dexp(-xmu*y0) + 4.*dexp(-xmu*y1) +
      1 dexp(-xmu*y2)) / 6.
    30 continue
c
       do 50 ir = 1,nr
       xl0 = xl2
c
       r = RP + ir * dr
       xl2 = 0.
       fi=0.
       call chemin (r,fi,teta,y2)
c      call chemin (r,0.,teta,y2)

c
       do 40 ifi = 1, nf
       y0 = y2
       fi = ifi * df
       call chemin (r,fi,teta,y2)
       fi = fi - 0.5 * df
       call chemin (r,fi,teta,y1)
       xl2 = xl2 + r * (dexp(-xmu*y0) + 4.*dexp(-xmu*y1) +
      1 dexp(-xmu*y2)) / 6.
    40 continue
c
       r = r - 0.5 * dr
       xl1 = 0.
       fi=0.
       call chemin (r,fi,teta,y2)
c      call chemin (r,0.,teta,y2)
c
       do 45 ifi = 1, nf
       y0 = y2
       fi = ifi * df
       call chemin (r,fi,teta,y2)
       fi = fi - 0.5 * df
       call chemin (r,fi,teta,y1)
       xl1 = xl1 + r * (dexp(-xmu*y0) + 4.*dexp(-xmu*y1) +
      1 dexp(-xmu*y2)) / 6.
    45 continue
c
       xxl = xxl + (xl0 + 4.*xl1 + xl2) / 6.
c
    50 continue
c
       xxl = xxl * dr * df / (pi * (1. - RP*RP))
c
c   98 print 10, tetaD, steta, xmu, xxl

       write(32,13)tetaD,Steta,xmu,xxl

    99 continue
	close(32)
       write(*,*) ' Do you want to continue (Y/N) ?  <CR>=N '
       read(*,14) rep
c      write(*,*),'rep=',rep
    14 format(a1)
       If(rep.EQ.' ')rep = 'N'
       if(rep.EQ.'Y'.OR.rep.EQ.'y')go to 200
        write(*,*),' Good Bye.  B.O.'

c

       stop
       end
c
       subroutine chemin (r, fi, teta, xl)
c

       IMPLICIT REAL*8 (A-H,O-Z)
C
       common /rayon/ RP, precis, PI
c
       xm = r * dcos(fi)
       ym = r * dsin(fi)
c
       xe = -dsqrt(1.D0 - ym*ym)
c
       if (teta*180./pi.EQ.45.) go to 8
       t2t = dtan(2.*teta)
       if (dabs(t2t).GT.(1./precis)) go to 8
c
       dsig = dsign(1.D0,t2t)
       par = dsin(fi) - dcos(fi)*t2t
       disc = 1.D0 + t2t*t2t - r*r*par*par
       xs = (-r*t2t*par + dsig * dsqrt(disc)) / (1.D0 + t2t*t2t)
       ys = xs * t2t + r*par
       go to 9
     8 xs = xm
       ys = dsqrt(1.D0 - xs*xs)
c
     9 xem = xm - xe
       xms = (xs-xm)*(xs-xm) + (ys-ym)*(ys-ym)
       xms = dsqrt(xms)
c
       xl = xem + xms
c
       if (RP.eq.0.) return

       ps = dabs(-xm * dsin(2.*teta) + ym * dcos(2.*teta))
       if (ps.GE.RP) go to 10
       ps2 = xm*xs + ym*ys - r*r
       if (ps2.GE.-precis) go to 10
c

       xv = 2.*dsqrt(rp*rp - ps*ps)
       xl = xl - xv
c

c
    10 if (dabs(ym).GE.RP) go to 20
       if (xm.LE.0.) go to 20
c

       xv2 = 2.*dsqrt(rp*rp - ym*ym)
       xl = xl - xv2
c
    20 if (dabs(xl).le.precis) xl=0.
       if (xl.LT.0.) then
         print *, 'erreur de chemin ', xl,r,fi
        endif
       return
       end