! Program for averaging the dust phase functions ! over the stellar disk with quadratic limb darkening. ! ! authors: J. Budaj ! contact: budaj@ta3.sk ! published: ! non original routins: ! ! input: ! diskaver.in (unit 2): ! u1,u2 (quad. limb dark coef), rstar[Rsol], dist[Rsol] ! diskaver.dat (unit 2, pf data table) ! diskaver_angles.in (unit 8, angles for output, optional) ! output: diskaver.out (unit 3) ! program diskaver implicit none ! IMPLICIT DOUBLE PRECISION (A-H,O-Z) integer, parameter:: i4=4,i8=8 integer, parameter:: mdim=2000,msize=21,mwave=400,mang1=65 real(i8), parameter:: clight=2.99792458d10 real(i8):: xp(mdim),yp(mdim),area(mdim,mdim) real(i8):: size(msize),wave(mwave),sumip(mang1) real(i8):: ang1(mang1),pf1(msize,mwave,mang1),pf(mang1) real(i8):: ang3(mang1),pf3(msize,mwave,mang1) real(i8):: pff(mang1),pffn(mang1) real(i8):: dsumi(mdim,mdim),r(mdim,mdim),wp(mdim,mdim) real(i8):: rad2deg,deg2rad,asize,dsize,awave real(i8):: rsol,pi,u1,u2,rstar,dist,rstar2,dist2 real(i8):: deltam,deltad,solang,resol,rm,rm2 real(i8):: step,stepx,stepy,r2,disksurf,zp,zp2,wp2 real(i8):: a,b,c,z,w,w2,domega,omega,ctheta,fld,sumi real(i8):: alpha0,calpha,alphad,pfxy,dsumip,gmean integer:: nxy,nbodx,nbody,nsize,nwave,nang1,i,j,k,ndat3 integer:: isize,iwave,kk WRITE(*,*) WRITE(*,*)'-------------------------------------------' WRITE(*,*)' Averaging phase functions over stellar disk' WRITE(*,*)'-------------------------------------------' WRITE(*,*) WRITE(*,*)' INPUT FILES: ' WRITE(*,*)' diskaver.in' WRITE(*,*)' diskaver.dat' WRITE(*,*)' diskaver_angles.in' WRITE(*,*)' OUTPUT FILE: diskaver.out' WRITE(*,*) rsol=6.9599d10 ! pi=4.d0*datan(1.d0) pi=dacos(-1.d0) deg2rad=pi/180.d0 rad2deg=180.d0/pi open(2,file='diskaver.in',status='old') read(2,*) read(2,*)u1,u2,rstar,dist,asize,awave close(2) rstar=rstar*rsol rstar2=rstar*rstar dist=dist*rsol dist2=dist*dist ! read original phase functions open(2,file='diskaver.dat',status='old') nsize=21 nwave=400 ! nsize=1 ! nwave=1 nang1=65 do i=1,nsize read(2,*) read(2,*)size(i) do j=1,nwave read(2,*)wave(j) do k=1,nang1 read(2,*)ang1(k),pf1(i,j,k) enddo enddo enddo close(2) dsize=dlog10(asize) call locate(size,nsize,dsize,i) if (i>=nsize)then isize=nsize elseif(dabs(dsize-size(i)) <= dabs(dsize-size(i+1)))then isize=i else isize=i+1 endif write(*,*)'log10size[mic]=',size(isize) awave=clight/(awave*1.d-4) call locate(wave,nwave,awave,i) if (i>=nwave)then iwave=nwave elseif(dabs(awave-wave(i)) <= dabs(awave-wave(i+1)))then iwave=i else iwave=i+1 endif write(*,*)'wave[mic]=',clight/wave(iwave)*1.d4,iwave write(*,'(a,es12.4)')'freq[HZ]=',wave(iwave) ! set up new angle values for output open(8,file='diskaver_angles.in',status='old') i=1 23 read(8,*,end=24)ang3(i) i=i+1 goto 23 24 ndat3=i-1 write(*,*)'Read No of angles from diskaver_angles.in',ndat3 if(ndat3.gt.mang1) stop'ndat3>mang1' close(8) open(3,file='diskaver.out',status='unknown') ! set the grid points x,y over the disk deltam=dasin(rstar/dist) deltad=rad2deg*deltam write(*,*)'angular radius of the star [deg]',deltad ! solid angle solang=2.d0*pi*(1.d0-dcos(deltam)) ! resolution 0.05 degrees resol=0.05d0 write(*,*)'resolution [deg]',resol if(resol>deltad)then stop'no need of calculations, star is already a point source' endif nxy=dint(deltad/resol+1.d0) ! rm -maximal radius visible rm2=rstar**2*(1.d0-(rstar/dist)**2) rm=dsqrt(rm2) step=rm/dble(nxy) nbodx=nint(2.d0*rm/step)+1 nbody=nint(rm/step)+1 ! make the number of points odd if(nbodx/2.eq.(nbodx+1)/2)nbodx=nbodx-1 write(*,*)'nbodx,nbody=',nbodx,nbody stepx=2.d0*rm/dble(nbodx-1) stepy=rm/dble(nbody-1) do i=1,nbodx xp(i)=-rm+stepx*dble(i-1) enddo do i=1,nbody yp(i)=0.d0+stepy*dble(i-1) enddo ! calculates area(i,j) of elements call surf(mdim,mdim,nbodx,nbody,xp,yp,area) disksurf=0.d0 omega=0.d0 sumi=0.d0 zp=dist-rstar2/dist zp2=zp*zp do i=1,nbodx do j=1,nbody dsumi(i,j)=0.d0 r2=xp(i)**2+yp(j)**2 r(i,j)=dsqrt(r2) if(r(i,j).le.rm)then disksurf=disksurf+area(i,j) wp2=r2+zp2 wp(i,j)=dsqrt(wp2) a=1.d0+r2/zp2 b=-2.d0*dist c=dist2-rstar2 z=(-b-dsqrt(b*b-4.d0*a*c))/2.d0/a w=z/zp*wp(i,j) w2=w*w domega=area(i,j)*zp/wp(i,j)/wp2 omega=omega+domega ctheta=(dist2-rstar2-w2)/2.d0/rstar/w fld=1.d0-u1*(1.d0-ctheta)-u2*(1.d0-ctheta)**2 dsumi(i,j)=fld*domega sumi=sumi+dsumi(i,j) endif enddo enddo write(*,*)'observed angle',omega/solang write(*,*)'observed disk',disksurf/pi/rm2 write(*,*)'disk averaged norm Intensity',sumi/omega ! uncomment the next 4 lines in case you want to calculate ! whole size x lambda grid ! do isize=1,nsize ! write(3,*) ! write(3,10)size(isize) ! do iwave=1,nwave do k=1,ndat3 alpha0=deg2rad*ang3(k) sumip(k)=0.d0 do i=1,nbodx do j=1,nbody if(r(i,j).le.rm)then calpha=(-xp(i)*dsin(alpha0)+zp*dcos(alpha0))/wp(i,j) alphad=rad2deg*dacos(calpha) ! call locate(ang1,nang1,alphad,kk) call intrp2(ang1,pf1,isize,iwave,nang1,mang1,msize, & mwave,alphad,pfxy) dsumip=dsumi(i,j)*pfxy sumip(k)=sumip(k)+dsumip endif enddo enddo pf3(isize,iwave,k)=sumip(k)/sumi enddo ! renormalize and recalculate mean cosine do k=1,ndat3 pff(k)=pf3(isize,iwave,k) enddo call pfnorm(ang3,pff,pffn,ndat3,mang1,gmean) write(3,20)wave(iwave),gmean write(*,20)wave(iwave),gmean do k=1,ndat3 write(3,40)ang3(k),pff(k),pffn(k) enddo ! uncomment the next 2 lines in case you want to calculate ! whole size x lambda grid ! enddo ! enddo 10 format(f5.1) 20 format(es13.5,f9.5) 40 format(f8.2,2es12.4) close(3) stop end !----------------------------------------------------------------------- subroutine surf(ndim1,ndim2,nbod1,nbod2,ar,at,area) ! to calculate area(i,j)-area of a projected surface element ! input: ar,at ! output: area implicit double precision (a-h,o-z) dimension ar(ndim1),at(ndim2),area(ndim1,ndim2) pi=3.1415926535897931d0 do i=2,nbod1-1 do j=2,nbod2-1 area(i,j)=(ar(i+1)-ar(i-1))*(at(j+1)-at(j-1))*0.25d0 enddo enddo pomi12=(ar(2)-ar(1))*0.5d0 pominn=(ar(nbod1)-ar(nbod1-1))*0.5d0 do j=2,nbod2-1 area(1,j)=pomi12*(at(j+1)-at(j-1))*0.5d0 area(nbod1,j)=pominn*(at(j+1)-at(j-1))*0.5d0 enddo do i=2,nbod1-1 dar=ar(i+1)-ar(i-1) area(i,1)=dar*(at(2)-at(1))*0.25d0 area(i,nbod2)=dar*(at(nbod2)-at(nbod2-1))*0.25d0 enddo area(1,1)=pomi12*(at(2)-at(1))*0.5d0 area(1,nbod2)=pomi12*(at(nbod2)-at(nbod2-1))*0.5d0 area(nbod1,1)=pominn*(at(2)-at(1))*0.5d0 area(nbod1,nbod2)=pominn*(at(nbod2)-at(nbod2-1))*0.5d0 return end !----------------------------------------------------------------------- subroutine intrp(x,y,ndat,mdat,xvalue,yvalue) ! general subroutine for linear interpolation ! x may be increasing or decreasing ! input: x, y, xvalue ! output: yvalue implicit double precision (a-h,o-z) dimension x(mdat),y(mdat) if(x(1).lt.x(ndat).and.xvalue.le.x(1))then yvalue=y(1) elseif(x(1).lt.x(ndat).and.xvalue.ge.x(ndat))then yvalue=y(ndat) elseif(x(1).gt.x(ndat).and.xvalue.ge.x(1))then yvalue=y(1) elseif(x(1).gt.x(ndat).and.xvalue.le.x(ndat))then yvalue=y(ndat) else call locate(x,ndat,xvalue,jj) der=(y(jj+1)-y(jj))/(x(jj+1)-x(jj)) yvalue=der*(xvalue-x(jj))+y(jj) endif return end !----------------------------------------------------------------------- subroutine intrp2(x,y,isize,iwave,ndat,mdat,msize,mwave, & xvalue,yvalue) ! general subroutine for linear interpolation ! x may be increasing or decreasing ! input: x, y, xvalue ! output: yvalue implicit double precision (a-h,o-z) dimension x(mdat),y(msize,mwave,mdat) if(x(1).lt.x(ndat).and.xvalue.le.x(1))then yvalue=y(isize,iwave,1) elseif(x(1).lt.x(ndat).and.xvalue.ge.x(ndat))then yvalue=y(isize,iwave,ndat) elseif(x(1).gt.x(ndat).and.xvalue.ge.x(1))then yvalue=y(isize,iwave,1) elseif(x(1).gt.x(ndat).and.xvalue.le.x(ndat))then yvalue=y(isize,iwave,ndat) else call locate(x,ndat,xvalue,jj) der=(y(isize,iwave,jj+1)-y(isize,iwave,jj))/(x(jj+1)-x(jj)) yvalue=der*(xvalue-x(jj))+y(isize,iwave,jj) endif return end !----------------------------------------------------------------------- SUBROUTINE LOCATE(XX,N,X,J) ! taken from Numerical recipes, ! search an ordered table by bisection ! given array xx and value x returns j such that ! x is between xx(j),xx(j+1) IMPLICIT DOUBLE PRECISION (A-H,O-Z) ! include 'param.inc' ! DIMENSION XX(mstarx) DIMENSION XX(N) JL=0 JU=N+1 10 IF(JU-JL.GT.1)THEN JM=(JU+JL)/2 IF((XX(N).GT.XX(1)).EQV.(X.GT.XX(JM)))THEN JL=JM ELSE JU=JM ENDIF GO TO 10 ENDIF J=JL RETURN END !----------------------------------------------------------------------- subroutine pfnorm(angf2,pff,pffn,nan,mdat,gmean) implicit double precision (a-h,o-z) dimension pff(mdat),pffn(mdat),ta(mdat),angf2(mdat) ! pi=3.1415926535897931d0 pi=dacos(-1.d0) do it=1,nan ta(it)=pi/180.d0*angf2(it) enddo sum=0.d0 gmean=0.d0 do it=2,nan sump1=pff(it-1)*dsin(ta(it-1)) gsum1=sump1*dcos(ta(it-1)) sump2=pff(it)*dsin(ta(it)) gsum2=sump2*dcos(ta(it)) dt=ta(it)-ta(it-1) sump=(sump1+sump2)*dt*0.5d0 sum=sum+sump gsum=(gsum1+gsum2)*dt*0.5d0 gmean=gmean+gsum enddo gmean=gmean/sum do it=1,nan pffn(it)=pff(it)*2.d0/sum enddo ! write(*,'(a,es13.5,f8.4,2i4)')'int pf sint dt, g=',sum,gmean return end