; file: ssb1_instruct.pro = IDL pro SSB1 "Gene Avrett" as instruction ; last: Oct 19 2018 Rob Rutten WeiHai ; site: https://robrutten.nl/rrweb/rjr-edu/exercises/ssb ; note: should be same as in instruction ; there is a STOP after every item; type .c to continue to next ; =========================================== START OF MAIN PROGRAM ; ========================= read FALC model close,1 openr,1,'falc.dat' nh=80 dummy='' for iskip=1,4 do readf,1,dummy falc=fltarr(11,nh) ; 11 FALC columns readf,1,falc h=reform(falc[0,*]) tau5=reform(falc[1,*]) colm=reform(falc[2,*]) temp=reform(falc[3,*]) vturb=reform(falc[4,*]) nhyd=reform(falc[5,*]) nprot=reform(falc[6,*]) nel=reform(falc[7,*]) ptot=reform(falc[8,*]) pgasptot=reform(falc[9,*]) dens=reform(falc[10,*]) ; much easier with astrolib routine readcol,'falc.dat',h,tau5,colm,temp,vturb,nhyd,nprot,nel,ptot,$ pgasplot,dens,skipline=4 ; ============================ FALC plots ; temperature stratification plot,h,temp,yrange=[3000,10000],$ xtitle='height [km]',ytitle='temperature [K]',$ title=' FALC temperature stratification' STOP ; ptot - colm plot,colm,ptot,$ xtitle='column mass [g cm-2]',ytitle='pressure [dyne cm-2]',$ title=' pressure and mass' plot,colm,ptot,/xlog,/ylog,$ xtitle='column mass [g cm-2]',ytitle='pressure [dyne cm-2]',$ title=' pressure and mass' STOP ; measure g print,' solar surface gravity g = ',total(ptot/colm)/nh STOP ; element mix mH=1.67352D-24 ; mass of hydrogen atom (AQ) mHe=3.97*mH ; mass helium atom plot,h,nhyd*mH/dens,/ynozero,$ xtitle='height [km]',ytitle='H mass density/ gas density',$ title=' hydrogen fraction' plot,h,(nhyd*mH + 0.1*nhyd*mHe)/dens,/ynozero,$ xtitle='height [km]',ytitle='H+He mass density/ gas density',$ title=' hydrogen + helium fraction' plot,h,1 - (nhyd*mH + 0.1*nhyd*mHe)/dens,/ynozero,$ xtitle='height [km]',ytitle='1 - (H+He mass)/gas',$ title=' remaining element fraction' print,' fraction mass density metals = ',$ total(1-(nhyd*mH + 0.1*nhyd*mHe)/dens)/nh STOP ; colm - h plot,h,colm,/ylog,$ xtitle='height [km]',ytitle='column mass [g cm-2]',$ title=' mass and height' STOP ; density scale height plot,h,dens,/ylog,$ xtitle='height [km]',ytitle='density [g cm-3]',$ title=' scale height' fit=fltarr(nh) Hdens=120. ; trial and error fit for i=0,nh-1 do fit[i] = 1.E-6*exp(-h[i]/Hdens) oplot,h,fit,linestyle=2 xyouts,1000,1E-9,'H = 120 km' STOP ; ideal gas law k=1.38065D-16 ; Boltzmann constant pgas=ptot*pgasptot ; gas pressure = total without turbulent plot,h,pgas,/ylog,$ xtitle='height [km]',ytitle='gas pressure, (N_H+N_e) KT',$ title=' pressure and height' oplot,h,(nhyd+nel)*temp*k,linestyle=1 plot,h,(nhyd+nel)*temp*k/pgas,/ynozero,$ xtitle='height [km]',ytitle='(N_H+N_e)kT/P_g',$ title=' ideal das test' plot,h,(nhyd*1.1+nel)*temp*k/pgas-1.,$ xtitle='height [km]',ytitle='(N_H+N_He+N_e)kT/P_g - 1',$ title=' ideal gas test' STOP ; ionization balance plot,h,nhyd,/ylog,$ xtitle='height [km]',ytitle='particle density [cm-3]',$ title=' hydrogen ionization' oplot,h,nprot,linestyle=1 oplot,h,nel,linestyle=2 oplot,h,nel-nprot,linestyle=1 xyouts,h[53],2*nhyd[53],'N!dH!n' xyouts,h[53],2*nel[53],'N!de!n' xyouts,h[53],0.2*nprot[53],'N!dp!n' xyouts,h[40],0.2*(nel[40]-nprot[40]),'N!de!n-N!dp!n' STOP ; hydrogen ionization fraction plot,h,nprot/nhyd,/ylog,$ xtitle='height [km]',ytitle='n_p / N_H',$ title=' hydrogen ionization fraction' STOP ; photon density print,' deep photon density = ', 20.*8D3^3 print,' high photon density = ', 20.*5.7D3^3/3 ; halve hemel + randverzwakking ; =========================================== EARTH plots STOP ; read earth.dat close,2 openr,2,'earth.dat' nhE=25 earth=fltarr(5,nhE) readf,2,earth hE=earth[0,*] pgasE=10.^earth[1,*] tempE=earth[2,*] densE=10.^earth[3,*] npartE=10.^earth[4,*] STOP ; plot Earth temperature stratification y=tempE plot,hE,tempE,$ xtitle='height [km]',ytitle='temperature [K]',$ title=' earth temperature stratification' xyouts,min(hE)+0.1*(max(hE)-min(hE)),$ min(y)+0.8*(max(y)-min(y)),$ 'EARTH ATMOSPHERE' STOP ; plot Earth density and pressure y=densE plot,hE,y,/ylog,$ xtitle='height [km]',ytitle='gas density [g cm-3]',$ title=' density and height' xyouts,min(hE)+0.5*(max(hE)-min(hE)),$ min(y)*10^(0.8*(alog10(max(y))-alog10(min(y)))),$ 'EARTH ATMOSPHERE',alignment=0.5 y=npartE plot,hE,y,/ylog,$ xtitle='height [km]',ytitle='particle density [cm-3]',$ title=' density and height' xyouts,min(hE)+0.5*(max(hE)-min(hE)),$ min(y)*10^(0.8*(alog10(max(y))-alog10(min(y)))),$ 'EARTH ATMOSPHERE',alignment=0.5 y=pgasE plot,hE,y,/ylog,$ xtitle='height [km]',ytitle='pressure [dyne cm-2]',$ title=' pressure and height' xyouts,min(hE)+0.5*(max(hE)-min(hE)),$ min(y)*10^(0.9*(alog10(max(y))-alog10(min(y)))),$ 'EARTH ATMOSPHERE',alignment=0.5 STOP ; Earth pressure and densities in one graph y=densE/densE(0) plot,hE,y,/ylog,$ xtitle='height [km]',ytitle='densities, pressure [normalized]',$ title=' normalized stratifications' oplot,hE,npartE/npartE(0),linestyle=2 oplot,hE,pgasE/pgasE(0),linestyle=3 xyouts,min(hE)+0.5*(max(hE)-min(hE)),$ min(y)*10^(0.8*(alog10(max(y))-alog10(min(y)))),$ 'EARTH ATMOSPHERE',alignment=0.5 STOP ; mean molecular weight y=densE/npartE/mH plot,hE,y,/ynozero,$ xtitle='height [km]',ytitle='mean molecular weight',$ title=' particle mass' xyouts,min(hE)+0.5*(max(hE)-min(hE)),$ min(y)+0.1*(max(y)-min(y)),$ 'EARTH ATMOSPHERE',alignment=0.5 ; plot extends below min(y) so this range fraction isn't exact STOP ; fit Earth density scale height y=densE plot,hE,y,/ylog,$ xtitle='height [km]',ytitle='gas density [g cm-3]',$ title=' scale height' xyouts,min(hE)+0.5*(max(hE)-min(hE)),$ min(y)*10^(0.9*(alog10(max(y))-alog10(min(y)))),$ 'EARTH ATMOSPHERE',alignment=0.5 fit=fltarr(nhE) HEdens=7. ; trial and error fit (well, listed in AQ too) for i=0,nhE-1 do fit[i] = densE[0]*exp(-hE[i]/HEdens) oplot,hE,fit,linestyle=2 xyouts,120,1E-10,'H = 7 km' STOP ; earth and solar particle density ratio print, ' sun-earth particle density ratio =',$ (nhyd[69]*1.1+nel[69])/npartE[0] end