;
; For HCCSSS calculation between Rcp and Rss
; H(n,..)={(n+1)/[Rcp^2(Rcp+a)^n] + n*(Rcp+a)^(n+1)/[Rcp^2(Rss+a)^(2n+1)]}^(-1)
;         {(a+r)^{-(n+1)}-(a+r)^n/(Rss+a)^(2n+1)}  n=0,1,...Nmax
; dHr(n,,..)=dH(n,..)/dr=
;    -{(n+1)/[Rcp^2(Rcp+a)^n] + n*(Rcp+a)^(n+1)/[Rcp^2(Rss+a)^(2n+1)]}^(-1)
;     {(n+1)/[r^2(r+a)^n] + n*(r+a)^(n+1)/[r^2(Rss+a)^(2n+1)}
; MODIFIED: 23APR2003 	Xuepu Zhao for nmax greater than 15
; H(n)=[Rcp^2/(r+a)][((rcp+a)/r+a))^n]*
;      [1-((r+a)/(Rss+a))^(2n+1)]/[n+1+n((Rcp+a)/(Rss+a))^(2n+1)]
; dHr(n)=-[Rcp^2/r^2][(Rcp+a)/(r+a))^n]*
;         [n+1+n((r+a)/(Rss+a))^(2n+1)]/[n+1+n((Rcp+a)/(Rss+a))^(2n+1)]
;
pro cssshf,Nmax,apar,Rcp,Rss,r,Hr,dHr,d2Hr
    Hr = fltarr(Nmax+1)
    dHr = fltarr(Nmax+1)
    d2Hr = fltarr(Nmax+1)
    arc = Rcp + apar
    ars = Rss + apar
    ar = r + apar
    arcoar=arc/ar
    aroars=ar/ars
    arcoars=arc/ars
    rc2 = Rcp * Rcp
    r2 = r * r
    for n = 0.0, Nmax do begin
        arcoarn=arcoar^n
        l=2*n+1
        aroarsl=aroars^l
        arcoarsl=arcoars^l
        Hr(n)=(1-aroarsl)/(n+1+n*arcoarsl)
        Hr(n)=Hr(n)*arcoarn*rc2/ar
        dHr(n)=-(n+1+n*aroarsl)/(n+1+n*arcoarsl)
        dHr(n)=dHr(n)*arcoarn*rc2/r2
;       arsn = ars^(2*n+1)
;       n1 = n+1
;       un = n1/(rc2*arc^n) + n*arc^n1/(rc2*arsn)
;       Hr(n) = (1.0/ar^n1 - ar^n/arsn)/un
;       dHr(n) = -(n1/(r2*ar^n) + n*ar^n1/(r2*arsn))/un
                     ; modified 10 March 1998
    endfor

end