;+ ; Contains the build_iso_fc_grid procedure ; and supporting functions ; ; :Author: ; Baptiste Cecconi ; ; :History: ; 2009/05/28: Created ; ; 2009/05/28: Last Edit ;- ; ; derivative function ;+ ; computes the tangent (sph coord)) to the local magnetic field line ; at the RTP point ;- 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 ; ------------------------------------------------------------------------------ ;+ ; build_iso_fc_grid is a procedure that <behavior desc here> ; ; :Uses: ; fi_freq, magnetic_field ; ; :Params: ; ff: in, required, type=sometype ; A parameter named ff ; iso_fc: in, required, type=sometype ; A parameter named iso_fc ; iso_fc_foot: in, required, type=sometype ; A parameter named iso_fc_foot ; ; :Keywords: ; save: in, optional, type=sometype ; A keyword named save ; verbose: in, optional, type=sometype ; A keyword named verbose ; nocrt: in, optional, type=sometype ; A keyword named nocrt ; harmonic: in, optional, type=sometype ; A keyword named harmonic ;- 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 ; converts Gauss to kHz eps=0.098 ; planetary flattening h0=0.02 ; step in planetary radius fmin=3. ; low freq limit in kHz if keyword_set(harmonic) then fmin /=2. nn = 1801 ; number of points in latitudinal grid lat_pole = 0.2 ; degree 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) ;h = fltarr(3,nn)+h0; = h0*(2*(theta_foot ge 90)-1) rtp(0,*,0) = radius_foot rtp(1,*,0) = theta_foot rtp(2,*,0) = 0 rtp_i = rtp(*,*,0) ;testing sign of h DERIV,rtp_i,drtp_i,btot_i,bvec_i,mfl_model=mfl_model,nocrt=nocrt ; calcul de la derive drtp/ds 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));,arret 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