subroutine intlog5a(a,b,g,d,e,x1,x2,ans1)

c - Returns the analytical solution to the integral:
c     Int_wrt_x_from_x1_to_x2 of Log[a*x + b + Sqrt[g*x**2 + 2*d*x + e]]

c - For the most part, the answer is "BB# - AA#" for #=1 through 7
c - except for #6 and #7 when the BB-AA term is combined into
c - a single term.

c - This subroutine is set up to explicitly exploit the
c - limited values of a,b,g,d,e which will be used
c - as part of the gravitational attraction integral.
c - A more generic form of this subroutine, called intlog5,
c - works for all a,b,g,d,e values

      implicit none
      integer*4 l 
      real*16 a,b,g,d,e,x1,x2

c - Included here are variables that are PROVED to never be complex
      real*16 insides1,s1
      real*16 a2,b2,g2,d2,e2,a3,a4,ab,abd,bdg,eg
      real*16 h01,h02,h03,h04,h05,h06,h10,h11,h12,h13,h14,h15,h16
      real*16 inatanA,inatanB
      real*16 rinlogA,rinlogB
      real*16 AA1,BB1,AA2,BB2

      complex*32 ans1


c - Included here are variables that are PROVED to at least
c - POSSIBLY be complex
      complex*32 i
      complex*32 qcsqrt,ccsqrt
      complex*32 s2,s3,s4A,s4B
      complex*32 h07,h08,h09
      complex*32 AA3,BB3,AA4,BB4
      complex*32 AA5,BB5,AA6,BB6,AA7,BB7
      complex*32 cinlogA,cinlogB
      complex*32 cmplxa,loga


      complex*32 Q01,Q02,Q03,Q04,Q05
      complex*32 C0Q06,C1Q06,Q06A,Q06B
      complex*32 C0Q07,C1Q07,Q07A,Q07B
      complex*32 Q08A,Q08B,Q09A,Q09B,Q10A,Q10B,Q11A,Q11B
      complex*32 Q12A,Q12B,Q13A,Q13B,Q14A,Q14B,Q15A,Q15B
      complex*32 C0Q16,C1Q16,Q16A,Q16B
      complex*32 C0Q17,C1Q17,Q17A,Q17B
      complex*32 Q18A,Q18B
      complex*32 p1Q02,p2Q02,p3Q02,p4Q02,p5Q02

      complex*32 BB6mAA6,BB7mAA7,BB3mAA3,BB4mAA4,BB5mAA5


      complex*32 clog1A,clog2A,inlog1A,inlog2A
      complex*32 clog1B,clog2B,inlog1B,inlog2B

      common /root/ l

c - Get some frequently used constants:

       i = (0.q0,1.q0)
      a2 = a*a
      b2 = b*b
      g2 = g*g
      d2 = d*d
      e2 = e*e

      a3 = a2*a
      a4 = a2*a2

      ab = a*b
      abd = ab*d
      bdg = b*d*g
      eg = e*g

c - Force s1 to be non-negative (it's impossible, but
c - round off errors might make it very slightly
c - negative near zero)
      insides1 = 2*abd - d2 - a2*e - b2*g + eg
      write(6,*) ' inside s1:',insides1
      if(abs(insides1).lt.1q-14)insides1=0.q0

      s1 = sqrt(insides1)
      s2 = qcsqrt(g)
c     s3 = qcsqrt(-2*abd + d2 + a2*e + b2*g - eg)
      s3 = i*s1

      s4A = qcsqrt(e + 2*d*x1 + g*x1**2)
      s4B = qcsqrt(e + 2*d*x2 + g*x2**2)

      h01 = ab - d
      h02 = a2 - g
      h03 = s1/h02
      h04 = b2 - e
      h05 = h01/2.q0/h02
      h06 = a*d - b*g
      h07 = h06/h02/s2
      h08 = 2*d/s2
      h09 = 2*s2
      h10 = -2*abd + a2*e + b2*g
      h11 = b*d - a*e
      h12 = -2*abd + d2 + a2*e + (b2-e)*g
      h13 = (a2 - g)**2
      h14 = -abd + d2 + a2*e - eg
      h15 = 3*a2*b*d - a3*e + bdg + a*(-2*d2 + (-2*b2 +e)*g)
      h16 = a2*b - 2*a*d + b*g

c - Get the easy ones first

cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c - Do the C1 term:
cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      AA1 = -x1
      BB1 = -x2    

      write(6,*) ' At C1:'
      write(6,*) '    C1(x2) = ',BB1
      write(6,*) '    C1(x1) = ',AA1
      write(6,*) '    diff   = ',BB1-AA1



cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c - Do the C2 term:
cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c     A2 = h03*atan((h01 + h02*x1)/s1)
c     B2 = h03*atan((h01 + h02*x2)/s1)
c     A2 = h03*atan(inatanA)
c     B2 = h03*atan(inatanB)
c     AA2 = h03*(i/2)*log((i+inatanA)/(i-inatanA))
c     BB2 = h03*(i/2)*log((i+inatanB)/(i-inatanB))


      if(s1.eq.0)then
        write(6,*) ' Precoded C2 term for zero a3 or A3'
        AA2 = 0.q0
        BB2 = 0.q0
      elseif(b.eq.0.and.d.eq.0.and.e.eq.0)then
        write(6,*) ' Precoded C2 term for b,d,e=0'
        AA2=0.q0
        BB2=0.q0
      else
        write(6,*) ' Normal C2 computation'
        inatanA = (h01 + h02*x1)/s1
        inatanB = (h01 + h02*x2)/s1
c       write(6,*) ' inatanA, AA2 = ',inatanA
c       write(6,*) ' inatanB, BB2 = ',inatanB
        AA2 = h03*atan(inatanA)
        BB2 = h03*atan(inatanB)
      endif

      write(6,*) ' At C2:'
      write(6,*) '    C2(x2) = ',BB2
      write(6,*) '    C2(x1) = ',AA2
      write(6,*) '    diff   = ',BB2-AA2


cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c - Do the C3 term:
cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c     AA3 = h05*log(h04 + 2*h01*x1 + h02*x1**2)
c     BB3 = h05*log(h04 + 2*h01*x2 + h02*x2**2)

c - Be prepared if a3 or A3 is zero
c - In either of those cases, s1=0
      if(s1.eq.0)then
        write(6,*) ' Precoded C3 term for a3 or A3 = 0'
        BB3mAA3 = 0.q0
      elseif(a.eq.0 .and. d.eq.0)then
        write(6,*) ' Precoded C3 term for a,d = 0'
        BB3mAA3 = 0.q0 
      elseif(b.eq.0 .and. d.eq.0)then
        write(6,*) ' Precoded C3 term for b,d=0'
        BB3mAA3 = 0.q0
      else
        write(6,*) ' Normal C3 computation'
        rinlogA = h04 + 2*h01*x1 + h02*x1**2
        rinlogB = h04 + 2*h01*x2 + h02*x2**2
c       write(6,*) ' rinlogA, AA3 = ',rinlogA
c       write(6,*) ' rinlogB, BB3 = ',rinlogB
        BB3mAA3 = h05*log(rinlogB/rinlogA)
      endif

      write(6,*) ' At C3:'
      write(6,*) '    diff   = ',BB3mAA3


cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c - Do the C4 term:
cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      if(x1.eq.0)then
        write(6,*) ' Precoded C4 term for x1=0'
        AA4 = 0.q0
        BB4 = x2 * log(a*x2 + b + s4B)
        BB4mAA4 = BB4 - AA4
      elseif(x2.eq.0)then
        write(6,*) ' Precoded C4 term for x2=0'
        AA4 = x1*log(b + a*x1 + s4A)
        BB4 = 0
        BB4mAA4 = BB4 - AA4
      else
        write(6,*) ' regular C4 computation'
        AA4 = x1*log(b + a*x1 + s4A)
        BB4 = x2*log(b + a*x2 + s4B)
        BB4mAA4 = BB4 - AA4
      endif

      write(6,*) ' At C4:'
      write(6,*) '    diff   = ',BB4mAA4

cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c - Do the C5 term:
cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c     AA5 = h07*log(h08 + h09*x1 + 2*s4A)
c     BB5 = h07*log(h08 + h09*x2 + 2*s4B)

      if(a.eq.0 .and. b.eq.0)then
        write(6,*) ' Precoded C5 term for a,b=0'
        BB5mAA5 = 0.q0 
      elseif(b.eq.0 .and. d.eq.0)then
        write(6,*) ' Precoded C5 term for b,d=0'
        BB5mAA5=0.q0
      else
        write(6,*) ' Normal C5 computation'
        cinlogA = h08+h09*x1+2*s4A
        cinlogB = h08+h09*x2+2*s4B
        write(6,*)' cinlogA, AA5 = ',cinlogA
        write(6,*)' cinlogB, BB5 = ',cinlogB
        BB5mAA5 = h07*log(cinlogB/cinlogA)
      endif

      write(6,*) ' At C5:'
      write(6,*) '    diff   = ',BB5mAA5


c - Now do the tough last 2 components (6 and 7)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c - Do the C6 and C7 terms:
cccccccccccccccccccccccccccccccccccccccccccccccccccccccc

c - Be prepared for if a3=0 or A3=0. In either
c - case, s1=0.
      if(s1.eq.0)then
        write(6,*) ' Precoded C6 term for a3 or A3 = 0'
        write(6,*) ' Precoded C7 term for a3 or A3 = 0'
        BB6mAA6=0
        BB7mAA7=0
      elseif(a.eq.0 .and. b.eq.0)then
        write(6,*) ' Precoded C6 term for a,b = 0'
        write(6,*) ' Precoded C7 term for a,b = 0'
        BB6mAA6 = 0.q0 
        BB7mAA7 = 0.q0 
      elseif(b.eq.0.and.d.eq.0.and.e.eq.0)then
        write(6,*) ' Precoded C6 term for b,d,e=0'
        write(6,*) ' Precoded C7 term for b,d,e=0'
        BB6mAA6=0.q0
        BB7mAA7=0.q0
      else
        write(6,*) ' Regular C6 computation'
        write(6,*) ' Regular C7 computation'
c - Start with the coefficients of the logarithms (clog terms)
        Q01 = -2*h02*h10
  
        p1Q02 = -2*a3*b*d
        p2Q02 = a4*e
        p3Q02 = -2*abd*g
        p4Q02 = b2*g2
        p5Q02 = a2*(2*d2+(b2-e)*g)
        Q02 = -2*a3*b*d + a4*e - 2*abd*g + b2*g2 + a2*(2*d2 + (b2-e)*g)

        Q03 = 2*a*h06*s3

        Q04 = h11

        Q05 = b*s3

c - No dependence on x in the clog terms
        clog1A=(-1/Q01)*ccsqrt(Q02-Q03)*(Q04-Q05)
        clog1B=clog1A

        clog2A=(-1/Q01)*ccsqrt(Q02+Q03)*(Q04+Q05)
        clog2B=clog2A

c - Now, move on to the guts of the logs (inlog terms)
        C0Q06 = 2*d*h13*h12
        C1Q06 = 2*g*h13*h12
        Q06A = C0Q06 + C1Q06*x1
        Q06B = C0Q06 + C1Q06*x2

        C0Q07 = 2*h13*h14*s3
        C1Q07 = 2*a*h13*h06*s3
        Q07A = C0Q07 + C1Q07*x1
        Q07B = C0Q07 + C1Q07*x2

c - No x dependence on the Q08 terms, and Q08=Q02! :)
        Q08A = Q02
        Q08B = Q08A

c - No x dependence on the Q09 terms, and Q09=Q03! :)
        Q09A = Q03
        Q09B = Q09A

c - No x dependence on the Q10 terms
        Q10A =  (a2*b - 2*a*d + b*g)*h12
        Q10B =  Q10A

c - No x dependence on the Q11 terms
        Q11A = s3*h15
        Q11B = Q11A


        Q12A = 4*h13*s3*s4A
        Q12B = 4*h13*s3*s4B

c - no x dependence on the Q13 terms
        Q13A = -2*h16*h12
        Q13B = Q13A

c - no x dependence on the Q14 terms
        Q14A = 2*s3*(-h15)
        Q14B = Q14A

c - no x depdndence on the Q15 terms
        Q15A = s3
        Q15B = Q15A

        C0Q16 = -h01
        C1Q16 = -h02
        Q16A = C0Q16 + C1Q16*x1
        Q16B = C0Q16 + C1Q16*x2

        C0Q17 = h01
        C1Q17 = h02
        Q17A = C0Q17 + C1Q17*x1
        Q17B = C0Q17 + C1Q17*x2

c - no x dependence on the Q18 terms
        Q18A = -s3
        Q18B = Q18A

c - Do the inside of the logarithms for C6
        inlog1A = (Q06A-Q07A)/ccsqrt(Q08A-Q09A)/(Q10A-Q11A)/(Q17A-Q18A)
     *          + Q12A/(Q13A-Q14A)/(Q15A-Q16A)
        inlog1B = (Q06B-Q07B)/ccsqrt(Q08B-Q09B)/(Q10B-Q11B)/(Q17B-Q18B)
     *          + Q12B/(Q13B-Q14B)/(Q15B-Q16B)

c - Do the inside of the logarithms for C7
        inlog2A = (Q06A+Q07A)/ccsqrt(Q08A+Q09A)/(Q10A+Q11A)/(Q17A+Q18A)
     *          + Q12A/(Q13A+Q14A)/(Q15A+Q16A)
        inlog2B = (Q06B+Q07B)/ccsqrt(Q08B+Q09B)/(Q10B+Q11B)/(Q17B+Q18B)
     *          + Q12B/(Q13B+Q14B)/(Q15B+Q16B)

c       AA6 = clog1A * log(inlog1A)
c       BB6 = clog1B * log(inlog1B)

c       AA7 = clog2A * log(inlog2A)
c       BB7 = clog2B * log(inlog2B)


        BB6mAA6 = clog1A * log(inlog1B/inlog1A)
        BB7mAA7 = clog2A * log(inlog2B/inlog2A)

      endif

      write(6,*) ' At C6:'
      write(6,*) '    diff   = ',BB6mAA6
      write(6,*) ' At C7:'
      write(6,*) '    diff   = ',BB7mAA7

c     ans1 = (BB1-AA1) + (BB2-AA2) + (BB3-AA3) + (BB4-AA4) 
c    *    + (BB5-AA5) + (BB6-AA6) + (BB7-AA7)
      ans1 = (BB1-AA1) + (BB2-AA2) + BB3mAA3 + BB4mAA4
     *    + BB5mAA5 + BB6mAA6 + BB7mAA7
      write(6,*) ' dif        :',ans1
      return
      end