Program 1. Use linear interpolation to separate each layer of the atmosphere into 200 layers.
while not eof(1) do begin
readf,1,format='(18x, f6.1, 3x, f5.0, 3x, f5.1, 4x, f4.0)', phpao, Ho, Tco, RHo
for kk=0, nn-2 ,0 do begin
Program 2. Use the result from program 1 to calibrate the refraction index of each layer.
pro refraction_index_20111028
; This program is used to calculate the index of refraction n from the observed Temperature, Pressure, and Height Data. reference is from ''Refractive index of air: new equations for the visible and near infrared, Phlip E. Ciddor, 1996'
close,/all
aaa=' '
phpa=double(0.0d0);observed pressure in hPa P=double(0.0d0) ;observed pressure in pa H=double(0.0d0) ;Height in meters
Tc=double(0.0d0) ;Temperature in centigrade Tk=double(0.0d0) ;Temperature in Kelvin RH=double(0.0d0) ;relative humidity in %
rho=double(0.0d0)
; density of atmosphere rhoaxs=double(0.0d0)
; The density of dry air under standared condition. (Tc=15, P=101325, 0% RH, 450ppm of CO2)
rhoa=double(0.0d0)
; The density of dry air (0% RH) under the observed conditions.
rhows=double(0.0d0)
; The density of pure water vapor under standared condition. (Tc=15, P=101325, 100%
RH, 450ppm of CO2) rhow=double(0.0d0)
; The density of water vapor under observed condition.
nprop=double(0.0d0)
; The index of refraction
nas=double(0.0d0) ; naxs=double(0.0d0)
; The index of refraction of dry air for a particular wavelength under standard condition nws=double(0.0d0)
; The index of refraction of pure water vapor for a particular wavelength under standared conditions.
f=double(0.0d0) ;f is the enhancement factor of water vapor in air svp=double(0.0d0) ;
Xw=double(0.0d0) ; Xw is the molar fraction of water vapor in moist air Z=double(0.0d0) ; Z is the compressibility of the moist air
;The density of dry air under standared condition. (Tc=15, P=101325, 0% RH, 450ppm of CO2) faxs=double(0.0d0) ;f is the enhancement factor of water vapor in air
svpaxs=double(0.0d0) ;
Xwaxs=double(0.0d0) ; Xw is the molar fraction of water vapor in moist air Zaxs=double(0.0d0) ; Z is the compressibility of the moist air
;The density of dry air (0% RH) under the observed conditions.
fa=double(0.0d0) ;f is the enhancement factor of water vapor in air svpa=double(0.0d0) ;
Xwa=double(0.0d0) ; Xw is the molar fraction of water vapor in moist air Za=double(0.0d0) ; Z is the compressibility of the moist air
; The density of pure water vapor under standared condition. (Tc=15, P=101325, 100% RH,
; The density of water vapor under observed condition.
fw=double(0.0d0) ;f is the enhancement factor of water vapor in air svpw=double(0.0d0) ;
Xww=double(0.0d0) ; Xw is the molar fraction of water vapor in moist air Zw=double(0.0d0) ; Z is the compressibility of the moist air
Tcs=double(15.0d0) ; Temperature at standared condition (Tc=15) in Centigrade Tks=double(Tcs+(273.15d0)) ; Temperature at standared condition in Kelvin Ps=double(101325d0) ; Pressure in Standared condition
RH0=double(0.0d0) ; Humidity of Dry air RH100=double(100.0d0); Humidity of 100% air
Ma=double((1d-3)*(28.9635d0+((12.011d-6)*(450.0d0-400.0d0))))
; The molar mass of dry air containing 450 ppm of CO2 Mw=double(0.018015d0) ; The molar mass of water vapor while not eof(1) do begin
readf,1,format='(f9.4, 2x, f10.4, 2x, f8.4, 2x, f8.4)', phpa, H, Tc, RH rho=double(((P*Ma)/(Z*R*Tk))*(1.0d0-(Xw*(1.0d0-(Mw/Ma)))))
rhoaxs=double(((Ps*Ma)/(Zaxs*R*Tks))*(1.0d0-(Xwaxs*(1.0d0-(Mw/Ma))))) rhoa=double(((P*Ma)/(Z*R*Tk))*(1.0d0-(Xw*(1.0d0-(Mw/Ma)))))
rhows=double(((Ps*Ma)/(Zws*R*Tks))*(1.0d0-(Xwws*(1.0d0-(Mw/Ma))))) rhow=double(((P*Ma)/(Z*R*Tk))*(1.0d0-(Xw*(1.0d0-(Mw/Ma)))))
nas=double((((k1/(k0-((sigma)^2.0d0)))+(k3/(k2-((sigma)^2.0d0))))/1d8)+1.0d0)
printf,2,format='(f9.4, 2x , f10.4, 2x, f8.4, 2x ,f8.4, 2x, f8.4, 2x, d15.5)', Phpa, H, Tc, Tk, RH, nprop*1.0d8
endwhile
close,/all print,'ok' end
Program 3. Calibrate the altitude angle of the real sun.
pro dii_20111028_4817
; This program is to calbrate the altitude angle (from the real horizon) of the real sun.
close,/all
openw,1,'G:\WORK\All_Observed_DatasX\20111028\dii_4817.txt'
;Rs=(1913.807d0/2.720353d0)/2.0d0 ;the solar angular diameter from JPL Horizon Re=6374248.3d0; Earth's r in meters
h=5.0d0 ; Height of the observe location in meters
diptheta=0.0d0 ; When h>o, the angle between "Horizon" to "real horizon"
phy=25.1851093d0 ; The longituge of the observation location.
; ===== The sun's RA and Dec during transit =====
RAt=210.140491d0 ; Transit RA in degree Dect=-12.976758d0 ; Transit Dec in degree
; ==== The sun's RA and Dec during take photo ====
RA=210.155095d0 ; Photoed RA in degree Dec=-13.053414d0 ; Photoed Dec in degree
;Dec=23+(8.12/60)
; ==== The time in hours ====
Tt=3.637015d0 ; Sun Transit time in hours T=9.143056d0 ; Take photo time
dT=0.0d0 ; Delta T dRA=0.0d0 ;Delta RA dDec=0.0d0 ; Delta Dec
dDegree=0.0d0 ; From Transit time to Photoed time, the angle to the zenith dT=T-Tt
RAp=RA-RAt dDec=Dec-Dect
diptheta=double((acos(Re/(Re+h)))/!dtor) ; dip angle of height in degree
L=double(asin(((sin((90.0d0)*!dtor))*(sin((Dec)*!dtor)))/(sin((90.0d0-phy)*!dtor))))
thtt=0.0d0 aaaaa=0.0d0 ThetaZ=0.0d0
alpha=double((asin((sin(phy*!dtor)) / (cos(dec*!dtor))))/!dtor) L=double((asin((sin(dec*!dtor)) / (cos(phy*!dtor))))/!dtor) thetaa=double((asin((sin(L*!dtor))*(sin(alpha*!dtor))))/!dtor) aaaaa=double(angle-thetaa)
thetaz=((asin(((sin(alpha*!dtor))*(cos(L*!dtor)))/(sin(phy*!dtor))))/!dtor) Beta=(acos((sin((thetaz+angle)*!dtor))*(cos(phy*!dtor))))/!dtor
Ll=double((asin((sin(aaaaa*!dtor)) / (sin(Beta*!dtor))))/!dtor)
Lx= (acos((1/(tan((thetaz+angle)*!dtor))) * (1/(tan((beta)*!dtor)))))/!dtor N=double(Lx+Dec-90.0d0)
ThetaX=double((asin((sin(N*!dtor))*(sin(Beta*!dtor))))/!dtor)
ThetaAa=double((asin(((cos(Beta*!dtor))*(sin(N*!dtor)))/(cos(ThetaX*!dtor))))/!dtor) ThetaAaa=double(360.0d0-(90.0d0+Ll-ThetaAa))
printf,1,format='(f11.7, 2x, f11.7)',ThetaX, ThetaAaa
print, 'ok' close,/all end
Program 4. Calibrate the modeled refracted sun
pro n2n_20111028_4817
; This program is to calibrate the modeled refracted sun close,/all
Re=6374248.3d0; Earth's R in meters without height
rsun=((1931.493d0/3600.d0)/2.0d0) ; the diameter of the sun, in degree h=5.0d0 ; the height of the observer
Reh=0.0d0 Reh=Re+h aaa=' '
dii=0.0d0 ; the altitude (from the real horizon) of the real sun
openr, 6,'G:\WORK\All_Observed_DatasX\20111028\dii_4817.txt'
deltain=dblarr(na-1.0d0) while not eof(1) do begin
readf,1,format='(11x , f10.4, 31x, d16.6)', Zo1, No1
for i=0.0d0, nnn do begin tt=(i*2.0d0*!dpi/nnn)/!dtor
ys(i)=(asin((sin((tt)*!dtor))*(sin((rsun)*!dtor))))/!dtor xs(i)=(asin((sin((90.0d0-tt)*!dtor))*(sin((rsun)*!dtor))))/!dtor dip=dii+ys(i)
bb=90.0d0+dip
alphaa=(asin(((Reh)*(sin((bb)*!dtor)))/(Re+Ri)))/!dtor acc=180.0d0-alphaa-bb
bx=(asin(((Re+Ri)*(sin((bb)*!dtor)))/(des)))/!dtor ax=180.0d0-bb-bx
printf,3,format='(f11.8, 2x, f11.8, 2x, f10.6, 2x, f10.6)',alphaa, acc, Xs(i), Ys(i) endfor
rxx=0.0d0
uu=uu+1
if uu eq 10000 then goto, jump2 if abs(aXX) le err then goto, jump2 goto, jump1
jump2:
printf,2,format='(f14.7, 2x, f10.2, 2x, f10.7, 2x, f10.7, 2x, f13.10)', XX, YY, dip, deltain(na-2), rxx
endwhile
;--- ---
print, 'ok' close,/all end
Program 5. Plot the real sun, modeled sun on the image of the observed sun.
pro RealSunLocation_20111028
; This program is used to plot the real sun, modeled sun on the image of the observed sun.
close,/all
file = 'G:\WORK\All_Observed_DatasX\20111028\images\_MG_4817_4836_4853X.jpg' image=read_image(file)
scale=2.232932d0; Image scale arcsec/pixel Rs=(1931.493d0/scale)/2.0d0
Re=6374248.3d0; Earth's r in meters
h=5.0d0 ; Height of the observe location in meters
diptheta=0.0d0 ; When h>o, the angle between "Horizon" to "real horizon"
phy=25.1851093d0 ; the longitude of the observation location RS=2652.0d0 ; The apparent sea level in the image (in pixel) RH=2567.0d0; The real horizon in the image (in pixel)
anglex=0.0d0 ; the altitude (fram the real horizon) of the real sun openr,4,'G:\WORK\All_Observed_DatasX\20111028\dii_4817.txt' readf,4,format='(f11.7, 2x, f11.7)', anglex, diptheta
close,4
;---draw the apparent sea level--- xfit2=fltarr(2)
;---draw the real sea level--- xfit3=fltarr(2)
yfit3=fltarr(2) xfit3(0)=0 xfit3(1)=3615
yfit3(0)=RH+(diptheta*3600/scale) yfit3(1)=RH+(diptheta*3600/scale)
;---draw the scale---
plot, xfit, yfit, linestyle=0, xrange=[0,3480], yrange=[5200,0], /isotropic , xtitle='pixel', ytitle='pixel', thick=2, xstyle=1, ystyle=1
TV,image , 2225 , 1415 , xsize=7075 , /true
oplot, xfit, yfit,color=1, linestyle=2 ; draw the real horizon oplot, xfit2, yfit2,color=3 ; draw the apparent sea level oplot, xfit3, yfit3,color=3,linestyle=1 ; draw the real sea level oplot, xfit4, yfit4,color=3 ;draw the scale
xyouts, 50, 2960,'Real sea level at altitude=5 m',color=3,CHARTHICK=2,CHARSIZE =0.5 xyouts, 50, 2800,'Real horizon',color=1,CHARTHICK=2,CHARSIZE =0.5
xyouts, 50, 2880,'Apparent sea level at altitude=5 m',color=3,CHARTHICK=2,CHARSIZE =0.5 xyouts, 80, 50+5*60/scale,'5 arcmin',color=3,CHARTHICK=2,CHARSIZE =0.5
;---- The first sun ----
while not eof(1) do begin
readf,1,format='(27x, f9.6, 3x, f9.6)', XXs, YYs xs(aa)=(XXs*3600.0d0/scale)+x
ys(aa)=-(YYs*3600.0d0/scale)+y aa=aa+1
oplot, xs, ys, color=1, thick=1 , psym=1, symsize=0.05, linestyle=0 xxc=fltarr(1)
XYouts, X-1200, 900, 'Unrefracted sun',color=1,CHARTHICK=2,CHARSIZE =0.5 XYouts, X-1200, 980, '09:08:35 UT',color=1,CHARTHICK=2,CHARSIZE =0.5 XYouts, X-1200, 300, 'Refracted sun',color=5,CHARTHICK=2,CHARSIZE =0.5 XYouts, X-1200, 380, '09:08:35 UT',color=5,CHARTHICK=2,CHARSIZE =0.5
Xxxx=X
ysr=fltarr(513) ac=0.0d0
while not eof(3) do begin
readf,3,format='(27x, f9.6)', XXsr xsr(ac)=(XXsr*3600.0d0/scale)+Xxxx ac=ac+1
endwhile aaa=' ' readf,2,aaa ab=0.0d0
while not eof(2) do begin
readf,2,format='(40x, f10.7)', YYsr ysr(ab)=-(YYsr*3600.0d0/scale)+RH ab=ab+1
endwhile
oplot, xsr,ysr, color=5, thick=1 , psym=1 , symsize=0.05 ,linestyle=0 close,/all
while not eof(1) do begin
readf,1,format='(27x, f9.6, 3x, f9.6)', XXs, YYs xs(aa)=(XXs*3600.0d0/scale)+x
ys(aa)=-(YYs*3600.0d0/scale)+y aa=aa+1
oplot, xs, ys, color=1, thick=1 , psym=1, symsize=0.05, linestyle=0 xxc=fltarr(1)
XYouts, X-1200, 1800, 'Unrefracted sun',color=1,CHARTHICK=2,CHARSIZE =0.5 XYouts, X-1200, 1880, '09:10:59 UT',color=1,CHARTHICK=2,CHARSIZE =0.5
XYouts, X-1200, 1200, 'Refracted sun',color=5,CHARTHICK=2,CHARSIZE =0.5 XYouts, X-1200, 1280, '09:10:59 UT',color=5,CHARTHICK=2,CHARSIZE =0.5
Xxxx=X
while not eof(3) do begin
readf,3,format='(27x, f9.6)', XXsr xsr(ac)=(XXsr*3600.0d0/scale)+Xxxx ac=ac+1
endwhile aaa=' ' readf,2,aaa ab=0.0d0
while not eof(2) do begin
readf,2,format='(40x, f10.7)', YYsr ysr(ab)=-(YYsr*3600.0d0/scale)+RH ab=ab+1
endwhile
oplot, xsr,ysr, color=5, thick=1 , psym=1 , symsize=0.05 ,linestyle=0 close,/all
while not eof(1) do begin
readf,1,format='(27x, f9.6, 3x, f9.6)', XXs, YYs xs(aa)=(XXs*3600.0d0/scale)+x
ys(aa)=-(YYs*3600.0d0/scale)+y aa=aa+1
oplot, xs, ys, color=1, thick=1 , psym=1, symsize=0.05, linestyle=0 xxc=fltarr(1)
yyc=fltarr(1) xxc(0)=X
yyc(0)=Y
XYouts, X-1200, 2450, 'Unrefracted sun',color=1,CHARTHICK=2,CHARSIZE =0.5 XYouts, X-1200, 2530, '09:13:17 UT',color=1,CHARTHICK=2,CHARSIZE =0.5 XYouts, X-1200, 2000, 'Refracted sun',color=5,CHARTHICK=2,CHARSIZE =0.5 XYouts, X-1200, 2080, '09:13:17 UT',color=5,CHARTHICK=2,CHARSIZE =0.5
Xxxx=X
while not eof(3) do begin
readf,3,format='(27x, f9.6)', XXsr xsr(ac)=(XXsr*3600.0d0/scale)+Xxxx ac=ac+1
endwhile aaa=' ' readf,2,aaa ab=0.0d0
while not eof(2) do begin
readf,2,format='(40x, f10.7)', YYsr ysr(ab)=-(YYsr*3600.0d0/scale)+RH ab=ab+1
endwhile
oplot, xsr,ysr, color=5, thick=1 , psym=1 , symsize=0.05 ,linestyle=0
device,/close close,/all print, 'ok' end
Appendix