;============================================================================= ;Programers: A. Quivers (Nat. Space Scholars Prog.) ; D. Berdichevsky (L3GSI/SSAI) and ; T. Narock (L3GSI/SSAI) ; ;Responsible NASA Representative: Dr. Ronald P. Lepping ; Ronald.P.Lepping@nasa.gov ; ;Organization: Code 696, NASA/GSFC - Laboratory for Extraterrestrial Physics ; ;Application Name: Interpolation for Magnetic Cloud Parameter Error ; Estimation Program ;============================================================================= ;Editing Information: July 30, 2003 version by AQ ;Revisions: dbb Sept 18, 2003; dbb Oct 11, 2003 ; ;Revisions TWN Oct-Nov, 2003 ;Quadratic Fitting Revision TWN Nov 14, 2003 ;Added calculation of new Sigma B0 TWN Nov 18, 2003 ;============================================================================= function noise_conversion, chi_input snoise = (155/8.3666)*chi_input return, snoise end ;============================================================================= Function avg_1, arr_3 n_arr3=N_ELEMENTS(arr_3) RETURN, TOTAL(arr_3)/FLOAT(n_arr3) end ;============================================================================= Function m_calc, arr_2 n_arr2=N_ELEMENTS(arr_2)/2 arr_2o= FLTARR(n_arr2-1) ;---------------------------------------------------- Midpoint Calcuation ---- FOR index_2=0, n_arr2-2 DO BEGIN ala=(arr_2(index_2+1,1)-arr_2(index_2,1))/(arr_2(index_2+1,0)-$ arr_2(index_2,0)) arr_2o(index_2)= ala ENDFOR RETURN, arr_2o end ;============================================================================= Function next_gen, arr_1 n_arr1=N_ELEMENTS(arr_1) arr_1o= FLTARR(n_arr1-1) ;---------------------------------------------------- Midpoint Calcuation ---- FOR index_1=0, n_arr1-2 DO BEGIN d = (arr_1(index_1)+arr_1(index_1+1))/2 arr_1o(index_1)=avg_1(d) ENDFOR RETURN, arr_1o END ;============================================================================= PRO fit_l, anal_arr, m, b, RMSE ;------------------------------------------- First Derivative Initializations n_anal=N_ELEMENTS(anal_arr)/2 n_delta_one=n_anal-1 delta_one=FLTARR(n_delta_one,2) delta_one(*,0)=next_gen(anal_arr(*,0)) delta_one(*,1)=m_calc(anal_arr) ;----------------------------------------------- 'm' Parameter Compilation --- m=avg_1(delta_one(*,1)) ;----------------------------------------------- 'b' Parameter Compilation --- anal_diff=FLTARR(n_anal) anal_diff=anal_arr(*,1)-(m*anal_arr(*,0)) b=avg_1(anal_diff) ;------------------------------------ Regresion Mean Square Error Compilation pre_RMSE=FLTARR(n_anal) pre_RMSE=anal_arr(*,1)-((m*anal_arr(*,0))+b) pre_RMSE=pre_RMSE^2 RMSE=avg_1(pre_RMSE) RETURN end PRO fit_q_tom, anal_arr, a, b, c, chisq coeff = poly_fit(anal_arr[*, 0], anal_arr[*, 1], 2) a = coeff[2] b = coeff[1] c = coeff[0] pre_chi= anal_arr[*,1]-( (anal_arr[*,0]^2)*a + anal_arr[*,0]*b + c ) pre_chisq= pre_chi^2 chi=avg_1(pre_chi) chisq=avg_1(pre_chisq) END ;============================================================================= PRO fit_q, anal_arr,a,b,c,chisq,chi ;-------------------------------------------- Delta length def --------------- n_anal=N_ELEMENTS(anal_arr)/2 n_delta_one= n_anal-1 n_delta_two= n_delta_one-1 ;----------------------------------------------------- Delta Setup ----------- delta_one= FLTARR(n_delta_one,2) delta_two= FLTARR(n_delta_two,2) delta_one(*,0)=next_gen(anal_arr(*,0)) delta_two(*,0)=next_gen(delta_one(*,0)) delta_one(*,1)= m_calc(anal_arr) delta_two(*,1)= m_calc(delta_one) ;----------------------------------------------- Fit Parameter Calculation --- a_2=avg_1(delta_two(*,1)) a=a_2/2. delta_one_diff=delta_one(*,1)-(delta_one(*,0)*a_2) b= avg_1(delta_one_diff) anal_diff=anal_arr(*,1)-((anal_arr(*,0)^2)*a+anal_arr(*,0)*b) c=avg_1(anal_diff) ;---------------------------------------------- Error Parameter Calculation -- pre_chi= anal_arr(*,1)-((anal_arr(*,0)^2)*a+anal_arr(*,0)*b+c) pre_chisq= pre_chi^2 chi=avg_1(pre_chi) chisq=avg_1(pre_chisq) RETURN end ;============================================================================= PRO series_1_anal, unit, unit2, n_v, n_CA, n_beta, v, close, beta, sigma, $ para_name, sigma_2, s_v, chi_in, s_beta, s_CA, b0, $ newb0 = newb0 ;---------------------------------------- Manual Quadratic Compilation Block fit_q_compilaiton: a=FLTARR(n_CA,n_beta) b=FLTARR(n_CA,n_beta) c=FLTARR(n_CA,n_beta) chisq=FLTARR(n_CA,n_beta) anal_send=FLTARR(n_v,2) ; THE INDEPENDENT VARIABLE dbb, 10-11-03 anal_send[0,0] = v[0] ; dbb, 10-11-03 anal_send[1,0] = v[1] ; dbb, 10-11-03 anal_send[2,0] = v[2] ; dbb, 10-11-03 anal_send[3,0] = v[3] ; dbb, 10-11-03 anal_send[4,0] = v[4] ; dbb, 10-11-03 FOR beta_i=0,n_beta-1 DO BEGIN FOR CA_i=0,n_CA-1 DO BEGIN ; THE DEPENDENT VARIABLE dbb, 10-11-03 anal_send[0,1] = sigma[0,CA_i,beta_i] ; dbb, 10-11-03 anal_send[1,1] = sigma[1,CA_i,beta_i] ; dbb, 10-11-03 anal_send[2,1] = sigma[2,CA_i,beta_i] ; dbb, 10-11-03 anal_send[3,1] = sigma[3,CA_i,beta_i] ; dbb, 10-11-03 anal_send[4,1] = sigma[4,CA_i,beta_i] ; dbb, 10-11-03 fit_q_tom, anal_send, as, bs, cs, chisqs a[CA_i,beta_i]=as b[CA_i,beta_i]=bs c[CA_i,beta_i]=cs chisq[CA_i,beta_i]=chisqs ENDFOR ENDFOR ;----------------------- Quadratic Fit Parameter Display to Terminal Screen -- plot_stat_display: PRINTF,UNIT,'================================================================' PRINTF,UNIT,'--------------------------- Quadratic Fit Parameters (Series 1) ' PRINTF,UNIT,'================================================================' PRINTF,UNIT,'' PRINTF,UNIT,'Quadratic Fit Parameters For Parameter:',para_name FOR beta_i=0,n_beta-1 DO BEGIN PRINTF,UNIT,'_______________---____________________________________________' PRINTF,UNIT,'' PRINTF,UNIT,'Beta Angle:',beta[beta_i], 'Deg' FOR CA_i=0,n_CA-1 DO BEGIN PRINTF,UNIT,'' PRINTF,UNIT,'Closest Approch:',close(CA_i) PRINTF,UNIT,'a=',a(CA_i,beta_i),'b=',b(CA_i,beta_i),'c=',c(CA_i,beta_i),$ 'chisq=',chisq(CA_i,beta_i) ENDFOR ENDFOR ;-------------------------------------------- Curve Array Processor -------- fit_plot_1a: cont=100 vp=INDGEN(cont) vp=FLOAT(vp) vmax=v(n_v-1) vp=(vp/FLOAT(cont-1))*vmax sigmap=FLTARR(cont,n_CA,n_beta) FOR beta_i=0,n_beta-1 DO BEGIN FOR CA_i=0,n_CA-1 DO BEGIN sigmap(*,CA_i,beta_i)=a(CA_i,beta_i)*(vp^2)+b(CA_i,beta_i)*vp+c(CA_i,beta_i) ENDFOR ENDFOR ;---------------------------------------------- Series 1 Plot Intiator ------- fit_plot_1b: ;; TWN 10/20/04 s_v=noise_conversion(chi_in) ; dbb, 10-11-03 s_v=noise_conversion(chi_in) ; dbb, 10-11-03 sigma_2=FLTARR(n_CA,n_beta) sigma_2=a*(s_v^2)+b*s_v+c ;---------------------------------------------- Plot Solutions Display ------- solution_display: PRINTF,UNIT,'===============================================-================' PRINTF,UNIT,'----------------------------- Sigmas of Selected Chi (Series 1) ' PRINTF,UNIT,'================================================================' PRINTF,UNIT,'' PRINTF,UNIT,'Series 1 sigmas for parameter:',para_name PRINTF,UNIT,'Selected Chi:',chi_in,'noise:',s_v,'n_T' FOR beta_i=0,n_beta-1 DO BEGIN PRINTF,UNIT,'____---_______________________________________________________' PRINTF,UNIT,'' PRINTF,UNIT,'Beta Angle:',beta[beta_i],'Deg' FOR CA_i=0,n_CA-1 DO BEGIN PRINTF,UNIT,'' PRINTF,UNIT,'Closest Approch:',close(CA_i) PRINTF,UNIT,'sigma(v)=',sigma_2[CA_i,beta_i] ENDFOR ENDFOR final_sigma0 = sigma_2[0, 0] + ( ( sigma_2[0, 2]-sigma_2[0, 0] ) * $ (s_CA/0.6) ) final_sigma1 = sigma_2[2, 0] + ( ( sigma_2[2, 2]-sigma_2[2, 0] ) * $ (s_CA/0.6) ) final_sigma = final_sigma0 + ( (final_sigma1-final_sigma0) * $ ((s_beta-90)/(30-90)) ) CASE para_name OF 'SigB0': label = 'Sigma(B0) (nT) ' 'SigR0': label = 'Sigma(R0) (x10^-2 AU) ' 'Sigt0': BEGIN label = 'Sigma(ASF) (%) ' final_sigma = (final_sigma/12.5)*100 END 'SigCA': BEGIN label = 'Sigma(CA) (%) ' final_sigma = final_sigma*100 END 'SigThetaA': label = 'Sigma(ThetaA) (degrees) ' 'SigCone': label = 'Sigma(Beta(cone angle)) (degrees) ' 'SigAbsThtA': label = 'Sigma(|ThetaA|) (degrees) ' ENDCASE printf, unit2, label + string(final_sigma, format = '(F6.3)') IF keyword_set(newb0) THEN BEGIN new_sigmab0, final_sigma, b0, new_sigma printf, unit2, 'Sigma(B0 (new)) (nT) ' + string(new_sigma, $ format ='(F6.3)') ENDIF ;============================================================================= ;----------------------------------------------------- Chi Input Confirmation ;============================================================================= confirm_noise: RETURN end PRO new_sigmab0, oldb0, b0, newb0 newb0 = oldb0*(b0/16.4) END ;;;;;; MAIN dbb 10-13-03 PRO interpolation_mc, cloud_num, b0, chi_in, s_CA, s_beta, file file1 = file + '.output' file2 = file + '.summary' openw, unit, file1, /get_lun openw, unit2, file2, /get_lun printf, unit2, 'Cloud Number: ' + strtrim(cloud_num, 2) printf, unit2, 'Chi: ' + strtrim(string(chi_in, $ format = (F5.3)'), 2) printf, unit2, 'Closest Approach: ' + strtrim(string(s_CA, $ format = '(F4.2)'), 2) printf, unit2, 'Beta (degrees): ' + strtrim(s_beta, 2) printf, unit2, '' n_v = 5 n_CA = 4 n_beta = 3 v=FLTARR(n_v) close=FLTARR(n_CA) beta=FLTARR(n_beta) sigma=FLTARR(n_v,n_CA,n_beta) ;------------------- definition of the 3 independent variables ---dbb 10-13-03 v[0] = 0.0 ; noise dbb 10-13-03 v[1] = 0.5 ; dbb 10-13-03 v[2] = 2.0 ; dbb 10-13-03 v[3] = 3.0 ; dbb 10-13-03 v[4] = 4.0 ; dbb 10-13-03 close[0] = 0.0 ; CA dbb 10-13-03 close[1] = 0.3 ; dbb 10-13-03 close[2] = 0.6 ; dbb 10-13-03 close[3] = 0.9 ; dbb 10-13-03 beta[0] = 90 ; Cone-angle dbb 10-13-03 beta[1] = 60 ; dbb 10-13-03 beta[2] = 30 ; dbb 10-13-03 ;--------------------------------- Call Common Data Procedure ---dbb 10-13-03 ;--------------------------------- for the dependent variable ---dbb 10-13-03 sigcommon, sigB0, sigR0, sigt0, sigCA, sigthetA, sigCONE, sigAbsthtA ;------------------------------------for sigB0 Series One Procedure --------- sigma = sigB0 series_1_anal, unit, unit2, n_v, n_CA, n_beta, v, close, beta, sigma, $ 'SigB0', sigma_2, s_v, chi_in, s_beta, s_CA, b0, /newb0 ;------------------------------------for sigR0 Series One Procedure --------- sigma = sigR0 series_1_anal, unit, unit2, n_v, n_CA, n_beta, v, close, beta, sigma, $ 'SigR0', sigma_2, s_v, chi_in, s_beta, s_CA ;---------------------------------------for sigt0 Series One Procedure ------ sigma = sigt0 series_1_anal, unit, unit2, n_v, n_CA, n_beta, v, close, beta, sigma, '' + $ 'Sigt0', sigma_2, s_v, chi_in, s_beta, s_CA ;---------------------------------------for sigCA Series One Procedure ------ sigma = sigCA series_1_anal, unit, unit2, n_v, n_CA, n_beta, v, close, beta, sigma, $ 'SigCA', sigma_2, s_v, chi_in, s_beta, s_CA ;---------------------------------------for sigthetA Series One Procedure --- sigma = sigthetA series_1_anal, unit, unit2, n_v, n_CA, n_beta, v, close, beta, sigma, $ 'SigThetaA', sigma_2, s_v, chi_in, s_beta, s_CA ;---------------------------------------for sigCONE Series One Procedure ---- sigma = sigCONE series_1_anal, unit, unit2, n_v, n_CA, n_beta, v, close, beta, sigma, $ 'SigCone', sigma_2, s_v, chi_in, s_beta, s_CA ;---------------------------------------for sigAbsthtA Series One Procedure - sigma = sigAbsthtA series_1_anal, unit, unit2, n_v, n_CA, n_beta, v, close, beta, sigma, $ 'SigAbsThtA', sigma_2, s_v, chi_in, s_beta, s_CA printf, unit2, '' printf, unit2, 'Version date: Nov. 18, 2003' free_lun, unit2 free_lun, unit end