pro GS_trends_test restore, 'C:\MyProjects\20170711_EOVSA\ovsainput_20170711.sav' ;, /verb cd,'C:\Users\gfleishm\Dropbox\Magnet School\Lecture7(2017)', CURRENT=old_dir ;cd,'C:\Users\Gelu\Dropbox\Active Projects\Magnet School\Lectures\Lecture7(2017)', CURRENT=old_dir ParmIn=fltarr(29) n=150 ;number of frequencies s=fltarr(5, n) ;s[0, *] - frequency ;s[1, *]=observed (at the Earth) intensity, O-mode ;s[2, *]=exp(-tau) O-mode ;s[3, *]=observed (at the Earth) intensity, X-mode ;s[4, *]=exp(-tau) X-mode ParmIn[0] =(0.72*725d5*10)^2 ;Area, cm^2 arcsec 11 ParmIn[1] =725d5*10*2 ;Depth, cm arcsec 11 ParmIn[2] =2.e7 *1e0 ;T_0, K ParmIn[3] =0.07+0.03 ;\eps (not used in this example) ParmIn[4] =4.2 ;\kappa (not used in this example) ParmIn[5] =32 ;16 ;number of integration nodes ParmIn[6] =0.02* 1 ;E_min, MeV ParmIn[7] =10. *1 ;E_max, MeV ParmIn[8] =1.0 ;E_break, MeV (not used in this example) ParmIn[9] =4.25 + 0;4.0 ;\delta_1 ParmIn[10]=6.0 ;\delta_2 (not used in this example) ParmIn[11]=2d10 *1;1d11*7 ;n_0 - thermal electron density, cm^{-3} ParmIn[12]=0.8e7 *1 ;n_b - nonthermal electron density, cm^{-3} ParmIn[13]=170 *1.2 ;B - magnetic field, G ParmIn[14]=80 ;-50 ;theta - the viewing angle, degrees ParmIn[15]=5e8 ;starting frequency to calculate spectrum, Hz ParmIn[16]=0.02 ;logarithmic step in frequency ParmIn[17]= 3 ;distribution over energy (PLW is chosen) 5 th+ac; 2 th; 3 ac ParmIn[18]= n ;number of frequencies (specified above) ParmIn[19]=1 ;4 ;distribution over pitch-angle (GLC is chosen) ParmIn[20]=60 ;loss-cone boundary, degrees ParmIn[21]=50+10 ;beam direction (degrees) in GAU and SGA (not used in this example) ParmIn[22]=0.4 ;\Delta\mu ParmIn[23]=1 ;a_4 in SGA (not used in this example) ParmIn[25]= 0 ;f^C_cr, smooth 0 ParmIn[26]=ParmIn[25] ;f^WH_cr ParmIn[27]=1 ;matching on ParmIn[28]=2 ;Q-optimization on (with 2 bisection steps before including the term \ln Q) ;libname='libGS_Std_HomSrc_CEH.dll' ;library name (MS Windows variant) ;libname='c:\MyFortran11_Projects\libGS_Std_HomSrc_C_T\libGS_Std_HomSrc_C_T\x64\Release\libGS_Std_HomSrc_C_64.dll' libname='libGS_Std_HomSrc_C_64.dll' res=call_external(libname, 'GET_MW', ParmIn, s, /f_value) ;calling the GS code nu=s[0, *] ;array of the emission frequencies, GHz ;Default spectrum: I_O=s[1, *] ;array of the O-mode intensities, sfu I_X=s[3, *] ;array of the X-mode intensities, sfu ;Test spectrum kk=0 ; 0 area; 1 depth; 2 T; 6 E_min; 7 E_max; [8 E_break for DPL]; 9 delta_1; [10 delta_2 for DPL]; 11 n_therm; 12 n_b; 13 B; 14 theta ParmIn[kk]= Parmin[kk]*0.8 ;ParmIn[13]=Parmin[13]*1.3 ;ParmIn[12]=Parmin[12]*0.5 res=call_external(libname, 'GET_MW', ParmIn, s, /f_value) ;calling the GS code I_Otest=s[1, *] ;array of the O-mode intensities, sfu I_Xtest=s[3, *] ;array of the X-mode intensities, sfu ;freq224=[2.3,2.6,2.8,3.2,3.6,4.2,4.8,5.6,6.6,7.8,8.7,10.1,13.2,15.7,19.9,22.9] ;freqNoRP=[1.0,2.0,3.75,9.4,17,35,80] ;sigNoRP=[119.,73., 74.,254., 612., 754., 271.] ;sig224=[69.,80.,76.,79.,87.,94.,107.,135.,171.,235.,285.,335.,549.,670.,869.,982.] ;polNoRP=[ -0.004, 0.045, -0.16, -0.05, -0.03, 0.14] ;pol224=[ -0.04,-0.058,-0.067,-0.07, -0.062, -0.11, -0.16, -0.15, -0.13, -0.085, -0.092, -0.098, -0.099,-0.043, -0.0037, -0.0064] ;degreeNoRP=polNoRP/sigNoRP[0:5] ;degree224=pol224/sig224 ;----------------------------------------------------------------- ;now plotting the spectrum and polarization: maxy=2*max(I_X+I_O) window, 0, title='Intensity';, retain=2 plot, nu, I_X+I_O, /xlog, /ylog, xr=[1,100],/xst,/yst, $ xtitle='!8f!3, GHz', $ ytitle='!8I!3, sfu',linest=0, yrange=[1e-2, maxy] ;oplot,freqNoRP,sigNORP,psym=4 oplot, nu, I_Xtest+I_Otest, thick=2 oplot, freq, flux, psym=2, thick=2 ;xyouts,freq224,sig224,'x' ;negative window,1, title='Polarization';, retain=2 plot, nu,(I_X-I_O)/(I_X+I_O), /xlog, yrange=[-.3, .3],xr=[1,80],$ xtitle='!8f!3, GHz', $ ytitle='!7g!3' oplot, nu,(I_Xtest-I_Otest)/(I_Xtest+I_Otest), thick=2 ;oplot,freqNoRP[0:5],polNoRP,psym=4 ;oplot,freqNoRP[0:5],freqNoRP[0:5]*0,linest=3 ;xyouts,freq224,pol224,'x' ;negative cd,old_dir end