pro Ratio ;--- 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 data ;+++++++++++++++++++++++++++++++++++ ;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 +++++++++++++++++++++++++++++++++++++ ;openr,1, 'h1y75d010h02m38s34_hdm.1' ;openr,1, 'h2y76d023h11m54s56_hdm.1' ;openr,1, 'h2y76d107h07m50s54_hdm.1' ;openr,1, 'h2y76d086h14m42s38_hdm.0' ;openr,1, 'h2y76d086h14m44s00_hdm.1' ;openr, 1, 'h2y76d107h07m44s44_hdm.0' openr,1, 'h2y76d066h00m01s45_hdm.0' ;openr,1, 'h2y76d077h13m32s13_hdm.0' ;openr,1, 'h2y76d067h22m32s33_hdm.1' ;;;openr,1, 'h2y76d072h10m04s38_hdm.0' ;openr,1, 'h2y76d107h07m43s23_hdm.0' ;openr,1,'h2y76d107h07m40s41_hdm.0' ;openr,1,'h2y76d107h07m42s02_hdm.0' ;openr,15, 'names.txt' ;openr,15, 'C:\Dokumente und Einstellungen\Marsch\Eigene Dateien\Helios\Filenames\H1-75-010.txt' ;while (not eof(15)) do begin ; fnametemp='' ;' Y76D107/h2y76d107h00m20s31_hdm.0' ; readf,15, fnametemp ; print, fnametemp ; filen=strmid(fnametemp,1,33) ; fname=strmid(fnametemp,9,18) ;print, fname ;strdata=string ;strmom =string ;strplot=string strdata='Ion_vdf_data.dat' print,strdata strmom='Ion_moments.dat' print,strmom strplot='Ion_plots.ps' print,strplot ;++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ; Open output file ;openr,1,'C:\Dokumente und Einstellungen\marsch\Eigene Dateien\Helios\'+filen 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) ;print, countp 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 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 kstart=maxk print, 'kstart ', kstart ;------------------------------------------------------------------------------------ knoal = kdimen khignoa = klow + kstart ratcrit = alog10(2) print, 'ratcrit', ratcrit kbegin = klow + kstart 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 k_min = knoal print,'proton upper and alpha lower limits found' print,'knoal= ', knoal,' k_min=', k_min ;-----------out put information and data for interpolation ----------------- openw, 10,'proton_vdf_3d.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 ------------------------------- print, 'knoal', knoal, 'kdimen', kdimen 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 print, 'k1 ',k1, 'kdimen ', kdimen 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. ;--------------------------------------------------------------------------- openw, 11,'alpha_vdf_3d.dat' ;information and data for interpolation k_a= kdimen - k_min print, 'k_a ', k_a if (k_a LE 1) then begin k_a = k_a + 1 k_min = k_min-1 endif 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 ;continue_alpha: 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) ;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) = v_a(km_a)*COS_A*COS_E vy_max(1) = v_a(km_a)*SIN_A*COS_E vz_max(1) = v_a(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 if (m_spc eq 1) then begin 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 k_spc= k_a 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) 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 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 close,3 close,2 ;--------------------- 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 ;**************************** 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 ;----------------- 2-d plot with alphas included ------------------ ; x radial ; z northward ;------------------------------------------------------------------ ;device, filename = strplot, XSIZE = 10.0 ;!x.range=[-200,800] ;plot,alog10(fp(imax+1,jmax-1,*)),XTICKS=kend,color=2 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)' ;nlevels=50 ;,/ISOTROPIC 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 ;plots,wx_1(2,0,*),wz_1(2,0,*),psym=1 ;plots,wx_1(2,1,*),wz_1(2,1,*),psym=1 ;plots,wx_1(2,2,*),wz_1(2,2,*),psym=1 ;plots,wx_1(2,3,*),wz_1(2,3,*),psym=1 ;plots,wx_1(2,4,*),wz_1(2,4,*),psym=1 ;plots,wx_1(2,5,*),wz_1(2,5,*),psym=1 ;plots,wx_1(2,6,*),wz_1(2,6,*),psym=1 endif ;******************************************************************************* ; speed range [-200,800] for protons and alpha particles ; ;******************************************************************************* ;device,/color,YSIZE=10.0,XSIZE=10.0,filename='ions_1D_data.ps' ;device, filename = strplot, XSIZE = 10.0 ;!x.range=[0,kdimen-1] !p.multi=[0,1,2] maxfp=max(fp) print,'maxfp: ', maxfp print, imax, jmax, imax+1, jmax+1, imax-1, jmax-1 ;print, fp 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 ;plots,wx_1(2,0,*),wz_1(2,0,*),psym=1 ;plots,wx_1(2,1,*),wz_1(2,1,*),psym=1 ;plots,wx_1(2,2,*),wz_1(2,2,*),psym=1 ;plots,wx_1(2,3,*),wz_1(2,3,*),psym=1 ;plots,wx_1(2,4,*),wz_1(2,4,*),psym=1 ;plots,wx_1(2,5,*),wz_1(2,5,*),psym=1 ;plots,wx_1(2,6,*),wz_1(2,6,*),psym=1 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 ;for ic= 29,31 do begin charge(ic)= alog10(6.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 ;******************************************************************************** ;device, filename = strplot, XSIZE = 10.0 !p.multi=[0,1,1] ;!x.range=[-300,300] ;!y.range=[-300,300] 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 ;plots,wx_1(2,0,*),wz_1(2,0,*),psym=1 ;plots,wx_1(2,1,*),wz_1(2,1,*),psym=1 ;plots,wx_1(2,2,*),wz_1(2,2,*),psym=1 ;plots,wx_1(2,3,*),wz_1(2,3,*),psym=1 ;plots,wx_1(2,4,*),wz_1(2,4,*),psym=1 ;plots,wx_1(2,5,*),wz_1(2,5,*),psym=1 ;plots,wx_1(2,6,*),wz_1(2,6,*),psym=1 device,/close ;device,/close, filename = strplot, YSIZE = 10.0, XSIZE = 10.0 set_plot,'win' ;endwhile close,15 ;---------- end loop over as many data sets as listed in names.txt -------------- end ;-------------------------- end programme ----------------------------------------