;
; PURPOSE: Calculate field line starting at Rss, phri, thrj using HCCSSS model
; NOTE: g,h and gc,hc should be in same unit, Gauss or mT
; WRITTEN: 20AUG2003	Xuepu Zhao 
; MODIFIED: 10FEB2006	Xuepu Zhao  Add 'fte' for Nagoya group
; MODIFIED: 15Mar2006   Xuepu Zhao  Add keyword r1
;
pro iflij_csss,g,h,gc,hc,phri,thrj,fte,np,xw,yw,zw,rw,tw,pw,$
               R1=r1,APAR=apar,RCP=rcp,RSS=rss,RBOTTOM=rbottom

    if N_params( ) lt 1 then begin
       print,'pro iflij_csss,g,h,gc,hc,phri,thrj,brss,brcp,br1,fteij,$'
       print,'               phd1,thd1,apar=,rcp=,rss=,/ss,rbottom='
       print,'default apar=0.2,rcp=2.25,rss=15.0,rbottom=1.0'
       return
    endif
;
    if not keyword_set(apar) then apar=0.2
    if not keyword_set(rcp) then rcp=2.25
    if not keyword_set(rss) then rss=15.0
    if not keyword_set(rbottom) then rbottom=1.0
;
    d2r=!DTOR
    nao=!VALUES.F_NAN
    dstep=0.01            ; to be modified to save computer time
    nstep=1000
    xv=fltarr(nstep)
    yv=xv
    zv=xv
    rv=xv & tv=xv & pv=xv
; 
    cph=cos(phri)
    sph=sin(phri)
    cth=cos(thrj)
    sth=sin(thrj)
    if not keyword_set(r1) then r1=rss
    x = r1*sth*cph
    y = r1*sth*sph
    z = r1*cth
    xv(0)=x & yv(0)=y & zv(0)=z
    rv(0)=r1 & tv(0)=thrj/!DTOR & pv(0)=phri/!DTOR
    if r1 gt rcp then begin
       pb0_csss,gc,hc,r1,phri,cth,sth,apar,rcp,rss,Br,Bt,Bp
       sig = -1.
       Br1=Br
       r=r1
    endif else begin
       ppb_hc,g,h,r1,phri,cth,sth,Br,Bt,Bp
       Br1=ABS(Br)
       is1=1
       r=r1
       goto,lbl1
    endelse
;  calculate the positive vector field at a point between Rcp and Rss
    for is1=1,nstep-1 do begin
        BB = sqrt(Br*Br+Bt*Bt+Bp*Bp)
        r0=r
        ds = sig * r * dstep / BB
        dx = ds*(Br*sth*cph+Bt*cth*cph-Bp*sph)
        dy = ds*(Br*sth*sph+Bt*cth*sph+Bp*cph)
        dz = ds*(Br*cth-Bt*sth)
        x = x + dx
        y = y + dy
        z = z + dz
        r = sqrt(x*x + y*y + z*z)
        if r ge r0 then begin 		; Occur often near cusp points
           ftpij=nao
           phd1=nao
           thd1=nao
           return
        endif
        if abs(z) eq 0.0 then begin
           thr=!pi/2.0
        endif else begin
           arg=sqrt(x*x +y*y)/z
           thr=atan(arg)
           if z lt 0.0 then thr=thr+!pi
        endelse
        phr=atan(y,x)
        if phr lt 0.0 then phr=phr+2*!pi
        xv(is1)=x
        yv(is1)=y
        zv(is1)=z
        rv(is1)=r
        tv(is1)=thr/d2r
        pv(is1)=phr/d2r
        cth=COS(thr)
        sth=SIN(thr)
        cph=COS(phr)
        sph=SIN(phr)
        if r le rcp then goto,lbl1
        pb0_csss,gc,hc,r,phr,cth,sth,apar,Rcp,Rss,Br,Bt,Bp
    endfor
 lbl1: print,'is1: ',is1 
;
; calculate field line between Rs and Rcp
;
    if r1 gt rcp then ppb_hc,g,h,r,phr,cth,sth,Br,Bt,Bp
    if Br gt 0 then sig = -1 else sig = 1
    for is2=is1,nstep-1 do begin
        r0=r
        BB = sqrt(Br*Br + Bt*Bt+ Bp*Bp)
        ds = sig * r * dstep / BB
        dx = ds*(Br*sth*cph+Bt*cth*cph-Bp*sph)
        dy = ds*(Br*sth*sph+Bt*cth*sph+Bp*cph)
        dz = ds*(Br*cth-Bt*sth)
        x = x + dx
        y = y + dy
        z = z + dz
        r = sqrt(x*x + y*y + z*z)
        if r ge r0 then begin
           br1=nao
           fteij=nao
           phd1=nao
           thd1=nao
           return
        endif
        if r le rbottom then begin
           np=is2
           fte=(Br/Br1)/rss^2   	; Br here is not exactly at r=1 
           goto, lbl2
        endif
        if z eq 0.0 then begin
           thr=!pi/2.0
        endif else begin
           arg=sqrt(x*x +y*y)/z
           thr=atan(arg)
           if z lt 0.0 then thr=thr+!pi
        endelse
        phr = atan(y,x)
        if phr lt 0.0 then phr=phr+2*!pi
        xv(is2)=x
        yv(is2)=y
        zv(is2)=z
        rv(is2)=r
        tv(is2)=thr/d2r
        pv(is2)=phr/d2r
        cth=COS(thr)
        sth=SIN(thr)
        cph=COS(phr)
        sph=SIN(phr)
        ppb_hc,g,h,r,phr,cth,sth,Br,Bt,Bp
    endfor
  lbl2:
    xw=xv(0:np-1)
    yw=yv(0:np-1)
    zw=zv(0:np-1)
    rw=rv(0:np-1)
    tw=tv(0:np-1)
    pw=pv(0:np-1)
    np=-1*sig*np

end