pro DERIV, rtp, drtp_ds, b, bv,mfl_model=mfl_model, nocrt=nocrt
nn = n_elements(rtp)/3
MAGNETIC_FIELD,rtp(0,*),rtp(1,*),rtp(2,*),bv,b,mfl_model=mfl_model,/deg,nocrt=nocrt
drtp_ds=[bv(0,*)/b, $
bv(1,*)/rtp(0,*)*!radeg/b, $
bv(2,*)/rtp(0,*)*!radeg/sin(rtp(1,*)/!radeg)/b]
return
end
PRO build_iso_fc_grid,ff,iso_fc,iso_fc_foot,save=save,verbose=verbose,$
nocrt=nocrt,harmonic=harmonic
if keyword_set(verbose) then verbose=1b else verbose=0b
mfl_model = 'SPV'
fsb=2.79924835996d3
eps=0.098
h0=0.02
fmin=3.
if keyword_set(harmonic) then fmin /=2.
nn = 1801
lat_pole = 0.2
theta_foot = findgen(nn)/float(nn-1)*180.
ww = where(theta_foot ge lat_pole and theta_foot le 180-lat_pole,nn)
theta_foot = theta_foot(ww)
index = lindgen(nn)
radius_foot = 1.d0/sqrt(1.d0+eps*(2.d0-eps)/(1.d0-eps)^2*(cos(theta_foot/!radeg))^2)
nsteps = 500
if keyword_set(harmonic) then nsteps *=2
rtp = dblarr(3,nn,nsteps)
fc = dblarr(nn,nsteps)
rtp(0,*,0) = radius_foot
rtp(1,*,0) = theta_foot
rtp(2,*,0) = 0
rtp_i = rtp(*,*,0)
DERIV,rtp_i,drtp_i,btot_i,bvec_i,mfl_model=mfl_model,nocrt=nocrt
h = rebin(h0*(2*(drtp_i(0,*,0) le 0)-1),3,nn)
fc (*,0) = btot_i*fsb
arret = btot_i*fsb ge fmin
index_i = index
ww_continue = where(arret, cnt_continue)
i=0
if verbose then print,i,cnt_continue,min(fc(*,i))
if cnt_continue ne 0 then begin
sgn = 2*(drtp_i(0,*) ge 0)-1
drtp = reform(drtp_i,3,nn)
btot = reform(btot_i, nn)
bvec = reform(bvec_i,3,nn)
h2=-h/2.d0 & h6=-h/6.d0
if verbose then plot,rtp(0,*),rtp(1,*)/!radeg,/iso,/polar,psym=3,yr=[0,10],xr=[-10,10]
repeat begin
index_i = index_i(ww_continue)
sgn_prev = sgn(index_i)
drtp_i = drtp(*,index_i)
h_i = h(*,index_i)
h2_i = h2(*,index_i)
h6_i = h6(*,index_i)
rtp_i = rtp(*,index_i,i)+h2_i*drtp_i
DERIV,rtp_i,dyt,bi,mfl_model=mfl_model,nocrt=nocrt
rtp_i=rtp(*,index_i,i)+h2_i*dyt
DERIV,rtp_i,dym,bi,mfl_model=mfl_model,nocrt=nocrt
rtp_i=rtp(*,index_i,i)+h_i*dym
dym=dym+dyt
DERIV,rtp_i,dyt,bi,mfl_model=mfl_model,nocrt=nocrt
rtp_i=rtp(*,index_i,i)+h6_i*(drtp_i+dyt+2.*dym)
DERIV,rtp_i,drtp_i,btot_i,bvec_i,mfl_model=mfl_model,nocrt=nocrt
sgn = 2*(drtp_i(0,*) ge 0)-1
arret_i = (btot_i*fsb ge fmin AND rtp_i(0,*) ge 1.d0/sqrt(1.d0+eps*(2.d0-eps)/(1.d0-eps)^2*(cos(rtp_i(1,*)/!radeg))^2) AND fc(index_i,i) gt btot_i*fsb)
ww_continue = where(arret_i, cnt_continue)
i++
if (cnt_continue ne 0 and i lt nsteps) then begin
rtp(*,index_i(ww_continue),i) = rtp_i(*,ww_continue)
btot(index_i(ww_continue)) = btot_i(ww_continue)
fc(index_i(ww_continue),i) = btot_i(ww_continue)*fsb
bvec(*,index_i(ww_continue)) = bvec_i(*,ww_continue)
drtp(*,index_i(ww_continue)) = drtp_i(*,ww_continue)
if verbose then oplot,rtp(0,*,i),rtp(1,*,i)/!radeg,/polar,psym=3
if verbose then print,i,cnt_continue,min(fc((where(fc(*,i) gt 0)),i))
endif
endrep until cnt_continue eq 0 or i eq nsteps-1
endif else begin
rtp = reform(rtp_i, 3,nn)
drtp = reform(drtp_i,3,nn)
btot = reform(btot_i, nn)
bvec = reform(bvec_i,3,nn)
endelse
fi_abc_08 = 800l+lindgen(8)
fi_abc_16 = 1600l+lindgen(16)
fi_abc_32 = 3200l+lindgen(32)
fi_a = [fi_abc_08,fi_abc_16,fi_abc_32]
fi_b = [fi_abc_08,fi_abc_16,fi_abc_32]+10000000l
fi_c = [fi_abc_08,fi_abc_16,fi_abc_32]+20000000l
fi_hf_1 = 100l + (1L+lindgen(641))*10000l
fi_hf_2 = rebin(reform(200l + lindgen(2),2,1),2,641)+rebin(reform((1L+lindgen(641))*10000l,1,641),2,641)
fi_hf_4 = rebin(reform(400l + lindgen(4),4,1),4,641)+rebin(reform((1L+lindgen(641))*10000l,1,641),4,641)
fi_hf_8 = rebin(reform(800l + lindgen(8),8,1),8,641)+rebin(reform((1L+lindgen(641))*10000l,1,641),8,641)
fi_hf = [reform(fi_hf_1,641),reform(fi_hf_2,641*2),reform(fi_hf_4,641*4),reform(fi_hf_8,641*8)]
fi_h1 = fi_hf+30000000l
fi_h2 = fi_hf+40000000l
fi = [fi_a,fi_b,fi_c,fi_h1,fi_h2]
ff = fi_freq(fi)
off = sort(ff)
fi = fi(off)
ff = ff(off)
ff = uniq_list(ff)
nf = n_elements(ff)
if keyword_set(harmonic) then ff = ff/2
iso_fc_tmp = fltarr(3,nn,nf)
if verbose then plot,[0],/iso,psym=3,yr=[0,10],xr=[-10,10]
for i=0,nn-1 do begin
ww = where(fc(i,*) ne 0, cnt_fc)
if verbose then print,i,cnt_fc
if cnt_fc gt 1 then begin
iso_fc_tmp(0,i,*) = interpol(rtp(0,i,ww),fc(i,ww),ff)
iso_fc_tmp(1,i,*) = interpol(rtp(1,i,ww),fc(i,ww),ff)
ww1 = where(ff gt max(fc(i,ww)) or ff lt min(fc(i,ww)),cnt1,complement=ww2,ncomplement=cnt2)
if cnt1 ne 0 then iso_fc_tmp(*,i,ww1)=-1
if verbose and cnt2 ne 0 then oplot,iso_fc_tmp(0,i,ww2),iso_fc_tmp(1,i,ww2)/!radeg,/polar,psym=3
endif
endfor
iso_fc = ptrarr(nf,/allocate)
iso_fc_foot = ptrarr(nf,/allocate)
for i=0l,nf-1 do begin
ww = where(iso_fc_tmp(0,*,i) gt 0,cnt)
if cnt gt 0 then begin
stop
*iso_fc(i) = reform(iso_fc_tmp(*,ww,i),3,cnt)
*iso_fc_foot(i) = fltarr(3,cnt)
(*iso_fc_foot(i))(0,*) = radius_foot(ww)
(*iso_fc_foot(i))(1,*) = theta_foot(ww)
endif else begin
ptr_free,iso_fc(i)
ptr_free,iso_fc_foot(i)
endelse
endfor
file_out = getenv('ROOT_RPWS') + '/bin/'
file_out += 'iso_fc_grid_'+mfl_model
if keyword_set(nocrt) then file_out += '_nocrt'
if keyword_set(harmonic) then file_out += '_harmonic'
file_out +='.sav'
ff_iso_fc = ff
if keyword_set(save) then save,file=file_out,ff_iso_fc,iso_fc,iso_fc_foot
return
end