pro Multidata_Ions_Analysis close,/all ;--- 13 March 2003 ------- Merged analysis and plot programme ---------------- ;Part A: Read the Helios data file ;Part B: Determine a restricted array of the 3-D VDF from (12,12,32) ;Part C: Determine proton distribution and separate alphas ;Part D: Determine alpha particle distribution ;Part E: Calculate ion moments ;Part F: Plot the measured datap ;+++++++++++++++++++++++++++++++++++ ;Part A: Read the Helios data file ;+++++++++++++++++++++++++++++++++++ ; ; read the data files produced by Mr. Xian-zhi Ao in 2002 with ; Scherer's program ; ; programme to read the Helios ion/electron velocity distribution ; functions, vdfp, vdfe and the corresponding velocity components. ; program for the 3-D proton (2-D electron) velocity distribution function ; and determination of the velocity ranges and maximum vdf value ;do i = 1, 2 !sc ; do k = 1, 2 ! perishift ; do l = 1, 2 ! shift ; do m = 1, 32 ! energy ; do n = 1, 9 ! elevation ; do o = 1, 16 !azimuth ; vx(o,n,m,l,k,i) = volt1a(m,l,i) ; 1 * cos(azimuth(o,l,k,i)) * cos(elevation(n,i)) ; vy(o,n,m,l,k,i) = volt1a(m,l,i) ; 1 * sin(azimuth(o,l,k,i)) * cos(elevation(n,i)) ; vz(o,n,m,l,k,i) = volt1a(m,l,i) ; 1 * sin(elevation(n,i)) ; end do ; end do ; end do ; end do ; end do ; end do ; V_x radial ; V_y ; V_z northward out of the ecliptic plane ;------------------- variables and arrays ------------------------------------ comment=strarr(1) vdfp =dblarr(12,12,32) ; protons/alpha-particles 3-D cubes countp =intarr(12,12,32) vpolar_1=dblarr(32) azi_1 =dblarr(12) el_1 =dblarr(12) vdfe =dblarr(8,16) ; electrons 2-D squares vxe =dblarr(8,16) vye =dblarr(8,16) counte=intarr(8,16) l = 12*12*32 ;= 4608 maximum proton array size ;print, 'lp=', l l = 8*16 ;= 128 maximum electron array size ;print, 'le=', l ndimp = 4608 ndime = 128 ;---------------------- definitions of variables ------------------------------ ; ; iearcarrot = Carrington rotation number from Earth ; ihelcarrot = Carrington rotation number from spacecraft ; i1a_3_on = instrument i1a on (=0) or i3 on (=1) ; index_shift ; minus = indicates a HDM file which contained bad ; data (frames), but could be handled as NDM file ; mode = HDM or NMD ; iperihelion_shift = on or of (0 or 1) ; idummy, iyesno = dummy variable ; indexi1. = position of maximum value of instrument ; ;integer, ianzahl_e, ianzahl_p, idim, iearcarrot, ihelcarrot, i1a_3_on ;integer, index_shift, minus, mode ;integer, iperihelion_shift, idummy, ivesno ;integer, index_i1a(1), index_i1b(1) ; ; earcarlat,earcarlng = Carrington longitude and latitude from Earth ; earsdissun = distance Earth-Sun (always zero, historical) ; helcarlat,helcarlng = Carrington longitude and latitude from Helios ; heldissun = distance Helios S/C to Sun ; vnorm, vrad = S/C tangential (in ecliptic), radial velocity(AU/day) ; angle = spacecraft-Sun-Earth angle ; np = proton density [cm^{-3}] ; vxs,vys,vzs = center of mass velocity (1-fluid velocity) ; x,y,z,zfp,zfe = dummy parameters used to read in the 3-D VDF ; sinfi, cosfi = sine and cosine of azimuth of proton velocity ; sintheta, costheta = sine and cosine of elevation of proton velocity ; background_i1a,b = background of instruments (here the first two ; and last two channels) ; ;real, earcarlat, earcarlng, earsdissun ;real, dissun, carlat, carlng, vnorm, vrad, angle ;real, np, vxs, vys, vzs, x, y, z, zfp, zfe ;real, sintheta, costheta, sinfi, cosfi, background_i1a, background_i1b ; b = 3 components of magnetic field [nT] ; bsig = stanadrad deviation of the 3 components ; of magnetic field [nT] ; i1aint = 32 onboard integrated energy channels of ; protons insrtument i1a ; i1bint = 32 onboard integrated energy channels of ; protons insrtument i1b ; i1av, i1bv = corresponding velocities to i1aint, i1bint ; i1aalp = 3 dimenional array: density [cm^{-3}], ; velocity [km/s] and temperature [K] of alphas ; i1apro = 5 dimenional array: density [cm^{-3}], ; velocity [km/s], temperature [K] of protons ; i1bpro = 3 dimenional array: density [cm^{-3}], ; velocity [km/s], temperature [K] of protons b =dblarr(3) bsig =dblarr(3) i1aint=dblarr(32) i1av =dblarr(32) i1bint=dblarr(32) i1bv =dblarr(32) i1apro=fltarr(5) i1aalp=fltarr(3) i1bpro=fltarr(3) ; ;++++++++++++++++++++ Open and read data files +++++++++++++++++++++++++++++++++++++ openw,4,'anisotropy.dat' openw,5,'betapa.dat' openw,6,'wind_vel.dat' for da=105,107 do begin day=strtrim(string(da),1) ;openr,15, 'gooddata_old.txt' openr,15,'~/Desktop/Analysis_Helios/h276(017-366)/Y76D'+day+'/gooddata.txt' ;openr,15, 'C:\Dokumente und Einstellungen\Marsch\Eigene ;Dateien\Helios\Filenames\H2_76_107.txt' while (not eof(15)) do begin fnametemp='' ;' Y76D107/h2y76d107h00m20s31_hdm.0' readf,15, fnametemp ;print, fnametemp filen=strmid(fnametemp,1,33) print,filen fname=strmid(fnametemp,9,18) ;print, fname ;strdata=string ;strmom =string ;strplot=string strdata='ion_vdf_data_'+fname+'.dat' ;print,strdata strmom='ion_moments_'+fname+'.dat' ;print,strmom strplot='ion_plots_'+fname+'.ps' ;print,strplot ;++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ; Open output file ;openr,1,'C:\Dokumente und Einstellungen\marsch\Eigene Dateien\Helios\'+filen openr,1,'~/Desktop/Analysis_Helios/h276(017-366)/'+filen ;;; Lorenzo ;;; ;openw,2,strdata ;openw,3,strmom ;printf,2,'*********' ; now read in the data. The first data set, the correction from ; Pizzo et al. and for the finite solar diameter will not be used ; (the first three records). ; ; All other parameters are read in. For a definition see top. ; ; A 'read' without a parameter reads of a 'comment' in the data file ; ;print, 'enter reading' ; ;--------------------------------------------------------------------------- readf,1, comment ;print, comment readf,1, comment ;print, comment readf,1, comment ;print, comment readf,1, imode, ishift, iperihelion_shift, minus, i1a_3_on ;print, imode, ishift, iperihelion_shift, minus, i1a_3_on readf,1, heldissun, helcarlng, helcarlat, helcarrot ;print, heldissun, helcarlng, helcarlat, helcarrot readf,1, eardissun, earcarlng, earcarlat, angle, earcarrot ;print, eardissun, earcarlng, earcarlat, angle, earcarrot readf,1, vrad, vnorm ;print, vrad, vnorm readf,1, i1apro ;print,'density [cm^{-3}],velocity [km/s], temperature [K] of protons' ;print, i1apro readf,1, i1aalp ;print, i1aalp readf,1, i1bpro ;print,'i1bpro', i1bpro readf,1, b, bsig ;print,'b=', b, bsig readf,1, comment ;print, comment readf,1, i1aint ;print, i1aint readf,1, comment ;print, comment readf,1, i1av ;print, i1av readf,1, comment ;print, comment readf,1, i1bint ;print, i1bint readf,1, comment ;print, comment readf,1, i1bv ;print, i1bv readf,1, comment ;print, comment readf,1, comment ;print, comment readf,1, comment ;print, comment ;print, 'Integrated data' ;print, i1aint ;print, i1av ;print, i1bint ;print, i1bv bx=b(0) by=b(1) bz=b(2) ;printf,3, 'bx,by,bz',bx,by,bz ;printf,2,'b=',bx,by,bz ;--------------------------------------------------------------------------- ; Define velocity with respect to center of mass (1-D velocity). ; Transform from degree to radian. Calculate center of mass ; velocity 'v.s' ; pi180 = 3.141592653589793238462643d0/180.d0 ;print, pi180 ; i1apro(3),i1apro(4) ( in degree),sin( radians) sinfi = sin(i1apro(3) * pi180) cosfi = cos(i1apro(3) * pi180) sintheta = sin(i1apro(4) * pi180) costheta = cos(i1apro(4) * pi180) vxs = i1apro(1) * cosfi * costheta vys = i1apro(1) * sinfi * costheta vzs = i1apro(1) * sintheta ;printf,3, 'solar wind (proton) velocity components' ;printf,3, 'vxs,vys,vzs',vxs,vys,vzs ;print, 'solar wind (proton) velocity components' ;print, 'vxs,vys,vzs',vxs,vys,vzs ;---------------------------- Read VDFs -------------------------------------- ; ; now loop through the distribution fucntion. ; First read in the indices i,j,k of distribution function and the velocity ; components. ; The indices are: i for azimuth, ; j for elevation ; k for voltage. ; The dummy variables zfp,zfe the distribution function, and ; z,x,y, the respective velocity components. ; The variable 'dummy' the count rate to be used for statistics. ; The dummy variables are then transformed in the center of mass system ; and written to ; the velocity arrays : wx,wy,wz. ; If an error occurs during read, the program will continue at label electrons. ; The next data will be the electrons, which are handled later. ; ;========================= protons =========================================== for mi=0,11 do begin for mj=0,11 do begin for mk=0,31 do begin vdfp(mi,mj,mk) = 0.d0 countp(mi,mj,mk) = 0 vpolar_1(mk)=0.d0 azi_1(mi)=0.d0 el_1(mj)=0.d0 endfor endfor endfor ;printf,2, 'input kk, ii, jj, dummy, zfp (cm/s)^-3 (cm)^-3' ;printf,2, 'vpolar_s (km/s),azi_s (radians),el_s(radians)' for n_p=0, ndimp do begin ; the loop parameter (after for) atomatically add 1 after one cicle on_ioerror, electrons ; jump to electron data readf,1, i, j, k, zfp, dummy, x, y, z,vpolar_s,azi_s,el_s ;print,i, j, k, zfp, dummy, x, y, z,vpolar_s,azi_s,el_s i=fix(i) j=fix(j) k=fix(k) dummy = fix(dummy) ii= i-1 jj= j-1 kk= k-1 ;printf,2, kk, ii, jj, dummy, zfp,vpolar_s,azi_s,el_s,format='(3i3,1i4,2x,e15.8,3f10.3)' vdfp(ii,jj,kk) = zfp countp(ii,jj,kk)= dummy vpolar_1(kk)=vpolar_s azi_1(ii)=azi_s el_1(jj)=el_s endfor ;============================ electrons =========================== ; ; Now read electron data. Procedure similar to that of protons. ; First print one coment line, then loop: Structure: i,j indices ; of the distribution function, zfe value of distribution function, ; dummy is the count rate, x,y are the electron velocities electrons1: readf,1, comment ;print, comment ; readf,1, comment ; print, comment for mi=0,7 do begin for mj=0,15 do begin vdfe(mi,mj) = 0.d0 vxe(mi,mj) = 0.d0 vye(mi,mj) = 0.d0 counte(mi,mj) = 0 endfor endfor for n_e=0, ndime do begin ; with the 'while...', the loop parameter n_e is constant 0. while not EOF(1) do begin readf,1, i, j, zfe, dummy, x, y i = fix(i) j = fix(j) dummy = fix(dummy) ii = i - 1 jj = j - 1 ;print, 'electrons ', ii, jj, dummy,' ', zfe, x, y vdfe(ii,jj) = zfe vxe(ii,jj) = x vye(ii,jj) = y counte(ii,jj) = dummy ; ;printf,2,'n_e=',n_e endwhile endfor electrons: ;close,1 ;----------------- Output spherical coordinates ---------------------------------- ;printf,3,'azi_1,i' for i=0, 11 do begin ;printf,3,azi_1(i),i,format='(f10.3,i5)' endfor ;printf,3,'el_1(j),j' for j=0, 11 do begin ;printf,3,el_1(j),j,format='(f10.3,i5)' endfor ;printf,3,'vpolar_1(k),k' for k=0, 31 do begin ;printf,3,vpolar_1(k),k,format='(f10.3,i5)' endfor ;---------------------- end reading the data from file -------------------- ;++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ;Part B: Determine a restricted array of the 3-D VDF from (12,12,32) ;+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ;------------ Determine the velocity ranges and maxima --------------------- ;------------ Limit the domain of the ion and electron --------------------- ;------------ velocity distribution functions ------------------------------ ;------------ The arrays are defined above --------------------------------- ;------------- as determined in the read subroutine ------------------------ ;------ Determine the velocity domain with VDF above 1-2-count level ------- ;------ Determine proton indices of the new boundaries --------------------- ptoti=intarr(12) ptotj=intarr(12) ptotk=intarr(32) for k=0,31 do begin ptotk(k) = total(countp(*,*,k)) ;print, k, ptotk(k) endfor for j=0,11 do begin ptotj(j) = total(countp(*,j,*)) ;print, j, ptotj(j) endfor for i=0,11 do begin ptoti(i) = total(countp(i,*,*)) ;print, i, ptoti(i) endfor ; electrons etoti=intarr(8) etotj=intarr(16) for i=0,7 do begin etoti(i) = total(counte(i,*)) ;print, i, etoti(i) endfor for j=0,15 do begin etotj(j) = total(counte(*,j)) ;print, j, etotj(j) endfor ;print, 'determine proton array limits' ;Set the minimum count rate, mk, mj, mi, to determine the restricted range ;------------------------------------------------------------------------- ; klow, khigh, jlow, jhigh, ilow, ihigh ;---------------------------------------- mk=3 ;15 ;10 ;5 mj=2 ;10 ;2 mi=2 ;10 ;2 for k=0,31 do begin if ptotk(k) GT mk then GOTO, continue_k1 endfor continue_k1: ;print, k, ptotk(k) klow = k for k=0,31 do begin kk = 31-k if ptotk(kk) GT mk then GOTO, continue_k2 endfor continue_k2: ;print, kk, ptotk(kk) khig = kk for j=0,11 do begin if ptotj(j) GT mj then GOTO, continue_j1 endfor continue_j1: ;print, j, ptotj(j) jlow = j for j=0,11 do begin jj = 11-j if ptotj(jj) GT mj then GOTO, continue_j2 endfor continue_j2: ;print, jj, ptotj(jj) jhig = jj for i=0,11 do begin if ptoti(i) GT mi then GOTO, continue_i1 endfor continue_i1: ;print, i, ptoti(i) ilow = i for i=0,11 do begin ii = 11-i if ptoti(ii) GT mi then GOTO, continue_i2 endfor continue_i2: ;print, ii, ptoti(ii) ihig = ii ;print, 'ilow=', ilow, ' ihig=', ihig, ' jlow=', jlow, ' jhig=', jhig, $ ;' klow=', klow, ' khig=', khig ;------ Determine the restricted array size (idimen,jdimen,kdimen) -------------- idimen = ihig - ilow + 1 jdimen = jhig - jlow + 1 kdimen = khig - klow + 1 number= idimen*jdimen*kdimen ;print, 'idimen=', idimen, ' jdimen=', jdimen, ' kdimen=', kdimen ;--- The arrays dimensions may be negative in some cases, when the count rates are ---- ;--- zero which happens in some odd spectra. Then jump over this spectrum! ---------- if idimen lt 0 then goto, continue_spectrum if jdimen lt 0 then goto, continue_spectrum if kdimen lt 0 then goto, continue_spectrum fp=dblarr(idimen,jdimen,kdimen) ; protons/alpha-particles reduced 3-D cubes cp=intarr(idimen,jdimen,kdimen) vpolar_r=dblarr(kdimen) azi_r =dblarr(idimen) el_r =dblarr(jdimen) m0=0 ; the number of data gaps for i=0, idimen-1 do begin ii = ilow + i for j=0, jdimen-1 do begin jj = jlow + j for k=0, kdimen-1 do begin kk = klow + k fp(i,j,k) = 0.d0 cp(i,j,k) = 0 vpolar_r(k)=0.d0 azi_r(i) =0.d0 el_r(j) =0.d0 ;---------------------------------------------------------------------------------- fp(i,j,k) = vdfp(ii,jj,kk) cp(i,j,k) = countp(ii,jj,kk) vpolar_r(k)=vpolar_1(kk) azi_r(i) =azi_1(ii) el_r(j) =el_1(jj) continue_10: endfor endfor endfor ;------------------------------------- r0=float(m0)/float(number) ;print,m0,r0 ;------------------------------------- cptot=dblarr(kdimen) faint=dblarr(kdimen) va =dblarr(kdimen) fbint=dblarr(kdimen) vb =dblarr(kdimen) ;--------------------- Reduce the size of the integrated spectra -------------- ratiomax = max(i1aint)/max(i1bint) ;print, 'ratiomax= ', ratiomax for k=0,kdimen-1 do begin kk = klow + k cptot(k)= ptotk(kk) faint(k)= i1aint(kk) va(k) = i1av(kk) fbint(k)= i1bint(kk)*ratiomax ; to set the same peak value vb(k) = i1bv(kk) ;print,faint(k),va(k),fbint(k),vb(k) endfor ;----------- Determine the summed-up (over i and j) VDF of protons ----------- ;fsum=dblarr(kdimen) ;for k=0,kdimen-1 do begin ; fsum(k) = total(fp(*,*,k)) ; ;print, k, fsum(k) ;endfor ;print, faint ;print, va ;print, fbint ;print, vb ;++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ;Part C: Determine proton distribution and separate alphas ;++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ;print, 'enter alpha separation' ;------------- find the separation index in the I1A/B ratio spectrum -------------- i1amax=max(faint,kamx) i1bmax=max(fbint,kbmx) ;print,' maxa ', i1amax, kamx, ' maxb ', i1bmax, kbmx if (kbmx GT kamx) then maxk=kbmx else maxk=kamx intratio = dblarr(32) intratio = alog10(i1bint)-alog10(i1aint) ;print, intratio intrmax = intratio(maxk+klow) ;print, 'ratio ', intrmax, ' at fabint maximum index', maxk minrat= min(intratio,imin) mink=imin-klow ;print, 'ratiomin ', minrat, ' mink ', mink if (maxk GT mink) then kstart=maxk else kstart=mink ;print, 'kstart ', kstart ;------------------------------------------------------------------------------------ knoal = kdimen ratcrit = alog10(2) ;print, ratcrit kbegin= klow + kstart khignoa = kbegin ;print, 'kbegin =', kbegin, 'khignoa =', khignoa, 'khig =', khig for kk = kbegin,khig do begin khignoa = kk ;print, 'ratio ', intratio(kk), kk if (intratio(kk) GE ratcrit) then goto, continue endfor continue: knoal = khignoa - klow + 2 k_min = knoal - 2 if (knoal LT kdimen) then goto, continue_ions else knoal = kdimen k_min = kdimen print, 'no alphas separated' ;print,'knoal= ', knoal,' k_min=', k_min goto, continue_protons continue_ions: ;print,'proton upper and alpha lower limits found' ;print,'knoal= ', knoal,' k_min=', k_min ;-----------out put information and data for interpolation ----------------- continue_protons: ;openw, 10,'proton_vdf_3d_'+fname+'.dat' ;printf,10,idimen,jdimen,knoal,format='(3I8)' ;printf,10,vxs,vys,vzs,bx,by,bz,format='(6f15.8)' ;-------------- find fpmax=fp(imax,imax,kmax) ------------------------------ fpmax=0 imax=0 jmax=0 kmax=0 for k=0,kdimen-1 do begin for j=0,jdimen-1 do begin for i=0,idimen-1 do begin if (fpmax lt fp(i,j,k) ) then begin fpmax=fp(i,j,k) imax=i jmax=j kmax=k endif endfor endfor endfor COS_A_m = COS(azi_r(imax)) COS_E_m = COS(el_r(jmax)) SIN_A_m = SIN(azi_r(imax)) SIN_E_m = SIN(el_r(jmax)) vx_m = vpolar_r(kmax)*COS_A_m*COS_E_m vy_m = vpolar_r(kmax)*SIN_A_m*COS_E_m vz_m = vpolar_r(kmax)*SIN_E_m ;printf,10,fp(imax,jmax,kmax),imax,jmax,kmax,format='(e15.8,3I8)' ;print, fp(imax,jmax,kmax),imax,jmax,kmax,format='(e15.8,3I8)' ;------------------reduce 3-D proton cube ------------------------------- fp_noa =dblarr(idimen,jdimen,knoal) cp_noa =intarr(idimen,jdimen,knoal) vpol =dblarr(knoal) azimuth =dblarr(idimen) elevation=dblarr(jdimen) faint_noa=dblarr(knoal) fbint_noa=dblarr(knoal) va_noa =dblarr(knoal) vb_noa =dblarr(knoal) for k1=0,knoal-1 do begin k=k1 vpol(k) = vpolar_r(k1) faint_noa(k)= faint(k1) fbint_noa(k)= fbint(k1) va_noa(k) = va(k1) vb_noa(k) = vb(k1) for j=0,jdimen-1 do begin for i=0,idimen-1 do begin fp_noa(i,j,k)= fp(i,j,k1) cp_noa(i,j,k)= cp(i,j,k1) azimuth(i) = azi_r(i) elevation(j) = el_r(j) ;;printf,2, k, i, j, dummy,fp_noa,vpol(k),azimuth (i),elevation(j),format='(4i3,2x,e15.8,3f10.3)' ;printf,10,fp_noa(i,j,k),vpol(k),azimuth(i),elevation(j),i+1,j+1,k+1,format='(e15.8e2,3d15.8,3i8)' endfor endfor endfor close,10 ;printf,2, 'vpol', vpol ;printf,2, 'azimuth', azimuth ;printf,2, 'elevation', elevation ;printf,2, 'fpnoa',fp_noa ;print,'protons' ;print,'j,elevation',elevation ;print, 'vpol',vpol ;+++++++++++++++++++++++++++++++++++++++++++++++++++ ;Part D: Determine alpha particle distribution ;+++++++++++++++++++++++++++++++++++++++++++++++++++ ; Distribution function for alpha particles: f_a ; ; The part of fp (idimen,jdimen,kdimen) with k ranging between kmin ; and (kdimen-1), is considered as belonging to the alpha particles. ; The velocity V must be divided by sqrt(2) and the VDF fp multiplied by a factor of 4 ; to get alpha velocity distribution f_a(idimen,jdimen,k_a), where k_a =kdimen-kmin. ;--------------------------------------------------------------------------- continue_alphas: k_a= kdimen - k_min ;print, 'no alphas ', 'k_a= ', k_a ;print,'ciao-2' if (k_a EQ 0) then goto, continue_end ;print,'ciao-1' ;openw, 11,'alpha_vdf_3d_'+fname+'.dat' ;information and data for interpolation ;printf,11,idimen,jdimen,k_a,format='(3I8)' ;printf,11,vxs,vys,vzs,bx,by,bz,format='(6f15.8)' f_a = dblarr(idimen,jdimen,k_a) c_a = intarr(idimen,jdimen,k_a) v_a = dblarr(k_a) azi_a = dblarr(idimen) ele_a = dblarr(jdimen) for k=0,k_a-1 do begin kk = k + k_min for j=0,jdimen-1 do begin for i=0,idimen-1 do begin f_a(i,j,k)=fp(i,j,kk)*4. c_a(i,j,k)=cp(i,j,kk) v_a(k)=vpolar_r(kk)/sqrt(2.) azi_a(i)=azi_r(i) ele_a(j)=el_r(j) endfor endfor endfor ;------------------ find the maximum f_a=f_a(im_a,j_ma,km_a) -------------------- fm_a=0 im_a=0 jm_a=0 km_a=0 for k=0,k_a-1 do begin for j=0,jdimen-1 do begin for i=0,idimen-1 do begin if (fm_a lt f_a(i,j,k) ) then begin fm_a=f_a(i,j,k) im_a=i jm_a=j km_a=k endif endfor endfor endfor ;printf,11,f_a(im_a,jm_a,km_a),im_a,jm_a,km_a,format='(e15.8,3I8)' for k=0,k_a-1 do begin for j=0,jdimen-1 do begin for i=0,idimen-1 do begin ;printf,11,f_a(i,j,k),v_a(k),azi_a(i),ele_a(j),i+1,j+1,k+1,format='(e15.8e2,3d15.8,3i8)' endfor endfor endfor close,11 ;------ end of putting out the distribution of alpha particles: f_a ------------ ;+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ; Part E: Calculate ion moments ;+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ; ; m_spc=1 for protom moments (spc species) ; k_spc= knoal ; f_spc= fp_noa ; vo_spc= vpol ; az_spc= azimuth ; el_spc=elevation ; ; m_spc=2 for alpha moments ; k_spc= k_a ; f_spc= f_a ; vo_spc=v_a ; az_spc=azi_a ; el_spc=ele_a ; ; jdimen and kdimen do not change ;---------------------------------------------------------------------------------- ;Number density n_d=fltarr(2) ;print,'ciao0' ;The velocity at the observed maximum VDF_spc in the satellite frame vx_max=fltarr(2) vy_max=fltarr(2) vz_max=fltarr(2) ;The group velocity corrected to the v_max by moment integration vc_x=fltarr(2) vc_y=fltarr(2) vc_z=fltarr(2) ;The final corrected group velocity in the satellite frame vs_x=fltarr(2) vs_y=fltarr(2) vs_z=fltarr(2) ; ;The mean thermal velocity in the v_max frame, inlcuding the corrections vox,voy,voz ;The parallel thermal velocity ;The perpendicular thermal velocity v_th_t=fltarr(2) v_ll =fltarr(2) v_perp=fltarr(2) ;The perpendicular temperature ;The parallel temperature ;The Anisotropy mu =fltarr(2) T_ll =fltarr(2) T_perp=fltarr(2) T_t =fltarr(2) A =fltarr(2) fact =fltarr(2) ; factor to set the lower limit of integration domain f_max =fltarr(2) qx =fltarr(2) ; heat flux components qy =fltarr(2) qz =fltarr(2) vx_max(0)=vx_m vy_max(0)=vy_m vz_max(0)=vz_m COS_A = COS(azi_a(im_a)) SIN_A = SIN(azi_a(im_a)) COS_E = COS(ele_a(jm_a)) SIN_E = SIN(ele_a(jm_a)) vx_max(1)=va(km_a)*COS_A*COS_E vy_max(1)=va(km_a)*SIN_A*COS_E vz_max(1)=va(km_a)*SIN_E mu(0) = 1. mu(1) = 4. fact(0) = 0.001 fact(1) = 0.03 f_max(0)= fpmax f_max(1)= fm_a ;----- Calculate first proton moments and then alpha moments -------------------- for m_spc = 1,2 do begin ;print,'ciao1' if (m_spc eq 1) then begin ;print,'ciao2' k_spc= knoal f_spc=dblarr(idimen,jdimen,k_spc) vo_spc=dblarr(k_spc) az_spc=dblarr(idimen) el_spc=dblarr(jdimen) for i=0, idimen-1 do begin for j=0, jdimen-1 do begin for k=0, k_spc-1 do begin f_spc(i,j,k)= fp_noa(i,j,k) if (f_spc(i,j,k) lt f_max(m_spc-1)*fact(m_spc-1)) then f_spc(i,j,k)=0. vo_spc(k)= vpol(k) az_spc(i)= azimuth(i) el_spc(j)= elevation(j) endfor endfor endfor endif ;---------------------- alpha moments ------------------------------------ ; ; k_spc= k_a ; f_spc= f_a ; vo_spc=v_a ; az_spc=azi_a ; el_spc=ele_a if (m_spc eq 2) then begin ;print,'ciao3' k_spc= k_a if (k_spc eq 0) then goto, continue_end f_spc=dblarr(idimen,jdimen,k_spc) vo_spc=dblarr(k_spc) az_spc=dblarr(idimen) el_spc=dblarr(jdimen) for i=0, idimen-1 do begin for j=0, jdimen-1 do begin for k=0, k_spc-1 do begin f_spc(i,j,k)= f_a(i,j,k) if (f_spc(i,j,k) lt f_max(m_spc-1)*fact(m_spc-1)) then f_spc(i,j,k)=0. vo_spc(k)= v_a(k) az_spc(i)= azi_a(i) el_spc(j)=ele_a(j) endfor endfor endfor endif ;++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ; ; for core proton only (m_spc=1) ; . ;++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ; calculate proton velocity components ;--------------------------------------------------------------- ; in the frame of v at maximum vdf wx=dblarr(idimen,jdimen,k_spc) wy=dblarr(idimen,jdimen,kdimen) wz=dblarr(idimen,jdimen,kdimen) ;---- in the satellite frame ------------------------------------ xp=dblarr(idimen,jdimen,k_spc) yp=dblarr(idimen,jdimen,kdimen) zp=dblarr(idimen,jdimen,kdimen) nt=0 for i=0, idimen-1 do begin for j=0, jdimen-1 do begin for k=0, k_spc-1 do begin ;-------------- az,el in radans, sin(radians) --------------------- COS_A = COS(az_spc(i)) COS_E = COS(el_spc(j)) SIN_A = SIN(az_spc(i)) SIN_E = SIN(el_spc(j)) xp(i,j,k) = vo_spc(k)*COS_A*COS_E yp(i,j,k) = vo_spc(k)*SIN_A*COS_E zp(i,j,k) = vo_spc(k)*SIN_E ;----- Transform velocity to the maximum frame at velocity: v_x_max,v_x_max,v_x_max wx(i,j,k) = xp(i,j,k) - vx_max(m_spc-1) wy(i,j,k) = yp(i,j,k) - vy_max(m_spc-1) wz(i,j,k) = zp(i,j,k) - vz_max(m_spc-1) endfor endfor endfor ;------------ calculate moments by simple summation ----------------------------- ;print,'----------------- enter moments ---------------------------------------' pn =0. v_x=0. v_y=0. v_z=0. d_el =abs(el_spc(jdimen-1)-el_spc(0))/float(jdimen-1) d_az =abs(az_spc(idimen-1)-az_spc(0))/float(idimen-1) ;print, 'd_el=', d_el, ' d_az=', d_az for i=0, idimen-1 do begin for j=0, jdimen-1 do begin for k=0, k_spc-1 do begin IF k EQ 0 THEN d_v=vo_spc(1) - vo_spc(0) ELSE $ IF k EQ k_spc-1 THEN d_v=vo_spc(k_spc-1)- vo_spc(k_spc-2) ELSE $ d_v=(vo_spc(k+1)-vo_spc(k-1))/2. ;print,'d_v', d_v unit=(1.e5) ; from km to cm dn= f_spc(i,j,k)*vo_spc(k)*vo_spc(k)*cos(el_spc(j))*d_v*d_el*d_az*unit^3 pn = pn + dn v_x= v_x + wx(i,j,k)*dn v_y= v_y + wy(i,j,k)*dn v_z= v_z + wz(i,j,k)*dn endfor endfor endfor n_p=i1apro(0) ;print,'density=', n_p ;print, 'moments', pn,v_x/pn,v_y/pn,v_z/pn n_d(m_spc-1)=pn vc_x(m_spc-1) = v_x/pn vc_y(m_spc-1) = v_y/pn vc_z(m_spc-1) = v_z/pn ;-------- Correct for the slight shift from maximum to real bulk frame ------------- vs_x(m_spc-1) = vx_max(m_spc-1) + vc_x(m_spc-1) vs_y(m_spc-1) = vy_max(m_spc-1) + vc_y(m_spc-1) vs_z(m_spc-1) = vz_max(m_spc-1) + vc_z(m_spc-1) ;----- magnetic field unit vector --------------------------------------------- e_b=b/sqrt(bx^2 + by^2 + bz^2) ;print, e_b w1=0. w2=0. w3=0. wt=0. q=[0,0,0] ; auxiliary heat flux vector e_x = [1,0,0] ; should always be in the direction: az=0, el=0 e_by = crossp (e_x , e_b) ; e_by = e_x X e_b e_bz = crossp (e_by, e_b) ; e_bz = e_by X e_b e_by=e_by/sqrt(e_by(0)^2 + e_by(1)^2 + e_by(2)^2) e_bz=e_bz/sqrt(e_bz(0)^2 + e_bz(1)^2 + e_bz(2)^2) unit=(1.e5) ; to change velocity units from km/s to cm/s vox=vc_x(m_spc-1) voy=vc_y(m_spc-1) voz=vc_z(m_spc-1) for i=0, idimen-1 do begin for j=0, jdimen-1 do begin for k=0, k_spc-1 do begin IF k EQ 0 THEN d_v= vo_spc(1)-vo_spc(0) ELSE $ IF k EQ k_spc-1 THEN d_v= vo_spc(k_spc-1)-vo_spc(k_spc-2) ELSE $ d_v= (vo_spc(k+1)-vo_spc(k-1))/2. dn= f_spc(i,j,k)*vo_spc(k)*vo_spc(k)*cos(el_spc(j))*d_v*d_el*d_az*unit^3 v_w=[(wx(i,j,k)-vox),(wy(i,j,k)-voy),(wz(i,j,k)-voz)] ; in the direction e_b c1= v_w*e_b w1= w1+ (c1(0)+c1(1)+c1(2))^2*dn ; in the direction e_by c2= v_w*e_by w2= w2+ (c2(0)+c2(1)+c2(2))^2*dn ; in the direction e_bz c3= v_w*e_bz w3= w3+ (c3(0)+c3(1)+c3(2))^2*dn ; total thermal speed ct= v_w*v_w wt= wt+ (ct(0)+ct(1)+ct(2))*dn q = q + (ct(0)+ct(1)+ct(2))*v_w*dn ; heat flux endfor endfor endfor qx(m_spc-1) = 0.5*q(0) qy(m_spc-1) = 0.5*q(1) qz(m_spc-1) = 0.5*q(2) v_th_by =sqrt(w2/pn) v_th_bz =sqrt(w3/pn) v_th_per=sqrt(0.5*(w2+w3)/pn) v_th_t(m_spc-1)= sqrt(wt/pn) v_ll(m_spc-1) = sqrt(w1/pn) v_perp(m_spc-1)= sqrt(0.5*(wt-w1)/pn) T_ll(m_spc-1) = mu(m_spc-1)*(v_ll(m_spc-1) /91.)^2 ; in unit of 10^6 degree T_perp(m_spc-1)= mu(m_spc-1)*(v_perp(m_spc-1)/91.)^2 T_t(m_spc-1) = mu(m_spc-1)*(v_th_t(m_spc-1)/91.)^2 A(m_spc-1) = (v_perp(m_spc-1)/v_ll(m_spc-1))^2 - 1. ;print,v_ll(m_spc-1),v_perp(m_spc-1),v_th_t(m_spc-1),v_th_by,v_th_bz,v_th_per ;if m_spc eq 1 then printf,3,' ########## Proton Moments ########### ' ;if m_spc eq 2 then printf,3,' ########## Alpha Moments ########### ' ;printf,3,' ' ;printf,3,'n_p(i1apro(0))=',n_p,'np=',n_d(m_spc-1) ;printf,3,' ' ;printf,3,'Maximum VDF_spc, correction of v_max, bulk velocity in satellite frame' ;printf,3,vx_max(m_spc-1),vc_x(m_spc-1),vs_x(m_spc-1) ;printf,3,vy_max(m_spc-1),vc_y(m_spc-1),vs_y(m_spc-1) ;printf,3,vz_max(m_spc-1),vc_z(m_spc-1),vs_z(m_spc-1) ;printf,3,' ' ;printf,3,'v_ll,..v_perp,...v_th_t,...v_th_by,...v_th_bz,...v_th_per' ;printf,3,v_ll(m_spc-1) ,v_perp(m_spc-1),v_th_t(m_spc-1),v_th_by,v_th_bz,v_th_per ;printf,3,' ' ;printf,3,'T_ll=',T_ll(m_spc-1) ,' T_perp=',T_perp(m_spc-1),' A=Tperp/T_ll-1=',A(m_spc-1) if (m_spc eq 1) then begin printf,4,A(m_spc-1) ;betapa=T_ll(m_spc-1)*n_d(m_spc-1)*4.03*10^(-11)*86.1/() ;printf,5,T_ll(m_spc) endif ;printf,3,'qx,qy,qz', q endfor ; end of calculation of proton moments and then alpha moments ;---------------------------------------------------------------------------- ;printf,3,' ########## Proton alpha comparison ####################' ;printf,3,' ' B_gaus = sqrt(bx^2 + by^2 + bz^2)*1.d-1*1.d-5 ; b in 0.1 nT v_A = 2.18d11*B_gaus*n_d(0)^(-0.5)/unit ; km/s, 1 gaus=1.d5 nT beta_ll = 2.* (v_ll(0)/v_A)^2 printf,5,beta_ll printf,6,vxs abundance = n_d(1)/n_d(0) ;--------- Alpha velocity relative to the proton velocity v_drift= [(vs_x(1)-vs_x(0)),(vs_y(1)-vs_y(0)),(vs_z(1)-vs_z(0))] v_d=sqrt(v_drift(0)^2 + v_drift(1)^2 + v_drift(2)^2) ;--------- Alpha velocity relative along magnetic direction c4= v_drift*e_b v_d_b=abs((c4(0)+c4(1)+c4(2))) ;--------- Angle between the alpha drift velocity and magnetic field ang_cos= (v_d_b/v_d) ang=acos(v_d_b/v_d) ;in radians ;--------- Ratio between the alpha total temperature and the proton temperature ratio_t= T_t(1)/T_t(0) ;print,v_A,beta_ll ;printf,3,'e_b=',e_b,' B (nT)=',B_gaus*1.e5 ;printf,3,'v_A=',v_A, ' beta_ll=',beta_ll ;printf,3,'alpha abundance= ',abundance, ' T_a_t/T_p_t= ',ratio_t ;printf,3,'v_d= ',v_d ,' v_d_b= ',v_d_b,' v_d_b/V_A=',v_d_b/V_A ;printf,3,'ang_pa(in degree)= ',ang*180./3.1416, ' v_d_b/V_d= ',v_d_b/V_d continue_end: ;--------------------- end of moments ------------------------------------------ ;print, '-------------------------------enter plotting---------------------------' ;++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ; Part F: Plot the measured data ;++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ; Calculate the distribution function (protons and alphas) in the plane with ; minimum azimuth angle (I=2); x-z plane, z northward, x radial ; ************************ arrays ************************************** wx_i=fltarr(kdimen,jdimen) wz_i=fltarr(kdimen,jdimen) wbx=fltarr(kdimen) wbz=fltarr(kdimen) fpz=fltarr(kdimen,jdimen) cpz=fltarr(kdimen,jdimen) wx_1=fltarr(idimen,jdimen,kdimen) wy_1=fltarr(idimen,jdimen,kdimen) wz_1=fltarr(idimen,jdimen,kdimen) for i=0, idimen-1 do begin for j=0, jdimen-1 do begin for k=0, kdimen-1 do begin ; azi,ei in radans, sin(radians) COS_A = COS(azi_r(i)) COS_E = COS(el_r(j)) SIN_A = SIN(azi_r(i)) SIN_E = SIN(el_r(j)) wx_1(i,j,k) = vpolar_r(k)*COS_A*COS_E - vxs wy_1(i,j,k) = vpolar_r(k)*SIN_A*COS_E - vys wz_1(i,j,k) = vpolar_r(k)*SIN_E - vzs endfor endfor endfor i= 2 for k=0,kdimen-1 do begin for j=0,jdimen-1 do begin COS_A = COS(azi_r(i)) COS_E = COS(el_r(j)) SIN_A = SIN(azi_r(i)) SIN_E = SIN(el_r(j)) vx_p = vpolar_r(k)*COS_A*COS_E vy_p = vpolar_r(k)*SIN_A*COS_E vz_p = vpolar_r(k)*SIN_E wx_i(k,j)=vx_p - vxs wz_i(k,j)=vz_p - vzs wbx(k)=wx_i(k,j) wbz(k)=wbx(k)*bz/bx fpz(k,j)=fp(i,j,k)/max(fp) cpz(k,j)=cp(i,j,k) endfor endfor ;****************** prints ************************************************** ;print, 'vxs=',vxs,' knoal=', knoal ;print,'ilow=',ilow,' ihig=',ihig, ' jlow=',jlow,' jhig=',jhig, ' klow=',klow, ' khig=',khig for m=0,idimen-1 do begin ;print,'azi_r(m)= ', azi_r(m) endfor for m=0,jdimen-1 do begin ;print,'el_r(m)= ', el_r(m) endfor if (3 eq 4) then begin ;**************************** plots *************************************** ;set_plot,'win' set_plot,'ps' device, filename = strplot, YSIZE = 10.0, XSIZE = 10.0 ;********************************************************************************* ; integrated energy spectrum with alphas included ;********************************************************************************* !p.multi=[0,1,2] ;!x.range=[300,1500] plot, va, alog10(faint), LINESTYLE=0, $ xstyle=1, xrange=[300,1500], xtitle='speed /(km/s)', ytitle='Integral VDF' oplot, vb, alog10(fbint), LINESTYLE=2 zero = dblarr(32) zero = zero-18.3 oplot, va, zero, PSYM=7 ;------------------------------ with alphas included --------------------------- ;plot, va-vxs, alog10(faint), LINESTYLE=0 ;oplot, vb-vxs, alog10(fbint), LINESTYLE=2 ;-------------------------------without alphas -------------------------------- ;plot, va_noa-vxs, alog10(faint_noa), LINESTYLE=0 ;oplot, vb_noa-vxs,alog10(fbint_noa), LINESTYLE=2 ;zero = dblarr(32) ;zero = zero-18.5 ;plot, i1av-vxs, alog10(i1aint), LINESTYLE=0 ;oplot, i1bv-vxs, alog10(i1bint), LINESTYLE=2 ;oplot, i1av-vxs, zero,PSYM=7 cforij=1 if (cforij eq 2) then begin contour, cpt,wxp,wyp,nlevels=10 ;,/ISOTROPIC plots,wxp,wyb endif if (cforij eq 1) then begin contour, cpz/max(cpz),wx_i,wz_i,levels=[0.01,0.032,0.1,0.2,0.4,0.6,0.8],$ C_LINESTYLE=[1,1,0,0,0,0,0],$ C_ANNOTATION=['0.01','0.032','0.1','0.2','0.4','0.6','0.8'], xrange=[-200,800], $ xtitle='Vx /(km/s)', ytitle='Vz /(km/s)' plots,wbx,wbz for j=0, jdimen-1 do begin for k=0, kdimen-1 do begin plots,wx_1(2,j,k),wz_1(2,j,k),psym=1 endfor endfor endif ;******************************************************************************* ; speed range [-200,800] for protons and alpha particles ; ;******************************************************************************* !p.multi=[0,1,2] maxfp=max(fp) ;print,'maxfp: ', maxfp ;print, imax, jmax, imax+1, jmax+1, imax-1, jmax-1 plot, alog10(fp(imax,jmax,*)/maxfp),XTICKS=kdimen-1,linestyle=0,xrange=[0,kdimen-1],$ xtitle='voltage channel number', ytitle='1-D VDFs' oplot,alog10(fp(imax+1,jmax,*)/maxfp),linestyle=1 oplot,alog10(fp(imax-1,jmax,*)/maxfp),linestyle=2 oplot,alog10(fp(imax,jmax+1,*)/maxfp),linestyle=3 oplot,alog10(fp(imax,jmax-1,*)/maxfp),linestyle=4 ;oplot,alog10(cptot/max(cptot)),colar=1,linestyle=2 ;plot, alog10(fp(1,2,*)/maxfp,XTICKS=kdimen-1,linestyle=0 ;oplot,alog10(fp(2,2,*)/maxfp,linestyle=1 ;oplot,alog10(cp(1,4,*)/max(cpz)),linestyle=0 ;oplot,alog10(cp(2,4,*)/max(cpz)),linestyle=1 ;plot, alog10(faint),XTICKS=kdimen-1,LINESTYLE=0 ;oplot,alog10(fbint),LINESTYLE=2 ; ;va-vxs, ;vb-vxs, ;!x.range=[-200,800] if (cforij eq 2) then begin contour, cpt,wxp,wyp,nlevels=10 ;,/ISOTROPIC plots,wxp,wyb endif maxcpz = max(cpz) ;print, 'maxcpz: ',maxcpz if (cforij eq 1) then begin contour, cpz/maxcpz,wx_i,wz_i,levels=[0.01,0.032,0.1,0.2,0.4,0.6,0.8],$ C_LINESTYLE=[1,1,0,0,0,0,0], xrange=[-200,800], $ xtitle='Vx /(km/s)', ytitle='Vz /(km/s)' ; C_ANNOTATION=['0.01','0.032','0.1','0.2','0.4','0.6','0.8'] plots,wbx,wbz for j=0, jdimen-1 do begin for k=0, kdimen-1 do begin plots, wx_1(2,j,k), wz_1(2,j,k), psym=1 endfor endfor endif ;------------------------------------------------------------------------------ !p.multi=[0,1,3] zero=fltarr(32) zero=zero-18.5 plot, i1av-vxs, alog10(i1aint), LINESTYLE=0, xstyle=1, ystyle=1, $ xtitle='speed in wind frame /(km/s)', ytitle='I1a VDF (-- I1b)' oplot, i1bv-vxs, alog10(i1bint), LINESTYLE=2 oplot, i1av-vxs, zero, PSYM=7 charge=fltarr(32) for ic= 0, 20 do begin charge(ic)= alog10(1.0) endfor for ic= 21,31 do begin charge(ic)= alog10(2.0) endfor plot, (i1av+i1bv)/2.-vxs, alog10(i1bint)-alog10(i1aint), $ LINESTYLE=0, xstyle=1, ystyle=1, xtitle='thermal speed /(km/s)', ytitle='Log(I1b/I1a)' oplot, (va_noa+vb_noa)/2.-vxs, alog10(fbint_noa)-alog10(faint_noa), LINESTYLE=2 oplot, (i1av+i1bv)/2.-vxs, charge, PSYM=2 plot, va_noa-vxs, alog10(faint_noa), LINESTYLE=0, xstyle=1, ystyle=1, $ xtitle='thermal speed /(km/s)', ytitle='I1a VDF (-- I1b)' oplot, vb_noa-vxs, alog10(fbint_noa), LINESTYLE=2 ;******************************************************************************** ; Contour plot of combined proton data [-300,300] including proton core and beam ;******************************************************************************** !p.multi=[0,1,1] contour,[[0,0],[1,1]],/nodata,xstyle=4,ystyle=4 px=!x.window*!d.x_vsize py=px contour, fpz,wx_i,wz_i,/noerase,/xst,/yst,pos=[px(0),py(0),px(1),py(1)],/dev,$ levels=[0.01,0.032,0.1,0.2,0.4,0.6,0.8],$ C_LINESTYLE=[1,1,0,0,0,0,0],$ C_ANNOTATION=['0.01','0.032','0.1','0.2','0.4','0.6','0.8'],$ xrange=[-300,300], yrange=[-300,300], $ xtitle='thermal speed /(km/s)', ytitle='thermal speed /(km/s)' plots,wbx,wbz for j=0, jdimen-1 do begin for k=0, knoal-1 do begin plots, wx_1(2,j,k), wz_1(2,j,k), psym=1 endfor endfor device,/close endif ;device,/close, filename = strplot, YSIZE = 10.0, XSIZE = 10.0 ;set_plot,'win' continue_spectrum: idimen = 12 jdimen = 12 kdimen = 32 ;print, idimen, jdimen, kdimen close,1 ;close,2 ;close,3 endwhile close,15 endfor ;---------- end loop over as many data sets as listed in names.txt -------------- close,4 close,5 close,6 end ;-------------------------- end programme ----------------------------------------