pro runsample

;1. Assign phd,thd,lad,cth,sth based on, say, iph=72, jth=30
    pzgrid,72,30,phd,thd,lad,cth,sth
;2. Calcualte ga & ha using nmax=9 and 'mdigh_1968.asc'
  ; pgh_mbrsc,'MDI.1968.fits',9,ga,ha,ghmf='mdigh_1968.asc',/ifl
    pgh_mbrsc,'MDI.1951.fits',9,ga,ha,ghmf='mdigh_1951.asc',/ifl
    print,'ga(0:2,0:2):',ga(0:2,0:2)
;3. Calculate Br,Bt,Bp at a point located at, say, r=2.0,lad(10),phd(20)
    ppb_hc,ga,ha,2.0,phd(20)*!DTOR,cth(10),sth(10),Br,Bt,Bp
    print,'Br,Bt,Bp at 2.0, ..:',Br,Bt,Bp
;4. Calculate Br,Bt,Bp on a spherical surface of r <= rcp, say r=2.0
    csmph,phd,9,cmph,smph
    pdp,cth,9,P,dP
    psb_hc,ga,ha,2.0,sth,cmph,smph,P,dP,bbr,bbt,bbp
;5. Calculate gca,hca and produce 'ghfcsss_1968.asc'  based on ga,ha,nmax=9
  ; pghfcsss_gaha,ga,ha,gca,hca,'ghfcsss_1968.asc',nmax=9
    pghfcsss_gaha,ga,ha,gca,hca,'ghfcsss_1951.asc',nmax=9
    print,'gca(0:2,0:2): ',gca(0:2,0:2)
;6. Calculate br,bt,bp at r=5,thd(10),phd(20) using ga,ha,gca,hca
    phri=phd(20)*!DTOR & thrj=thd(10)*!DTOR
    ppb_csss,ga,ha,gca,hca,5.,phri,thrj,br,bt,bp
    print,'br,bt,bp at r=5: ',br,bt,bp
;7. Calculate bbr,bbt,bbp on a spherical surface of, say, r=5
    psb_csss,ga,ha,gca,hca,5,phd,thd,bbr,bbt,bbp
;8. Calculate FTE & inward field line from 15 Rs at, say, phd(20),thd(10)
  ; iflij_csss,ga,ha,gca,hca,phri,thrj,fte,np,xw,yw,zw,rw,tw,pw
  ; print,'fte,np: ',fte,np
    phri=phd(21)*!DTOR & thrj=thd(3)*!DTOR
    iflij_csss,ga,ha,gca,hca,phri,thrj,fte,np,xw,yw,zw,rw,tw,pw
    print,'fte,np: ',fte,np
    phri=phd(30)*!DTOR & thrj=thd(14)*!DTOR
    iflij_csss,ga,ha,gca,hca,phri,thrj,fte,np,xw,yw,zw,rw,tw,pw
    print,'fte,np from r=15.0: ',fte,np
    iflij_csss,ga,ha,gca,hca,phri,thrj,fte,np,xw,yw,zw,rw,tw,pw,r1=6.0
    print,'fte,np from r=6.0: ',fte,np
    iflij_csss,ga,ha,gca,hca,phri,thrj,fte,np,xw,yw,zw,rw,tw,pw,r1=2.5
    print,'fte,np from r=2.5: ',fte,np

end