[mesa-users] Error in Vink wind formula

[Bill Paxton]

Hi Roni,

Okay -- so I've been using the Galactic values for vinf/vesc instead of correcting that ratio for Z.
I've made the correction, so it will show up in the next release.   Thanks!!!

Cheers,
Bill

p.s. if anyone wants to patch their current version of mesa for this bug,
here's the new eval_vink_wind routine that lives in star/private/winds.f
replace your version of that routine with the following, then
cd mesa/star/test
./mk
./ck
./export


         subroutine eval_vink_wind(w)
            real(dp), intent(inout) :: w
            real(dp) :: alfa, w1, w2, Teff_jump, logMdot, dT, vinf_div_vesc

            ! alfa = 1 for hot side, = 0 for cool side
            if (T1 > 27500d0) then
               alfa = 1
            else if (T1 < 22500d0) then
               alfa = 0
            else ! use Vink et al 2001, eqns 14 and 15 to set "jump" temperature
               Teff_jump = 1d3*(61.2d0 + 2.59d0*(-13.636d0 + 0.889d0*log10(Z/Zsolar)))
               dT = 100d0
               if (T1 > Teff_jump + dT) then
                  alfa = 1
               else if (T1 < Teff_jump - dT) then
                  alfa = 0
               else
                  alfa = (T1 - (Teff_jump - dT)) / (2*dT)
               end if
            end if
            
            if (alfa > 0) then ! eval hot side wind (eqn 24)
               vinf_div_vesc = 2.6d0 ! this is the hot side galactic value
               vinf_div_vesc = vinf_div_vesc*(Z/Zsolar)**0.13d0 ! corrected for Z
               logMdot = &
                  - 6.697d0 &
                  + 2.194d0*log10(L1/Lsun/1d5) &
                  - 1.313d0*log10(M1/Msun/30) &
                  - 1.226d0*log10(vinf_div_vesc/2d0) &
                  + 0.933d0*log10(T1/4d4) &
                  - 10.92d0*log10(T1/4d4)**2 &
                  + 0.85d0*log10(Z/Zsolar)
               w1 = 10**logMdot
            else
               w1 = 0
            end if
            
            if (alfa < 1) then ! eval cool side wind (eqn 25)
               vinf_div_vesc = 1.3d0 ! this is the cool side galactic value
               vinf_div_vesc = vinf_div_vesc*(Z/Zsolar)**0.13d0 ! corrected for Z
               logMdot = &
                  - 6.688d0 &
                  + 2.210d0*log10(L1/Lsun/1d5) &
                  - 1.339d0*log10(M1/Msun/30) &
                  - 1.601d0*log10(vinf_div_vesc/2d0) &
                  + 1.07d0*log10(T1/2d4) &
                  + 0.85d0*log10(Z/Zsolar)
               w2 = 10**logMdot
            else
               w2 = 0
            end if
            
            w = alfa*w1 + (1 - alfa)*w2
            
            if (dbg) write(*,*) 'vink wind', w
            
         end subroutine eval_vink_wind





On Apr 16, 2012, at 7:53 AM, Roni Waldman wrote:

> Hi Bill,
> In eqns 24, 15 you have the term log(vinf/vesc/2).
> For this the Galactic ratio vinf/vesc is 2.6 or 1.3 for hot and cool side respectively, as in mesa.
> However (see eqns 21, 22 and related paragraphs in the paper) vinf has also a Z dependence, namely, according to Leitherer et al. 1992, vinf ~ Z**0.13.
> This will give a shallower dependence on Z, Mdot ~ Z**0.69 and Z**0.64 for hot and cool stars respectively, rather than Z**0.85 as in mesa.
> To fix this, you need to multiply the vinf/vesc term by (Z/Zsolar)**0.13.
> Roni
>  
> From: Bill Paxton [mailto:paxton at kitp.ucsb.edu] 
> Sent: יום ב 16 אפריל 2012 17:41
> To: Roni Waldman
> Cc: mesa-users at lists.sourceforge.net
> Subject: Re: Error in Vink wind formula
>  
> Hi Roni,
>  
> Q: I checked again, and what is currently in mesa matches the Vink et al eqns 24 and 25,
> so the terms you are providing are a correction to the paper?   Or is there something else
> in the paper that I missed?  Or do I have the wrong paper?

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://lists.mesastar.org/pipermail/mesa-users/attachments/20120416/807cadab/attachment.html>