;
;  pdpj:        Calculate Schmidt-type associated Legendre function using
;               P_m^m=(2-dm0)^0.5 ((2m-1)!!/2m!!)^0.5 (1-x^2)^{m/2}
;                     dm0=1 when m=0 and dm0=0 when m>0
;  P_{l+1}^m=((l+1)^2-m^2)^{-0.5}[(2l+1) x P_l^m - (l^2-m^2)^0.5 P_{l-1}^m]
;               dP_m^m=m (2-dm0)^0.5 ((2m-1)!!/2m!!)^0.5 x (1-x^2)^{(m-1)/2}
;  dP_{l+1}^m=((l+1)^2-m^2)^{-0.5} [(2l+1)[-(1-x^2)^0.5 P_l^m + x dP_l^m] -
;                                    (l^2-m^2)^0.5 dP_{l-1}^m]
;  input:       cthj, Nmax
;  output:      Pj(Nmax+1,Nmax+1), dPj(Nmax+1,Nmax+1) for specified thetad
;               including (Nmax+1)(Nmax+2)/2 elements
;  wrottten:    Xuepu Zhao      Feb. 28, 1996
;
pro pdpj06,cthj,Nmax,pj,dpj

    sthj=(1 - cthj*cthj)^0.5
    nm1=Nmax+1
    pj=fltarr(nm1,nm1)
    pj=double(pj)
    dpj=pj
; get pj, dpj
    for m=0.0,Nmax do begin
        if m eq 0 then t1=1.0 else t1=2.0^0.5
        dt1=t1*m*cthj/sthj
        for mi=0, m-1 do begin
            mm=((mi+mi+1.0)/(mi+mi+2.0))^0.5
            t1=t1*sthj*mm  ; get P_m^m
            t1=double(t1)
            dt1=dt1*sthj*mm ; get dP_m^m
            dt1=double(dt1)
        endfor
        pj(m,m)=t1
        dpj(m,m)=dt1
        if m lt Nmax then begin
           m2=m*m
           t0=0
           dt0=0
           s1=0
           for n=m, Nmax-1 do begin
               l2=n+1.0
               l22=l2*l2
               s2=(l22-m2)^0.5
               f1=2*n+1.0
               t2=(f1*cthj*t1 - s1*t0)/s2  ; get P_m+1^m,...,P_n^m
               dt2=(f1*(cthj*dt1 - sthj*t1) - s1*dt0)/s2  
                                           ;get dP_m+1^m,.,dP_n^m
               pj(n+1,m)=t2
               dpj(n+1,m)=dt2
               s1=s2
               t0=t1
               t1=t2
               dt0=dt1
               dt1=dt2
           endfor
        endif
    endfor

end