I spend three hours on this code yesterday
Here is the wrong one
pro plot_solar_wind
openr, file, 'files/solar_wind.dat', /get_lun
header = ''
for i = 0, 60 do readf, file, header
month_hours = fltarr(31 * 24)
Vsw = fltarr(31 * 24)
Nsw = fltarr(31 * 24)
Psw = fltarr(31 * 24)
Rmp = fltarr(31 * 24)
day = 0
hour = 0
Vsw_in = 0.0
Nsq_in = 0.0
pts = 0
while not eof(file) do begin
readf, file, day, hour, Nsw_in, Vsw_in, format = '(i2, 9x, i2, 10x, d14.5, d14.3)'
month_hours[pts] = (day - 1) * 24 + hour
Vsw[pts] = Vsw_in
Nsw[pts] = Nsw_in
; print, Vsw_in
pts += 1
endwhile
free_lun, file
openw, file, 'files/some_kind_of_data.dat', /get_lun
for i = 0, pts - 1 do begin
if Vsw le 0 then Vsw = !Values.F_NaN
if Nsw le 0 then Nsw = !Values.F_NaN
if Vsw gt 0 and Nsw gt 0 then begin
Vsw *= 1E+3
Nsw *= 1E-6
; print, Vsw, Nsw
Psw[i] = (1.67262E-27 * Nsw[i]) * (Vsw[i] ^ 2)
Rmp[i] = ((31000E-9)^2 / (1.2566E-6 * Psw[i]))^(1 / 6.0)
print, Vsw[i], Nsw[i], Psw[i], Rmp[i]
if Rmp[i] gt 6.6 then begin
printf, file, i, Vsw[i], Rmp[i], ' Magnetosphere'
endif else begin
printf, file, i, Vsw[i], Rmp[i], ' Solar wind'
endelse
endif else begin
Psw[i] = !Values.F_NaN
Rmp[i] = !Values.F_NaN
printf, file, i, Vsw[i], Rmp[i], ' No data'
endelse
endfor
free_lun, file
window, 1
plot, Vsw, ytitle = 'Speed (m/s)', min = 0
window, 2
plot, Nsw, ytitle = 'Density (#/m^3)', min = 0
window, 3
plot, Psw, ytitle = 'Pressure (#/m^3)', min = 0
end
And here is the right one
pro plot_solar_wind
openr, file, 'files/solar_wind.dat', /get_lun
header = ''
for i = 0, 60 do readf, file, header
month_hours = fltarr(31 * 24)
Vsw = fltarr(31 * 24)
Nsw = fltarr(31 * 24)
Psw = fltarr(31 * 24)
Rmp = fltarr(31 * 24)
day = 0
hour = 0
Vsw_in = 0.0
Nsq_in = 0.0
pts = 0
while not eof(file) do begin
readf, file, day, hour, Nsw_in, Vsw_in, format = '(i2, 9x, i2, 10x, d14.5, d14.3)'
month_hours[pts] = (day - 1) * 24 + hour
Vsw[pts] = Vsw_in
Nsw[pts] = Nsw_in
; print, Vsw_in
pts += 1
endwhile
free_lun, file
openw, file, 'files/some_kind_of_data.dat', /get_lun
for i = 0, pts - 1 do begin
if Vsw[i] le 0 then Vsw[i] = !Values.F_NaN
if Nsw[i] le 0 then Nsw[i] = !Values.F_NaN
if Vsw[i] gt 0 and Nsw[i] gt 0 then begin
Vsw[i] *= 1E+3
Nsw[i] *= 1E+6
; print, Vsw[i], Nsw[i]
Psw[i] = (1.67262E-27 * Nsw[i]) * (Vsw[i] ^ 2)
Rmp[i] = ((31000E-9)^2 / (1.2566E-6 * Psw[i]))^(1 / 6.0)
print, Vsw[i], Nsw[i], Psw[i], Rmp[i]
if Rmp[i] gt 6.6 then begin
printf, file, i, Vsw[i], Rmp[i], ' Magnetosphere'
endif else begin
printf, file, i, Vsw[i], Rmp[i], ' Solar wind'
endelse
endif else begin
Psw[i] = !Values.F_NaN
Rmp[i] = !Values.F_NaN
printf, file, i, Vsw[i], Rmp[i], ' No data'
endelse
endfor
free_lun, file
window, 1
plot, Vsw, ytitle = 'Speed (m/s)', min = 0
window, 2
plot, Nsw, ytitle = 'Density (#/m^3)', min = 0
window, 3
plot, Psw, ytitle = 'Pressure (#/m^3)', min = 0
end
Who see the difference without cheating