work.pro

用idl进行二项式拟合得例子,通过不同参数,可以拟合不同得二项式,还有ps与eps制图,适合初学者

pro work,f,f_out,M=m,P=ps,T_OUT=t_out
  ;manually or files
  if  KEYWORD_SET(m) then BEGIN
    print,'please enter the "Omega_M" "Lambda0" "H0" "Zmin" "Zmix"!'
    A1=1.0d-1
    A2=1.0d-1
    A3=1.0d-1
    A4=1.0d-1
    A5=1.0d-1
    Read,A1,A2,A3,A4,A5
    V1=A1
    V2=A2
    V3=A3
    Zmin=A4
    Zmix=A5
    
  endif else begin
    ;======================================================================================
    ;Enter the path and the name,the format is like "E:\xxx\xxx.xxx"
  
    infiles=strtrim(f,2)
    readcol,infiles,format='D,D,D,D,D',skipline=1,V1,V2,V3,V4,V5
    Zmix=V4
    Zmin=V5
    
  endelse
  file_name=f_out
 
  D=Dblarr(N_elements(V1),10000 )
  Z=dblarr(N_elements(V1),10000)
  For i=0,n_elements(V1)-1 do begin
    if (Zmin[i] LE 0) then begin
      Zmin[i]=1d-1
    endif
     
    Z[i,*]=Zmin[i]+ (Zmix[i]-Zmin[i])*Dindgen(10000)/9999
    D[i,*]=lumdist(Z[i,*],H0=100*V3[i],Omega_M=V1[i],Lambda0=V2[i])
  endfor
  x_min=min(Z)
  y_min=min(D)
  x_max=max(Z)
  y_max=max(D)
  ax=x_min-(x_min-x_max)*5d-2
  ay=10^[(alog10(y_max)+alog10(y_min)*0.1)/1.1]
  ;======================================================================
  if  KEYWORD_SET(ps) then BEGIN
  
    For i=0,n_elements(V1)-1 do begin
    
      set_plot,'ps'
      entry_device=!d.name
      device,filename=STRTRIM(file_name,2)+STRTRIM(string(i)+'.ps',2)
      device,xsize=16.0,ysize=16.0
      plot,Z[i,*],D[i,*],/ylog,xtitle='z',ytitle='!18d!l!17L!N !17(Mpc)!3',xrange=[x_min,x_max],yrange=[y_min,y_max],xstyle=1,ystyle=1
      xyouts,ax,ay,'!7(X!S!R!lK!N,!7X!S!R!LL!N,!8h!17)  =  ('+STRCOMPRESS(String(V1[i],format='(d5.2)'))+','+STRCOMPRESS(string(V2[i],format='(d5.2)'))+','+ STRCOMPRESS(String(V3[i],format='(d5.2)'))+')'
      device,/close_file
      set_plot,entry_device
    endfor
  ;==============================================================================
  endif else begin
    For i=0,n_elements(V1)-1 do begin
    
      entry_device=!d.name
      set_plot,'ps'
      DEVICE, /ENCAPSUL
      device,filename=STRTRIM(file_name,2)+STRTRIM(string(i)+'.eps',2)
      device,xsize=16.0,ysize=16.0
      plot,Z[i,*],D[i,*],/ylog,xtitle='z',ytitle='!18d!l!17L!N !17(Mpc)!3',xrange=[x_min,x_max],yrange=[y_min,y_max],xstyle=1,ystyle=1
      xyouts,ax,ay,'!7(X!S!R!lK!N,!7X!S!R!LL!N,!8h!17)  =  ('+STRCOMPRESS(String(V1[i],format='(d5.2)'))+','+STRCOMPRESS(string(V2[i],format='(d5.2)'))+','+ STRCOMPRESS(String(V3[i],format='(d5.2)'))+')'
      device,/close_file
      set_plot,entry_device
    endfor
  endelse
  
  ;=======================================================================
  ;text or not
  if KEYWORD_SET(t_out) then begin
    if KEYWORD_SET(m) then begin
      forprint,Z[0,*],D[0,*],TEXT=STRTRIM(file_name,2)+STRTRIM(string(1,format='(I5)'),2)+'.txt',format='(2(d15.6,15X))',comment=string('Redshift','Distance',format='(2(A15,15X))')
    endif else begin
      for j=0,N_ELEMENTS(t_out)-1 do begin
      
        forprint,Z[t_out[j],*],D[t_out[j],*],TEXT=STRTRIM(file_name,2)+STRTRIM(string(t_out[j],format='(I5)'),2)+'.txt',format='(2(d15.6,15X))',comment=string('Redshift','Distance',format='(2(A15,15X))')
      endfor
    endelse
  endif else begin
    return
  endelse
  
end