;ex2.pro: ; ; Program: This program shows an example of the regression routine ; ZREGR.PRO. The output is 6 pages of postcript plots. Solar ; flux measurements at 10.7 cm wavelength integrated over ; the disk of the sun have been taken in Canada since 1947, ; while sun spot numbers have been measured from Belgium and ; elsewhere since even before the 1700's. These two distinctly ; different measurements have a strikingly similar temporal ; behavior for time scales of a month or longer. This program ; takes advantage of this similar temporal dependence to ; create a modeled solar flux time series dating back to the ; 1700's using the regression routine ZREGR.PRO. The solar ; flux time series is modeled in the regression to have a ; cubic polynomial relationship with the sun spot number time ; series. ; To run: ; 1) Save this procedure as ex2.pro, and in same directory ; save the routine ZREGR.PRO. ; 2) In the same directory save the following ASCII data ; files: ; ; monthssn.dat <== Monthly mean sun spot number ; (SSN) data from the Sunspot ; Index Data Center (SSDC) at ; the Royal Observatory of ; Belgium. Time coverage: Jan ; 1749 through December 2000. ; yearssn.dat <== Annual average SSN data from ; SSDC. Time coverage: Year 1700 ; through year 2000. ; solar47-00.mon <== Monthly mean 10.7 cm solar ; flux data from the National ; Geophysical Data Center (NGDC) ; at National Oceanic and ; Atmospheric Administration ; (NOAA). Time coverage: Jan ; 1947 through December 2000. ; 3) Get into IDL and type ; ; .r ex2 ; ; To see plots type: ; ghostview idl.ps (or use ghostscript, etc.) month_ssn=fltarr(3024) ; 3024: Jan 1749, Feb 1749, ..., Dec 2000 year_ssn=fltarr(301) ; 301: 1700, 1701, ..., 2000 month_f107=fltarr(648) ; 648: Jan 1947, Feb 1947, ..., Dec 2000 dummy=fltarr(4,3024) dummy2=fltarr(2,301) openr,2,'monthssn.dat' readf,2,dummy & month_ssn(*)=dummy(2,*) & close,2 openr,2,'yearssn.dat' readf,2,dummy2 & year_ssn(*)=dummy2(1,*) & close,2 openr,2,'solar47-00.mon' readf,2,month_f107 & close,2 print,'_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/' print,' ' print,'Data arrays:' print,' ' print,' month_ssn(3024) <== Monthly mean sun spot number' print,' 3024: Jan 1749, Feb 1749, ..., Dec 2000' print,' year_ssn(301) <== Annually averaged sun spot number' print,' 301: Year 1700, 1701, ..., 2000' print,' month_f107(648) <== Monthly mean f10.7 cm solar flux' print,' 648: Jan 1947, Feb 1947, ..., Dec 2000' print,' ' set_plot,'ps' device,/landscape !p.font=1 xstf=['1947',' ',' ',$ '50',' ','52',' ','54',' ','56',' ','58',' ',$ '60',' ','62',' ','64',' ','66',' ','68',' ',$ '70',' ','72',' ','74',' ','76',' ','78',' ',$ '80',' ','82',' ','84',' ','86',' ','88',' ',$ '90',' ','92',' ','94',' ','96',' ','98',' ',$ '00',' '] q=month_ssn(2376:2376+647) q0=q plot,q,/xs,/ys,yrange=[0,300],min_value=0.01,$ thick=2,xthick=4,ythick=4,$ charsize=1.3,xcharsize=1.3,ycharsize=1.3,$ xtitle='Year',xticks=54,xtickname=xstf,$ ytitle='F10.7 (10^-22 W m^-2 Hz^-1) or Sun Spot Number',$ title='These 1947-2000 Data Are Used In The Regression Model' oplot,month_f107,thick=7,min_value=0.01 oplot,[290,320],[270,270],thick=7 oplot,[290,320],[250,250],thick=3 xyouts,325,270,'Monthly Solar F10.7',font=1,charsize=1.5 xyouts,325,250,'Monthly Sun Spot Number',font=1,charsize=1.5 ;REGRESSION MODEL FIT: T=month_f107 IERROR=1 N=648 & M=4 F=FLTARR(N,M) TMIN=0.01 & TMAX=300 FOR I=1,N DO BEGIN F(I-1,0)=1. F(I-1,1)=q(i-1) F(I-1,2)=q(i-1)^2 F(I-1,3)=q(i-1)^3 ENDFOR ZREGR,IERROR,N,M,T,TMIN,TMAX,F,WK,A,AINV,S,C,RESID,COV print,'REGRESSION COEFFICIENTS (SEE ROUTINE ZREGR.PRO):' print,'C(0):',C(0) print,'C(1):',C(1) print,'C(2):',C(2) print,'C(3):',C(3) print,'COVARIANCE MATRIX "COV":' print,COV ;SCATTER DIAGRAM WITH REGRESSION CUBIC POLYNOMIAL FIT: plot,q0,month_f107,/xs,/ys,xrange=[0,300],yrange=[0,300],psym=2,symsize=0.5,$ thick=3,xthick=4,ythick=4,$ charsize=1.3,xcharsize=1.3,ycharsize=1.3,$ xtitle='Monthly Sun Spot Number',$ ytitle='Monthly F10.7 (10^-22 W m^-2 Hz^-1)',$ title='Scatter Plot Of The F10.7 And SSN Monthly Data For Years 1947-2000' for i=0,299 do begin x1=float(i) & x2=i+1. y1=c(0)+c(1)*x1+c(2)*x1*x1+c(3)*x1*x1*x1 y2=c(0)+c(1)*x2+c(2)*x2*x2+c(3)*x2*x2*x2 oplot,[x1,x2],[y1,y2],thick=5 endfor xyouts,160,150,'Regression Model:',charsize=1.7,font=1 xyouts,120,120,'F10.7 = C0 + C1*SSN + C2*SSN^2 + C3*SSN^3',charsize=1.5,font=1 xyouts,130,90,'C0 = 65.2',charsize=1.5,font=1 xyouts,130,80,'C1 = 0.645',charsize=1.5,font=1 xyouts,130,70,'C2 = 0.00298',charsize=1.5,font=1 xyouts,130,60,'C3 = -8.88E-06',charsize=1.5,font=1 ;PLOT ORIGINAL AND REGRESSION MODELED F10.7 FOR 1947-2000 plot,T,/xs,/ys,yrange=[0,300],min_value=0.01,$ thick=1,xthick=4,ythick=4,$ charsize=1.3,xcharsize=1.3,ycharsize=1.3,$ xtitle='Year',xticks=54,xtickname=xstf,$ ytitle='F10.7 (10^-22 W m^-2 Hz^-1)',$ title='The Original And Modeled F10.7 Data Are Remarkably Similar' oplot,T-RESID,thick=6,min_value=0.01,linestyle=1 oplot,[290,320],[270,270],thick=1 oplot,[290,320],[250,250],thick=6,linestyle=1 xyouts,325,270,'Original F10.7 Time Series',font=1,charsize=1.4 xyouts,325,250,'Regression Model Fit',font=1,charsize=1.4 plot,RESID,/xs,/ys,yrange=[-150,150],min_value=-50,max_value=50,$ thick=2,xthick=4,ythick=4,$ charsize=1.3,xcharsize=1.3,ycharsize=1.3,$ xtitle='Year',xticks=54,xtickname=xstf,$ ytitle='F10.7 Model Residual (10^-22 W m^-2 Hz^-1)',$ title='These Are The Month-To-Month Errors In The Modeled Time Series' oplot,[0,647],[0,0],thick=4 ;REGRESSION F10.7 1750-2000 MONTHLY MEANS: q=month_ssn f107_regress=fltarr(3024) for mon=0,3023 do begin x1=month_ssn(mon) f107_regress(mon)=c(0)+c(1)*x1+c(2)*x1*x1+c(3)*x1*x1*x1 endfor b=where(month_ssn lt 0.01) f107_regress(b)=0 !p.multi=[0,1,2] xstf=['Jan1750',' ',' ',' ',' ',$ 'Jan1800',' ',' ',' ',' ',$ 'Jan1850',' ',' ',' ',' ',$ 'Jan1900',' ',' ',' ',' ',$ 'Jan1950',' ',' ',' ',' ',$ 'Jan 2001'] plot,f107_regress(1:*),/xs,/ys,min_value=0.01,thick=2,xthick=3,ythick=3,$ charsize=1.1,xcharsize=1.1,ycharsize=1.1,yrange=[0,300],$ xticks=25,xtickname=xstf,$ xtitle='Month',ytitle='Model F10.7 (10^-22 W m^-2 Hz^-1)',$ title=$ 'Monthly Mean F10.7 For 1750-2000 Using Coefficients C0,C1,C2,C3' ;openw,2,'f107_regress_monthly.dat' ;for i=1,3023 do begin ; printf,2,f107_regress(i),format='(f6.1)' ;endfor ;close,2 ;REGRESSION F10.7 1700-2000 ANNUAL MEANS: q=year_ssn f107_regress=fltarr(301) for iyear=0,300 do begin x1=year_ssn(iyear) f107_regress(iyear)=c(0)+c(1)*x1+c(2)*x1*x1+c(3)*x1*x1*x1 endfor b=where(year_ssn lt 0.01) f107_regress(b)=0 !p.multi=[0,1,2] xstf=['1700','10','20','30','40','50','60','70','80','90',$ '1800','10','20','30','40','50','60','70','80','90',$ '1900','10','20','30','40','50','60','70','80','90',$ '2000'] plot,f107_regress,/xs,/ys,min_value=0.01,thick=5,xthick=3,ythick=3,$ charsize=1.1,xcharsize=1.1,ycharsize=1.1,yrange=[0,300],$ xticks=30,xtickname=xstf,$ xtitle='Year',ytitle='Model F10.7 (10^-22 W m^-2 Hz^-1)',$ title=$ 'Annual Mean F10.7 For 1700-2000 Using Coefficients C0,C1,C2,C3' ;openw,2,'f107_regress_annual.dat' ;for i=0,300 do begin ; printf,2,f107_regress(i),format='(f6.1)' ;endfor ;close,2 device,/close end