c c=============================================================| c Comparison of the UV continuum of symbiotic binaries | c by | c stellar and nebular components of radiation | c (Eq. (1)): | c F_obs = kh*B(Th)*exp[-NH*sigma_Ray] + kn*e(Te,a(HeII)) | c where | c model parameters are: kh, Th, NH, kn, Te, a(HeII)) | c | c-------------------------------------------------------------| c Input: lambda(A), F_obs(erg/s/cm2/A), F_err(erg/s/cm2/A) | c file, e.g., 'example.in' | c Input: ranges of parameters: begin_par, end_par, step_par | c prepare the file, 'ranges' | c Input: oscilator strengths of the Lyman series H lines | c file, 'l_f' | c Output: grid of all model parameters | c give the name of the file, e.g., 'output1' | c Output: grid of the 200 best model parameters | c give the name of the file, e.g., 'output2' | c Output: selected SED model - the result, | c give the name of the file, e.g., 'model' | c=================================================================| c version v2: rangers of model parameters from the file 'ranges' | c Augustin Skopal, May 10, 2021 | c=================================================================| c program uvsed_v2 implicit double precision(a-h,o-z) character*30 in1,out1,out2,out3 character*2 yesno character*70 head1,head2,head3,head4 c dimension x(3000),f_obs(3000),err(3000),coef(3000) dimension all(40),fij(40),flux(3000),e_hheii(3000) dimension a(1000000),b(1000000),c(1000000),d(1000000) dimension e(1000000),f(1000000),g(1000000) common /a/eion,cc,ak,h data tcs/6.65246E-25/ c head1='# No. xi^2r NH a(HeII) Th' head2=' kh Te kn' head3='# xi^2r NH a(HeII) Th' head4=' kh Te kn' c c = H I hh1=9.12d-6 hh2=3.646d-5 hh3=8.204d-5 hh4=1.459d-4 hh5=2.2795d-4 hh6=3.2825d-4 hh7=4.4677d-4 hh8=5.8354d-4 hh9=7.3854d-4 hh10=9.1178d-4 hh11=1.10326d-3 hh12=1.31297d-3 c = HeII hhe1=2.28d-6 hhe2=9.12d-6 hhe3=2.051d-5 hhe4=3.646d-5 hhe5=5.698d-5 hhe6=8.204d-5 hhe7=1.1169d-4 hhe8=1.4590d-4 hhe9=1.8464d-4 hhe10=2.2795d-4 hhe11=2.7582d-4 hhe12=3.2825d-4 hhe13=3.8523d-4 hhe14=4.4677d-4 hhe15=5.1288d-4 hhe16=5.8354d-4 hhe17=6.5876d-4 hhe18=7.3854d-4 hhe19=8.2288d-4 hhe20=9.1178d-4 hhe21=1.00524d-3 hhe22=1.10326d-3 hhe23=1.20583d-3 hhe24=1.31297d-3 c-------------------------- c8=1.0d-8 c1=2.160259e-32 c2=6.84194e-38 c11=1.191d-5 c22=1.4388d0 eion=2.17866d-11 ! ionization of H (=13.598 eV) ak=1.380622e-16 h=6.626196e-27 cc=2.997925e+10 amax=1.0d6 ! maximum number of models c ate=15001.0d0 ! a correction for e(He++) coefficient z12=1.5 c 15 format(a30) c c Output: xi^2r NH abun Th kh Te kn 160 format(e12.4,1x,e11.4,1x,f5.2,f8.0,e11.3,f8.0,e12.4) c j xi^2r NH abun Th kh Te kn 161 format(i4,e12.4,1x,e11.4,1x,f5.2,f8.0,e11.3,f8.0,e12.4) c Output SED: 150 format(f8.1,1x,3e12.4) c c==================================================! c----- Testing points: lambda, flux, flux_err -----! c--------------------------------------------------! c write(*,16) 16 format(1x,'Input - file with the measured fluxes: ',$) read(*,15)in1 open(1,file=in1,status='old') c ndat=0 do 25 i=1,3000 read(1,*,end=29)x(i),f_obs(i),err(i) f_obs(i)=f_obs(i)/1E-13 ! in 1E-13 err(i)=err(i)/1E-13 ! in 1E-13 x(i)=c8*x(i) coef(i)=1.d0 25 continue c 29 close(1) ndat=i-1 ndf=ndat-5 ! number of degree of freedom write(*,'(2x,a,i5)')'Number of flux points = ',ndat c write(*,'(2x,a,i5)')'Number of n.d.f. = ',ndf c c========================================! c---- reading the [lambda,f_ij] file ----! c----------------------------------------! c open(2,file='l_f',status='old') c do 11 i=1,40 read(2,*,end=10)all(i),fij(i) ! reading file 'l_f' all(i)=c8*all(i) 11 continue 10 close(2) nn=i-1 c c========================================! c-----> Input of the "n_H" range --------! c-----> Input of the "abun" range -------! c-----> Input of the "T_h" range -------! c-----> Input of the "k_h" range --------! c-----> Input of the "T_e" range --------! c----------------------------------------! c open(3,file='ranges',status='old') c i=1 50 read(3,*,end=13)par_b,par_e,par_step ! reading 'ranges' goto(51,52,53,54,55) i 51 dnh_b=par_b dnh_e=par_e d_nh=par_step pnh=(dnh_e-dnh_b)/d_nh+1.0 if(pnh.lt.2.0)pnh=1.0 i=i+1 goto 50 52 abun_b=par_b abun_e=par_e d_abun=par_step pabun=(abun_e-abun_b)/d_abun+1.0 if(pabun.lt.2.0)pabun=1.0 i=i+1 goto 50 53 th_b=par_b th_e=par_e dt_h=par_step pth=(th_e-th_b)/dt_h+1.0 if(pth.lt.2.0)pth=1.0 i=i+1 goto 50 54 ak_b=par_b ak_e=par_e dak=par_step pkh=(ak_e-ak_b)/dak+1.0 if(pkh.lt.2.0)pkh=1.0 i=i+1 goto 50 55 te_b=par_b te_e=par_e dt_e=par_step pte=(te_e-te_b)/dt_e+1.0 if(pte.lt.2.0)pte=1.0 13 close(3) c pmodels=pnh*pabun*pth*pkh*pte c write(*,'(a,f9.0)')'Number of models: ',pmodels if(pmodels.gt.amax)then write(*,'(a)')'Too dense grid: number of models > 1 milion' goto 1003 endif c c========================================! c-- Cycling for n_H, a, T_h, k_h, T_e ---! c----------------------------------------! c k=1 ! the first model SED c dnh=dnh_b-d_nh 70 dnh=dnh+d_nh ! "n_H" if(dnh.gt.dnh_e)goto 1000 ! -----> all SEDs calculated abun=abun_b-d_abun 7 abun=abun+d_abun ! "a" if(abun.gt.abun_e)goto 70 ! -----> new "n_H" t_h=th_b-dt_h 77 t_h=t_h+dt_h ! "t_h" if(t_h.gt.th_e)goto 7 ! -----> new "a" ak_h=ak_b-dak 777 ak_h=ak_h+dak ! "k_h" if(ak_h.gt.ak_e)goto 77 ! -----> new "t_h" c c-------------------------------------------! c--- flux(i) = k_h*B(x(i),T_h)*exp[-tau] ---! c===========================================! c c bb for T_h c ========== do 22 i=1,ndat x5=x(i)**5 ca=ak_h*c11/x5 cb=c22/x(i) help=dexp(cb/t_h)-1.d0 bb=ca/help c c exp[-n_H*sigma_R] c =================== suma=0 do 222 ii=1,nn aa=fij(ii)/((x(i)/all(ii))**2.-1.0d0) suma=suma+aa 222 continue suma2=suma*suma sigmar=tcs*suma2 tau=dnh*sigmar c flux(i)=bb*dexp(-tau) 22 continue c c----------------------------------------! c t_e=te_b-dt_e 7777 t_e=t_e+dt_e ! "t_e" if(t_e.gt.te_e)goto 777 ! -----> new "k_h" c c===================================================================! c-- Volume emission coefficient 'e_hheii(i)' at wavelength 'x(i)' --! c-------------------------------------------------------------------! c do 21 i=1,ndat al=x(i) z=1.0d0 101 if(z.lt.z12)then ! calculate el_hi ! if(al.le.hh1)an=1.0d0 if((al.gt.hh1).and.(al.le.hh2))an=2.0d0 if((al.gt.hh2).and.(al.le.hh3))an=3.0d0 if((al.gt.hh3).and.(al.le.hh4))an=4.0d0 if((al.gt.hh4).and.(al.le.hh5))an=5.0d0 if((al.gt.hh5).and.(al.le.hh6))an=6.0d0 if((al.gt.hh6).and.(al.le.hh7))an=7.0d0 if((al.gt.hh7).and.(al.le.hh8))an=8.0d0 if((al.gt.hh8).and.(al.le.hh9))an=9.0d0 if((al.gt.hh9).and.(al.le.hh10))an=10.0d0 if((al.gt.hh10).and.(al.le.hh11))an=11.0d0 if((al.gt.hh11).and.(al.le.hh12))an=12.0d0 if(al.gt.hh12)an=13.0d0 else if(al.le.hhe1)an=1.0d0 ! calculate el_heii ! if((al.gt.hhe1).and.(al.le.hhe2))an=2.0d0 if((al.gt.hhe2).and.(al.le.hhe3))an=3.0d0 if((al.gt.hhe3).and.(al.le.hhe4))an=4.0d0 if((al.gt.hhe4).and.(al.le.hhe5))an=5.0d0 if((al.gt.hhe5).and.(al.le.hhe6))an=6.0d0 if((al.gt.hhe6).and.(al.le.hhe7))an=7.0d0 if((al.gt.hhe7).and.(al.le.hhe8))an=8.0d0 if((al.gt.hhe8).and.(al.le.hhe9))an=9.0d0 if((al.gt.hhe9).and.(al.le.hhe10))an=10.0d0 if((al.gt.hhe10).and.(al.le.hhe11))an=11.0d0 if((al.gt.hhe11).and.(al.le.hhe12))an=12.0d0 if((al.gt.hhe12).and.(al.le.hhe13))an=13.0d0 if((al.gt.hhe13).and.(al.le.hhe14))an=14.0d0 if((al.gt.hhe14).and.(al.le.hhe15))an=15.0d0 if((al.gt.hhe15).and.(al.le.hhe16))an=16.0d0 if((al.gt.hhe16).and.(al.le.hhe17))an=17.0d0 if((al.gt.hhe17).and.(al.le.hhe18))an=18.0d0 if((al.gt.hhe18).and.(al.le.hhe19))an=19.0d0 if((al.gt.hhe19).and.(al.le.hhe20))an=20.0d0 if((al.gt.hhe20).and.(al.le.hhe21))an=21.0d0 if((al.gt.hhe21).and.(al.le.hhe22))an=22.0d0 if((al.gt.hhe22).and.(al.le.hhe23))an=23.0d0 if((al.gt.hhe23).and.(al.le.hhe24))an=24.0d0 if(al.gt.hhe24)an=25.0d0 endif c e_fb=0. 200 e_fb=e_fb+e_n(t_e,al,an,z) an=an+1.0d0 if(an.le.15.0)goto 200 c c----------------------------e_fb------------------------------------ c e_fb_ny=c1*e_fb*z**4 !............. erg cm+3 s-1 Hz-1 e_fb_l=c8*e_fb_ny*cc/al/al !............. erg cm+3 s-1 A-1 e_fb_cm=e_fb_l/c8 !..............erg cm+3 s-1 cm-1 al_aa=al/c8 c c----------------------------e_ff------------------------------------ c any=cc/al aid=dexp((-h*any)/(ak*t_e)) ttt=dsqrt(t_e) e_ff_ny=z*z*c2*aid/ttt !..............erg cm+3 s-1 Hz-1 e_ff_l=c8*e_ff_ny*cc/al/al !..............erg cm+3 s-1 A-1 e_ff_cm=e_ff_l/c8 !..............erg cm+3 s-1 cm-1 c c-------------------------------------------------------------------- c e_l=(e_fb_l+e_ff_l)/1.0d-28 if(z.lt.z12)then el_hi=e_l z=z+1.0d0 goto 101 else el_heii=e_l if((t_e.lt.ate).and.(al.lt.h4))el_heii=0.96d0*el_heii endif aaa=(1.0d0+2.0d0*abun) c--------------------------------------------------- e_hheii(i)=aaa*el_hi+abun*aaa*el_heii c--------------------------------------------------- c f_obs(i)-flux(i) = F_N = k_N x e_hheii(i) --> coef(i)=(f_obs(i)-flux(i))/e_hheii(i) ! = k_N c 21 continue c s_coef=0 do 23 i=1,ndat s_coef=s_coef+coef(i) 23 continue ak_n=s_coef/ndat ! = averaged k_N c chi2=0 do 24 i=1,ndat weight=err(i) chi2=chi2+((f_obs(i)-flux(i)-ak_n*e_hheii(i))/weight)**2 24 continue c chi2_red=chi2/ndf c c SED model parameters: a(k)=chi2_red b(k)=dnh c(k)=abun d(k)=t_h e(k)=ak_h f(k)=t_e g(k)=ak_n c k=k+1 ! the next model SED c goto 7777 ! goto for a new "T_e" ! c 1000 nmod=k-1 ! number of SED models write(*,'(/1x,a,i8)')'Number of SED models = ',nmod c c------------------------------------------------------------------------ c write(*,28) 28 format(' Write all models to a file? [y/n] ',$) read(*,30)yesno 30 format(a2) if(yesno.ne.'y')goto 31 c write(*,32) 32 format(1x,'Output - grid of all parameters: ',$) read(*,15)out1 open(4,file=out1,status='new') c write(4,'(a40,a30)')head3,head4 do 33 i=1,nmod write(4,160)a(i),b(i),c(i),d(i),e(i),f(i),g(i) 33 continue close(4) c c------------------------------------------------------------------ c c Arrarging SED models by 'xi^2_red' c ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ 31 iend=200 ! the firtst 200 SED models c do 26 i=1,iend do 26 j=nmod,i+1,-1 if(a(j).gt.a(j-1))goto 26 aa=a(j) bb=b(j) ccc=c(j) dd=d(j) ee=e(j) ff=f(j) gg=g(j) c a(j)=a(j-1) b(j)=b(j-1) c(j)=c(j-1) d(j)=d(j-1) e(j)=e(j-1) f(j)=f(j-1) g(j)=g(j-1) c a(j-1)=aa b(j-1)=bb c(j-1)=ccc d(j-1)=dd e(j-1)=ee f(j-1)=ff g(j-1)=gg c 26 continue c c write(*,35) 35 format(1x,'Output - grid of the 200 best model parameters: ',$) read(*,15)out3 open(7,file=out3,status='new') c write(*,'(1x,a40,a30)')head1,head2 write(7,'(1x,a40,a30)')head1,head2 c j=1 27 write(7,161)j,a(j),b(j),c(j),d(j),e(j),f(j),g(j) if(j.lt.11)write(*,161)j,a(j),b(j),c(j),d(j),e(j),f(j),g(j) j=j+1 if(j.gt.iend)goto 1001 goto 27 1001 close(7) c c A new grid? [yes/no] c------------------------------------------------------------------ c c Calculation of the selected model SED: c ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ write(*,18) 18 format(2x,'Calculate the SED model No.: ',$) read(*,'(i4)')im c chi2_red=a(im) dnh=b(im) abun=c(im) t_h=d(im) ak_h=e(im) t_e=f(im) ak_n=g(im) c write(*,19) 19 format(1x,'Output - file with the selected SED model: ',$) read(*,15)out2 open(6,file=out2,status='new') c c --- Writing SED parameters to the head of the file '5' --- c write(6,'(a,a30)')'# Measured flux points: ',in1 write(6,'(a,f8.4)')'# xi^2r =',chi2_red write(6,'(a,e12.4)')'# NH = ',dnh write(6,'(a,f6.2)')'# a(H++) = ',abun write(6,'(a,f8.0)')'# Th = ',t_h write(6,'(a,e12.4)')'# kh = ',ak_h write(6,'(a,f8.0)')'# Te = ',t_e write(6,'(a,e12.4)')'# kn = ',ak_n write(6,'(a)')'# SED = kh*$2 + kn*$4 (in 1E-13 erg/s/cm2/A)' write(6,'(a)')'# -------------------------------------------' write(6,'(a)')'# lambda B(Th)e(-tau) B(Th) e(HHeII)' write(6,'(a)')'# -------------------------------------------' c alb=c8*1000.0 ale=c8*5001.0 dal=c8*2.0d0 c c==============================================! c-- Nebula: em. coeff. 'e_hheii(i)' for T_e --! c----------------------------------------------! c al=alb-dal 88 al=al+dal if(al.gt.ale)goto 1002 ! stop, end c z=1.0d0 100 if(z.lt.z12)then ! calculate el_hi ! if(al.le.hh1)an=1.0d0 if((al.gt.hh1).and.(al.le.hh2))an=2.0d0 if((al.gt.hh2).and.(al.le.hh3))an=3.0d0 if((al.gt.hh3).and.(al.le.hh4))an=4.0d0 if((al.gt.hh4).and.(al.le.hh5))an=5.0d0 if((al.gt.hh5).and.(al.le.hh6))an=6.0d0 if((al.gt.hh6).and.(al.le.hh7))an=7.0d0 if((al.gt.hh7).and.(al.le.hh8))an=8.0d0 if((al.gt.hh8).and.(al.le.hh9))an=9.0d0 if((al.gt.hh9).and.(al.le.hh10))an=10.0d0 if((al.gt.hh10).and.(al.le.hh11))an=11.0d0 if((al.gt.hh11).and.(al.le.hh12))an=12.0d0 if(al.gt.hh12)an=13.0d0 else if(al.le.hhe1)an=1.0d0 ! calculate el_heii ! if((al.gt.hhe1).and.(al.le.hhe2))an=2.0d0 if((al.gt.hhe2).and.(al.le.hhe3))an=3.0d0 if((al.gt.hhe3).and.(al.le.hhe4))an=4.0d0 if((al.gt.hhe4).and.(al.le.hhe5))an=5.0d0 if((al.gt.hhe5).and.(al.le.hhe6))an=6.0d0 if((al.gt.hhe6).and.(al.le.hhe7))an=7.0d0 if((al.gt.hhe7).and.(al.le.hhe8))an=8.0d0 if((al.gt.hhe8).and.(al.le.hhe9))an=9.0d0 if((al.gt.hhe9).and.(al.le.hhe10))an=10.0d0 if((al.gt.hhe10).and.(al.le.hhe11))an=11.0d0 if((al.gt.hhe11).and.(al.le.hhe12))an=12.0d0 if((al.gt.hhe12).and.(al.le.hhe13))an=13.0d0 if((al.gt.hhe13).and.(al.le.hhe14))an=14.0d0 if((al.gt.hhe14).and.(al.le.hhe15))an=15.0d0 if((al.gt.hhe15).and.(al.le.hhe16))an=16.0d0 if((al.gt.hhe16).and.(al.le.hhe17))an=17.0d0 if((al.gt.hhe17).and.(al.le.hhe18))an=18.0d0 if((al.gt.hhe18).and.(al.le.hhe19))an=19.0d0 if((al.gt.hhe19).and.(al.le.hhe20))an=20.0d0 if((al.gt.hhe20).and.(al.le.hhe21))an=21.0d0 if((al.gt.hhe21).and.(al.le.hhe22))an=22.0d0 if((al.gt.hhe22).and.(al.le.hhe23))an=23.0d0 if((al.gt.hhe23).and.(al.le.hhe24))an=24.0d0 if(al.gt.hhe24)an=25.0d0 endif c e_fb=0. 201 e_fb=e_fb+e_n(t_e,al,an,z) an=an+1.0d0 if(an.le.15.0)goto 201 c c----------------------------e_fb------------------------------------ c e_fb_ny=c1*e_fb*z**4 !............. erg cm+3 s-1 Hz-1 e_fb_l=c8*e_fb_ny*cc/al/al !............. erg cm+3 s-1 A-1 e_fb_cm=e_fb_l/c8 !..............erg cm+3 s-1 cm-1 al_aa=al/c8 c c----------------------------e_ff------------------------------------ c any=cc/al aid=dexp((-h*any)/(ak*t_e)) ttt=dsqrt(t_e) e_ff_ny=z*z*c2*aid/ttt !..............erg cm+3 s-1 Hz-1 e_ff_l=c8*e_ff_ny*cc/al/al !..............erg cm+3 s-1 A-1 e_ff_cm=e_ff_l/c8 !..............erg cm+3 s-1 cm-1 c c-------------------------------------------------------------------- c e_l=(e_fb_l+e_ff_l)/1.0d-28 if(z.lt.z12)then el_hi=e_l z=z+1.0d0 goto 100 else el_heii=e_l if((t_e.lt.ate).and.(al.lt.hhe4))el_heii=0.96d0*el_heii endif aaa=(1.0d0+2.0d0*abun) c c------------------------------------------------------------- ehheii=aaa*el_hi+abun*aaa*el_heii ! = em. coeff. c------------------------------------------------------------- c c k_h*B(al,T_h) c ============= x5=al**5 ca=c11/x5 cb=c22/al help=dexp(cb/t_h)-1.d0 bb=ca/help ! = B(al,T_h) c c exp[-n_H*sigma_R] c ================= suma=0 do 36 ii=1,nn aa=fij(ii)/((al/all(ii))**2.-1.0d0) suma=suma+aa 36 continue suma2=suma*suma sigmar=tcs*suma2 tau=dnh*sigmar flux_ray=bb*dexp(-tau) ! B(al,T_h)exp(-tau) c c------------------------------------------------ write(6,150)al/c8,flux_ray,bb,ehheii c------------------------------------------------ c goto 88 c 1002 close(6) c stop 1003 end c c------------------------------------------------------------------------ c function e_n(t_e,al,an,z) implicit double precision(a-h,o-z) common /a/eion,cc,ak,h an2=an*an ! n**2 an3=an2*an ! n**3 zz=z*z ain=zz*eion/an2 ! I_n any=cc/al help=dexp((ain-h*any)/(ak*t_e)) tt=t_e**1.5 e_n=help/an3/tt return end c c========================================================================| c-- END -- END -- END -- END -- END -- END -- END -- END -- END -- END --| c========================================================================|