PRO EDD ;************************************************************************ ; Exercise 6-4 in Mihalas ; Solution to the transfer equation using the method of ; Variable Eddington Factors. ; To see the results just look at file 'results' after running program. ; By Dean Richardson, Myra Blaylock and Thomas Gower ; ; ; this code doesn't work correctly but it is a start. ; so don't use it but go ahead and use it as a template ; -e baron ;************************************************************************ close,12 close,13 close,14 close,15 close,16 ;openw, 12, 'tgrid3' ;openw, 13, 'coef3' openw, 14, 'ufile3' openw, 15, 'jfile3' openw,16, 'results' layer=50L tk=dblarr(layer) taup=dblarr(layer) taum=dblarr(layer) dtau=dblarr(layer) kloop=5 tminlg=-3.0 tmaxlg=1.0 ;-- set up tau grid--------------------------------------------- ;-- compute tau grid, if layer .ne. 50 .or. tminlg .lt. tmaxlg if(not (tminlg lt tmaxlg)) then begin print,'main: *err* : bad tauscale given!' print,'layer=',layer,' tminlg=',tminlg,' tmaxlg=',tmaxlg stop endif tstlg = (tmaxlg-tminlg)/(layer-3) tk(0) = 0.e0 ishift = 0 for i=1,layer-2 do begin tk(i+ishift) = 10.0^(tminlg+(i-2)*tstlg) if(tk(i) ge 1.0 and tk(i-1) lt 1.0) then begin tk(i+1) = tk(i) tk(i) = 1.0 ; -- fix bug for tk(i)=tk(i+1)=1: (phh, 27/oct/92) if(tk(i) eq tk(i+1)) then tk(i) = 0.50*(tk(i+1)-tk(i-1)) ; -- end bug fix ishift = 1 endif endfor ;--Determine tau_d-1/2, tau_d+1/2 and delta tau--------------------- for i=1,layer-2 do begin taup(i)=tk(i+1)-tk(i) taum(i)=tk(i)-tk(i-1) dtau(i)=0.5*(taum(i)+taup(i)) endfor for i=0,layer-3 do begin taup(i)=taup(i+1) taum(i)=taum(i+1) dtau(i)=dtau(i+1) ; printf, 12, tk(i),taup(i),taum(i),dtau(i) endfor ;------------------------------------------------------------------- au=dblarr(layer-2) bu=dblarr(layer-2) cu=dblarr(layer-2) aj=dblarr(layer-2) bj=dblarr(layer-2) cj=dblarr(layer-2) Lu=dblarr(layer-2) Lj=dblarr(layer-2) s=dblarr(layer-2) sp=dblarr(layer-2) ua=dblarr(layer-2) jn=dblarr(layer-2) fd=dblarr(layer-2) ;--Loop to change the number of angles------------------------------ for im=1,3 do begin if (im eq 1) then m=2 if (im eq 2) then m=4 if (im eq 3) then m=5 u=dblarr(m,layer-2) mu=dblarr(m) w=dblarr(m) ;---mu=cos(theta) and w=gausian quadrature weights------------------ if(m eq 2) then begin mu(0)=0.2113248654 mu(1)=0.7886751346 w(0)=0.5 w(1)=0.5 endif if(m eq 4) then begin mu(0)=0.0694318442 mu(1)=0.3300094782 mu(2)=0.6699905128 mu(3)=0.9305681558 w(0)=0.1739274226 w(1)=0.3260725774 w(2)=0.3260725774 w(3)=0.1739274226 endif if(m ne 2) and (m ne 4)then begin m=5 mu(0)=0.0469100770 mu(1)=0.2307653449 mu(2)=0.5000000000 mu(3)=0.7692346551 mu(4)=0.9530899230 w(0)=0.1184634425 w(1)=0.2393143352 w(2)=0.2844444444 w(3)=0.2393143352 w(4)=0.1184634425 endif ;--Set s=1, initially for all depths-------------------------------- for i=0,layer-3 do begin s(i)=1.0 sp(i)=10.0 endfor ;--Loop to change epsilon------------------------------------------- for iep=1,3 do begin if (iep eq 1) then ep=0.1 if (iep eq 2) then ep=0.01 if (iep eq 3) then ep=0.0001 ;--Big loop for the convergence of s-------------------------------- print, 'For m =', m, ' and ep =', ep printf, 16, 'For m =', m, ' and ep =', ep print, 'The columns are: k, u(1,D), J(D), f(D), and maxdiff' printf, 16, 'The columns are: k, u(1,D), J(D), f, and maxdiff' for k=0, kloop do begin printf,14, 'k =', k printf,15, 'k =', k ;--Determine coef. of u for tridiagonal matrix---------------------- for j=0,m-1 do begin for i=1,layer-4 do begin au(i)=(mu(j)^2)/(dtau(i)*taum(i)) cu(i)=(mu(j)^2)/(dtau(i)*taup(i)) bu(i)=au(i)+cu(i)+1 Lu(i)=s(i) endfor ;--Set boundary conditions------------------------------------------ au(0)=0.0 bu(0)=2*(mu(j)^2)/(taup(0)^2) + 2*mu(j)/taup(0) + 1 cu(0)=2*(mu(j)^2)/(taup(0)^2) Lu(0)=s(0) Lu(layer-3)=(1.0 + 2*mu(j)/taum(layer-3))*s(layer-3) au(layer-3)=2*(mu(j)^2)/(taum(layer-3)^2) bu(layer-3)=au(layer-3) + 2*mu(j)/taum(layer-3) + 1 cu(layer-3)=0.0 for i=0,layer-3 do begin ; printf, 13, au(i),bu(i),cu(i) endfor ;--Solve tridiagonal matrix and put solutions, u, into a----------- ; matrix where we have one column of depth values for each angle. ;------------------------------------------------------------------ ua=trisol(-au,bu,-cu,Lu) for i=0,layer-3 do begin u(j,i)=ua(i) if (j eq 1) then printf, 14, ua(i) endfor endfor ;--Determine Eddington factors and h.-------------------------------- for i=0,layer-3 do begin x1=0.0 x2=0.0 x3=0.0 for j=0,m-1 do begin x1=x1+w(j)*(mu(j)^2)*u(j,i) x2=x2+w(j)*u(j,i) if (i eq 0) then x3=x3+w(j)*mu(j)*u(j,i) endfor fd(i)=x1/x2 if (i eq 0) then h1=x3/x2 endfor ; print, 'h =', h1 ;--Determine coef. of J for tridiagonal matrix------------------ for i=1,layer-4 do begin aj(i)=fd(i-1)/(dtau(i)*taum(i)) bj(i)=(fd(i)/dtau(i))*(1/taum(i) + 1/taup(i)) + ep cj(i)=fd(i+1)/(dtau(i)*taup(i)) Lj(i)=ep endfor ;--Set boundary conditions for coef. of J---------------------- aj(0)=0.0 bj(0)=fd(0)/taup(0) + h1 cj(0)=fd(1)/taup(0) Lj(0)=0.0 aj(layer-3)=fd(layer-4)/taum(layer-3) bj(layer-3)=fd(layer-3)/taum(layer-3) cj(layer-3)=0.0 Lj(layer-3)=0.0 ;--Solve the tridiagonal matrix for J, (=jn) and recalculate s----- jn=trisol(-aj,bj,-cj,Lj) for i=0,layer-3 do begin s(i)=(1-ep)*jn(i) + ep printf,15, jn(i) endfor diff=0.0 for i=0,layer-3 do begin x=abs(s(i)-sp(i)) if (x gt diff) then begin diff=x y1=s(i) y2=sp(i) endif sp(i)=s(i) endfor maxdiff=diff/(0.5*(y1+y2)) print, k,u(1,layer-3),jn(layer-3),fd(layer-3),h1,maxdiff print, k,u(1,0),jn(0),fd(0),h1,maxdiff printf, 16, k,u(1,layer-3),jn(layer-3),fd(layer-3),maxdiff endfor endfor endfor ; close, 12 ; close, 13 close, 14 close, 15 close, 16 end