C Program to sort out the slightly complicated issue of binary star mass-loss
C ML1 and ML2 are BC1 for stars 1 and 2 - Assumes input is set up with star 1 data
C first
SUBROUTINE MASSLOSS(I,ML1,ML2, BCHORB, BCHSPIN1, BCHSPIN2)
IMPLICIT REAL*8(A-H, L, M, O-Z)
INTEGER MAXMSH
PARAMETER (MAXMSH = 2000)
COMMON H(60,MAXMSH),DH(60,MAXMSH),EPS,DEL,DH0,NMESH,JIN,IW(200)
COMMON /TRANS / HT(26,MAXMSH,2)
COMMON /AUXIN / ICL, ION, JW, IOP, INUC, IBC, ICN, IML(2),ISGTH,
: IMO, IDIFF
COMMON /SODDS / ALPHA, RML, CMG, CSI, CFE, CT(10), AGE, DT, M1,
: EC, BM, ANG, CM, MTA, MTB, TM(2), T0, M0, TC(2), OS, AC, RCD,
: RMG, RHL, XF, DR, AK1 ,RMT, AK2, ITH, IX, IY, IZ, IB, ISX(45),
: TRB
COMMON /INF / VIN(60)
COMMON /DINF / DVIN(60)
COMMON /OP / ZS, LEDD, MM, DG, GRADT, ETH, RRLF, EGR, RR, Q
COMMON /MASLOS/ AIJ(6,5),baseN
COMMON /CNSTS / CPI, PI4, CLN10, CA, CB, CC, CD, CG, CR(2), CEVB,
& CEN, CPL, CMEVMU, CSECYR, LSUN, MSUN, RSUN, TSUNYR
COMMON /EVMODE/ IMODE
COMMON /INERTI/ VI(2)
COMMON /ANGMOM/ VROT1, VROT2, FMAC, FAM, IRAM, IRS1, IRS2
COMMON /WINDS / WINDML(2), FAKEWIND(2), BE
COMMON /TIDES / MENV(2), RENV(2)
COMMON /ZAMS / TKH(2)
DIMENSION T(2), AM(2), M(2), MT(2), AR(2), R(2), DAR(2), L(2), XH(2),
: XHE(2), XC(2), XO(2), HSPIN(2), DHSPIN(2), ML(2),
: RAT(2), RLF(2), DAM(2), WINDACC(2), MOMINER(2),
: BCHSPIN(2), DMOMINER(2), HSPINDT(2), R2O(2), TF(2), HTF(2), HSTF(2),
: OSPIN(2), ACCLIMIT(2), OSC(2)
CBRT(VX) = DEXP(DLOG(VX)/3.0D0)
PS(VX) = 0.5D0*(VX+DABS(VX))
RLOBE(VX) = 0.49D0*VX*VX/(0.6D0*VX*VX+DLOG(1.0D0+VX))
C Set important variables into convenient pairs and calculate other
C necessary quantities
DO ISTAR = 1,IMODE!2
T(ISTAR) = DEXP(VIN(2+15*(ISTAR - 1)))
AM(ISTAR) = VIN(4+15*(ISTAR - 1))
DAM(ISTAR) = DVIN(4+15*(ISTAR - 1))
M(ISTAR) = DEXP(AM(ISTAR))
MT(ISTAR) = (M(ISTAR) - DEXP(AM(ISTAR) - DAM(ISTAR)))/DT
AR(ISTAR) = VIN(7+15*(ISTAR - 1))
R(ISTAR) = DEXP(AR(ISTAR))
DAR(ISTAR) = DVIN(7+15*(ISTAR - 1))
L(ISTAR) = VIN(8+15*(ISTAR - 1))
XH(ISTAR) = VIN(5+15*(ISTAR - 1))
XHE(ISTAR) = VIN(9+15*(ISTAR - 1))
XC(ISTAR) = VIN(10+15*(ISTAR - 1))
XO(ISTAR) = VIN(3+15*(ISTAR - 1))
HSPIN(ISTAR) = VIN(14+15*(ISTAR - 1))
DHSPIN(ISTAR) = DVIN(14+15*(ISTAR - 1))
END DO
HORB = VIN(13)
DHORB = DVIN(13)
IF (IMODE.EQ.2) THEN
BM = M(1) + M(2)
ELSE
M(2) = BM - M(1)
END IF
SEP = (M(1)+M(2))*(HORB/(M(1)*M(2)))**2.0
C Orbital angular velocity
OORB = HORB*BM/(M(1)*M(2)*SEP**2.0)
C surface boundary conditions
DO ISTAR = 1,IMODE!2
IF (ISTAR.EQ.1) ISTAROTHER = 2
IF (ISTAR.EQ.2) ISTAROTHER = 1
RAT(ISTAR) = M(ISTAR)/M(ISTAROTHER)
RLF(ISTAR) = AR(ISTAR) - DLOG(SEP*RLOBE(CBRT(RAT(ISTAR))))
C FUDGE TO AVOID PRE-MS RLOF
IF (AGE.LT.1d3) THEN
RLF(ISTAR) = -1d-1
END IF
C Save RLF data for use in funcs2
HT(24, 1, ISTAR) = RLF(ISTAR)
C Different surface mass bc for *1 or *2 of binary
IF (IB.EQ.1) THEN
C Select mass loss - RJS 24/6/03
IF (IML(ISTAR).EQ.0) BC1 = 0d0
IF (IML(ISTAR).EQ.1) THEN
BC1 = RML*L(ISTAR)*R(ISTAR)/M(ISTAR)
END IF
IF (IML(ISTAR).EQ.2) THEN
BC1 = RML*L(ISTAR)*R(ISTAR)/M(ISTAR)
C Convert Reimers to Blocker mass loss - need to sort out how to do M_ZAMS = 1.5
BC1 = 4.83d-9*(1.5d0)**(-2.1)*(L(ISTAR)/LSUN)**2.7*BC1
END IF
IF (IML(ISTAR).EQ.3) THEN
PERIOD = -2.07 + 1.94*DLOG10(1.4577*R(ISTAR)) - 0.9*DLOG10(M(ISTAR)/2.0)
PERIOD = 10**PERIOD
BC1 = 10**(-11.4 + 0.0123*PERIOD)
BC1 = MSUN/CSECYR * BC1
WIND = BC1
C Superwind, limited to vexp = 15 kms^-1 for P=>500 days
C VEXP = DMAX1(15d0,-13.5d0+5.6d-2*PERIOD)
VEXP = DMAX1(3.0d0,DMIN1(15d0,-13.5d0+5.6d-2*PERIOD))
SWIND = L(ISTAR)*1d26/(CC*1d-2*VEXP*1d3)
C SWIND = 2d-12*CSECYR/(MSUN*1d30)*SWIND
SWIND = SWIND/(MSUN*1d30)
BC1 = DMIN1(BC1,SWIND)
END IF
IF (IML(ISTAR).EQ.4) THEN
C Different surface mass bc for *1 or *2 of binary
c Stuff taken from thesis of L.Dray 2003
COHe=(XC(ISTAR)/3.0+XO(ISTAR)/4.0)/XHE(ISTAR)
SURFXH=XH(ISTAR)
SURFXHe=XHE(ISTAR)/4.0
c if(XHE.lt.0.1d0) THEN
c if(PME.lt.3d18) THEN
c CP3 = CT(9)*3d18
c endif
c endif
c*****NL***************************************
cc c pre-WR: de Jager 1998
zml1=(log10(T(ISTAR))-4.05d0)/0.75d0
zml2=(log10(L(ISTAR)/LSUN)-4.6d0)/2.1d0
BC1=0d0
do n2=0,5
do i2=0,n2
j2=n2-i2
C RJS 16/12/05 - prevent array out of bounds
IF (J2+1.LE.5) THEN
BC1=BC1-AIJ(i2+1,j2+1)*dcos(i2*dacos(zml1))
: *dcos(j2*dacos(zml2))
END IF
enddo
enddo
BC1=(SQRT(ZS*50d0))*(10d0**BC1)/CSECYR
if(zml1.ge.1d0.or.zml2.ge.1d0) BC1=0d0 !-5d90
if(zml1.le.-1d0.or.zml2.le.-1d0) BC1=0d0 !-5d90
BC1 = 2.0*BC1*MSUN
C RJS 9/6/06 - Added some bits to make smooth transitions between WR mass loss
BCPREWR = BC1
c !!!!WOLF RAYET MASS-LOSS!!!!!
CWhen XH(surface)<0.4 and log T > 4.0 the star is in the WNL phase:
CUse a constant rate of: 8e-5 M(sun) yr^-1
BCWNL = 8d-5*MSUN/CSECYR
IF (SURFXH.LT.0.4.AND.log10(T(ISTAR)).GT.3.9) THEN
BC1 = 1d1*(BCWNL - BCPREWR)*(log10(T(ISTAR)) - 3.9) + BCPREWR
END IF
IF (SURFXH.LT.0.4.AND.log10(T(ISTAR)).GT.4.0) BC1 = 8d-5*MSUN/CSECYR
BCWC = 1d-7*(M(ISTAR)/MSUN)**2.5*MSUN/CSECYR
IF (SURFXH.LT.3d-3.AND.log10(T(ISTAR)).GT.4.0) THEN
BC1 = 5d2*(BCWNL - BCWC)*(SURFXH - 1d-3) + BCWC
END IF
CWhen XH(surface)<1e-3 and log T > 4.0 the star is in the WNE,WC or WO
Cphase:
CFor WNE use: 1.0e-7 (M(WR)/M(sun))**2.5 M(sun) yr^-1
IF (SURFXH.LT.1d-3.AND.log10(T(ISTAR)).GT.4.0) THEN
BC1 = 1d-7*(M(ISTAR)/MSUN)**2.5*MSUN/CSECYR
IF (COHe.GT.3d-2) BC1 = 0.6d-7*(M(ISTAR)/MSUN)**2.5*MSUN/CSECYR
END IF
CFor WC/WO use: 0.6e-7 (M(WR)/M(sun))**2.5 M(sun) yr^-1
CStars become WC when (C+O)/He > 3e-2
END IF
IF (IML(ISTAR).EQ.5) THEN
C Newerest massive star mass loss rates from JJE
C Different surface mass bc for *1 or *2 of binary
c Stuff taken from thesis of L.Dray 2003
COHe=(XC(ISTAR)/3.0+XO(ISTAR)/4.0)/XHE(ISTAR)
SURFXH=XH(ISTAR)
SURFXHe=XHE(ISTAR)/4.0
cc c pre-WR: de Jager 1988
zml1=(log10(T(ISTAR))-4.05d0)/0.75d0
zml2=(log10(L(ISTAR)/LSUN)-4.6d0)/2.1d0
BC1=0d0
do n2=0,5
do i2=0,n2
j2=n2-i2
C prevent array out of bounds
IF (J2+1.LE.5) THEN
BC1=BC1-AIJ(i2+1,j2+1)*dcos(i2*dacos(zml1))
: *dcos(j2*dacos(zml2))
END IF
enddo
enddo
BCPREWR1 = BC1
BC1=(SQRT(ZS*50d0))*(10d0**BC1)*MSUN/CSECYR
if(zml1.ge.1d0.or.zml2.ge.1d0) BC1=0d0 !-5d90
if(zml1.le.-1d0.or.zml2.le.-1d0) BC1=0d0 !-5d90
C No factor of 2 this time!
C Vink et al. 2001 rates for OB stars - taken from JJE's code
RVIN=(ZS*50d0)**0.13d0
C RMVA is for hot side of stability jump 27500 - 50000 K
RMVA=-6.697+2.194*LOG10(L(ISTAR)/(LSUN*1d5))-1.313*LOG10(M(ISTAR)/(MSUN*30d0))
: -1.226*LOG10(RVIN*2.6/2.0)+0.933*LOG10(T(ISTAR)/4d4)-10.92*(LOG10(T(ISTAR)/4d4))**2
: +0.85*LOG10(50d0*ZS)
C RMVB is for cool side of jump 12500 - 22500 K
RMVB=-6.688+2.210*LOG10(L(ISTAR)/(LSUN*1d5))-1.339*LOG10(M(ISTAR)/(MSUN*30d0))
: -1.601*LOG10(RVIN*1.3/2.0)+1.07*LOG10(T(ISTAR)/2d4)+0.85*LOG10(50d0*ZS)
C Am I missing an MSUN out of all this? I think so...
IF(T(ISTAR).LE.5d4.AND.T(ISTAR).GT.2.75d4) BC1 =(10**RMVA)*MSUN/CSECYR
IF(T(ISTAR).LE.2.25d4.AND.T(ISTAR).GE.1.25d4) BC1 =(10**RMVB)*MSUN/CSECYR
IF(T(ISTAR).GT.2.25d4.AND.T(ISTAR).LE.2.75d4) THEN
BC1=((T(ISTAR)-2.25d4)*(10**RMVA)+(2.75d4-T(ISTAR))*(10**RMVB))*MSUN/(CSECYR*5d3)
ENDIF
IF(T(ISTAR).LT.1.25d4.AND.T(ISTAR).GT.1d4) THEN
BC1=((T(ISTAR)-1d4)*(10**RMVB)+(1.25d4-T(ISTAR))*(10**BCPREWR1)
: *(50d0*ZS)**0.5d0)*MSUN/(CSECYR*2.5d3)
ENDIF
IF(T(ISTAR).GT.5d4.AND.T(ISTAR).LT.6d4) THEN
BC1=((6d4-T(ISTAR))*(10**RMVA)+(T(ISTAR)-5d4)*(10**BCPREWR1)
: *(50d0*ZS)**0.5d0)*MSUN/(CSECYR*1d4)
ENDIF
BCPREWR = BC1
C Note this smooths out the bistability jump... RJS
CWhen XH(surface)<0.4 and log T > 4.0 the star is in the WNL phase:
BCWNL = -13.6 + 1.63*log10(L(ISTAR)/LSUN)+2.22*log10(XHE(ISTAR))
C Now has metallicity scaling
BCWNL = ((ZS/0.02)**0.5d0)*(10.0**BCWNL)*MSUN/CSECYR
IF (SURFXH.LT.0.4.AND.log10(T(ISTAR)).GT.3.9) THEN
BC1 = 1d1*(BCWNL - BCPREWR)*(log10(T(ISTAR)) - 3.9) + BCPREWR
END IF
IF (SURFXH.LT.0.4.AND.log10(T(ISTAR)).GT.4.0) BC1 = BCWNL
CWhen XH(surface)<1e-3 and log T > 4.0 the star is in the WNE,WC or WO
Cphase:
BCWC = -8.3 + 0.84*log10(L(ISTAR)/LSUN) + 2.04*log10(XHE(ISTAR)) +
: 1.04*log10(1d0-XHE(ISTAR))
C Also now metallicity scaled
BCWC = ((ZS/0.02)**0.5d0)*(10.0**BCWC)*MSUN/CSECYR
IF (SURFXH.LT.1d-3.AND.log10(T(ISTAR)).GT.4.0) THEN
C WC rate
IF (COHe.GT.2d-2) THEN
BC1 = 1d2*(BCWC-BCWNL)*(COHe - 2d-2) + BCWNL
END IF
IF (COHe.GT.3d-2) BC1 = BCWC
C IF (COHe.GT.1d0) BC1 = 1.9d-5*MSUN/CSECYR
END IF
END IF
C Save wind mass loss for orbital angular momentum calculation
C Enhance wind mass loss by 1/(1-Omega/Omega_crit)
OSPIN(ISTAR) = HSPIN(ISTAR)/VI(ISTAR)*SQRT(CG)
OCRIT = DSQRT(6.67d-11*(M(ISTAR)*1d30)/((1d9*R(ISTAR))**3.0))
OSC(ISTAR) = OSPIN(ISTAR)/OCRIT
BC1 = BC1/DMAX1(1d-2,(1-OSC(ISTAR)/0.8))
WINDML(ISTAR) = BC1
END IF
ML(ISTAR) = BC1
END DO
C Now decide on what the BC really is - needs to be outside the loop to know about
C both Roche Lobes
C Messed up Common Envelope stuff
FAKEWIND(1) = 0d0
FAKEWIND(2) = 0d0
MASSLIMIT = 1d-2 !2.5d-4
C Set limit for mass accretion at M/kelvin-helmholtz timescale
DO ISTAR = 1,IMODE
ACCLIMIT(ISTAR) = M(ISTAR)/TKH(ISTAR)*DMAX1(0d0,(1-OSC(ISTAR)/0.8))
END DO
C Fakewind deals with mass-loss in CE systems - will interfer with normal evolution
DO ISTAR=1,IMODE
IF (ISTAR.EQ.1) ISTAROTHER = 2
IF (ISTAR.EQ.2) ISTAROTHER = 1
ML(ISTAR) = MT(ISTAR) - RMG*M(ISTAR) + ML(ISTAR)
: + DMIN1(RMT*(PS(RLF(ISTAR)))**3,MASSLIMIT*MSUN/CSECYR)
IF (IMODE.EQ.2) THEN
C Add (1-omega/omega_crit) to reduce accretion rate
ML(ISTAR) = ML(ISTAR) - DMIN1(FMAC*DMIN1(RMT*(PS(RLF(ISTAROTHER)))**3.0
: ,MASSLIMIT*MSUN/CSECYR),ACCLIMIT(ISTAR)*MSUN/CSECYR)
END IF
END DO
FAKEWIND(1) = DMAX1(0d0,DMIN1(RMT*(PS(RLF(1)))**3,MASSLIMIT*MSUN/CSECYR) -
C : DMIN1(FMAC*RMT*(PS(RLF(1)))**3.0,ACCLIMIT(1)*MSUN/CSECYR))
: DMIN1(FMAC*RMT*(PS(RLF(1)))**3.0,ACCLIMIT(2)*MSUN/CSECYR))
C Should ACCLIMIT be the other star? I think so...
C Fakewind is supposed to be total mass lost by star one minus mass accreted by 2
FAKEWIND(2) = DMAX1(0d0,DMIN1(RMT*(PS(RLF(2)))**3,MASSLIMIT*MSUN/CSECYR) -
C : DMIN1(FMAC*RMT*(PS(RLF(2)))**3.0,ACCLIMIT(2)*MSUN/CSECYR))
: DMIN1(FMAC*RMT*(PS(RLF(2)))**3.0,ACCLIMIT(1)*MSUN/CSECYR))
C Ignore wind accretion if RLOF occurs
WINDACC(1) = 0d0
WINDACC(2) = 0d0
IF ((PS(RLF(1)))**3.0.EQ.0d0.AND.(PS(RLF(2))).EQ.0d0) THEN
C This assumes you're only changing the mass by winds with no RLF
C Wind accretion from Hurley et al. 2002
C Accrete material from the greater mass loser
IF (WINDML(1).GT.WINDML(2)) THEN
IDONOR = 1
IACC = 2
ELSE
IDONOR = 2
IACC = 1
END IF
VORB2 = CG*BM/SEP
C BetaW depends on spectral type = 7 for O-type stars
BETAW = 7.0
VESC2 = 2.0*BETAW*CG*M(IDONOR)/R(IDONOR)
V2 = VORB2/VESC2
C Have assumed a circular orbit - with alpha_w = 3/2
MTACC = (CG*M(IACC)/VESC2)**2.0*(3.0/(4.0*SEP**2.0))/(1+V2)**(3.0/2.0)
MTACC = MTACC*WINDML(IDONOR)
MTACC = DMIN1(MTACC,0.8*WINDML(IDONOR))
WINDACC(IACC) = MTACC
ML(IACC) = ML(IACC) - MTACC
END IF
C Store WINDML + WINDACC for use in determining when variable composition accretion
C is needed - only if sum of one of them is -ve. Only do this if there's no RLOF
IF ((PS(RLF(1)))**3.0.EQ.0d0.AND.(PS(RLF(2))).EQ.0d0) THEN
HT(23,1,1) = 0d0 !WINDML(1) - WINDACC(1)
HT(23,1,2) = 0d0 !WINDML(2) - WINDACC(2)
ELSE
HT(23,1,1) = 0d0
HT(23,1,2) = 0d0
END IF
C Fudge CE stuff
C IF(RLF(1).GT.0d0) THEN
C ML(1) = 1d-2*MSUN/CSECYR + ML(1)
C WINDML(1) = 1d-2*MSUN/CSECYR
C IF (M(2)/MSUN.LT.17) THEN
C ML(2) = -5d-3*MSUN/CSECYR + ML(2)
C END IF
C END IF
ML1 = ML(1)
ML2 = ML(2)
C Tidal friction stuff. From Hut (1981), using Zahn (1977) for frictional timescale
DO ISTAR = 1,IMODE
C TF in years
C TF(ISTAR) = 3.5*(M(ISTAR)/MSUN*(R(ISTAR)/RSUN)**2.0
C : /(L(ISTAR)/LSUN))**1.0/3.0*CSECYR
C DkT = 0.1/(R(ISTAR)**3.0/(CG*M(ISTAR)*TF(ISTAR)))
QQ = 1/RAT(ISTAR)
FQ = QQ*(1+QQ)
RADSEP = R(ISTAR)/SEP
C RJS 14/1/08
C Alternative tidal prescription from Hurley, Trout & Pols (2002)
C Convective envelope stars
C Convective turnover time
ME = DMAX1(0d0,M(ISTAR) - MENV(ISTAR))
RE = DMAX1(0d0,R(ISTAR) - RENV(ISTAR))
TCONV = ME*RE*(R(ISTAR) - 0.5*RE)
TCONV = 0.4311*(TCONV/(3d0*L(ISTAR)))**(1d0/3d0)*CSECYR
TCONV = DMAX1(TCONV,1d0)
C Tidal pumping timescale
PID = 1d0/DABS(OORB - OSPIN(ISTAR))/SQRT(CG)
C Numerical factor supressing tide
FCONV = DMIN1(1d0,(0.5d0*PID/TCONV)**2d0) ! check units match, should be ok now...
C Form of k/T
DkT = (2d0/21d0)*(FCONV/TCONV)*(ME/M(ISTAR))
C Dynamical tide with radiative damping, also from HTP 2002
C Don't understand their formula, but from eq 28 & 41
DkTR = SQRT(CG*M(ISTAR)/R(ISTAR)**3d0)*(1d0+QQ)**(5d0/6d0)
C This is a fit to tabulated data for MS stars from Zahn '75 -- how valid is it
C for different metallicities or non-MS stars?
E2 = 1.592d-9*(M(ISTAR)/MSUN)**2.84d0
DkTR = DkTR*E2*RADSEP**(5d0/2d0)
C write (*,*) DkT, DkTR
DkT = DkT + DkTR
C HTF is rate of change of orbital angular momentum
C Neglecting eccentricity - for now...
OSPIN(ISTAR) = HSPIN(ISTAR)/VI(ISTAR)
HTF(ISTAR) = -3.0*DKT*FQ*RADSEP**8.0*HORB*(1 - (OSPIN(ISTAR)/OORB))
C Spin angular momentum - HSTF.
HSTF(ISTAR) = 3.0*DKT*QQ**2.0*RADSEP**6.0*M(ISTAR)*R(ISTAR)**2.0
: *OORB*(1 - (OSPIN(ISTAR)/OORB))
END DO
IF (AGE.LT.1d3) THEN
HTF(1) = 0d0
HTF(2) = 0d0
HSTF(1) = 0d0
HSTF(2) = 0d0
END IF
C IF (RLF(1).GT.0d0) THEN
C IF (R(1)/SEP.GT.0.45) THEN
C FACAM = DMAX1(0d0, (0.55 - (R(1)/SEP))/0.1)
C HTF(1) = 0d0
C HTF(2) = 0d0
C HSTF(1) = 0d0
C HSTF(2) = 0d0
C END IF
C Boundary condition for orbit
C Evolution of orbital Ang. Mom. from Hurley et al (2002)
IF (IMODE.EQ.2) THEN
BCHORB = DHORB/DT + (HORB/BM)*((WINDML(1)+FAKEWIND(1))*M(2)/M(1) - WINDACC(1)
C : + (WINDML(2)+FAKEWIND(2))*M(2)/M(1) - WINDACC(2)) - HTF(1) - HTF(2)
: + (WINDML(2)+FAKEWIND(2))*M(1)/M(2) - WINDACC(2)) - HTF(1) - HTF(2)
ELSE
BCHORB = DHORB/DT + HORB/BM*((WINDML(1)+FAKEWIND(1))*M(2)/M(1)) - HTF(1)
END IF
C IF (RLF(1).GT.0d0) THEN
C BCHORB = DHORB/DT + 1.5*HORB*WINDML(1)/(M(1)+M(2))
C END IF
C Boundary condition for spin period
C Assuming solid body rotation
R2O(1) = HSPIN(1)/VI(1)*R(1)**2.0
R2O(2) = HSPIN(2)/VI(2)*R(2)**2.0
DO ISTAR = 1,IMODE
IF (ISTAR.EQ.1) ISTAROTHER = 2
IF (ISTAR.EQ.2) ISTAROTHER = 1
C Assume mass lost/gained as a shell - hence factor of 2/3. Assume mu_w = 1
C Loss from star
C This has to be made consistent with mass loss & accretion limits above
C HSPINDT(ISTAR) = 2.0/3.0*(WINDML(ISTAR) + RMT*(PS(RLF(ISTAR)))**3)*R2O(ISTAR)
HSPINDT(ISTAR) = 2.0/3.0*(WINDML(ISTAR)/DMAX1(1d-2,(1-OSC(ISTAR)/0.8))
: + DMIN1(RMT*(PS(RLF(ISTAR)))**3
: ,MASSLIMIT*MSUN/CSECYR))*R2O(ISTAR)
IF (IMODE.EQ.2) THEN
HSPINDT(ISTAR) = HSPINDT(ISTAR) - 2.0/3.0*(WINDACC(ISTAROTHER) +
: FAM*DMIN1(FMAC*DMIN1(RMT*(PS(RLF(ISTAROTHER)))**3.0
: ,MASSLIMIT*MSUN/CSECYR),ACCLIMIT(ISTAR)*MSUN/CSECYR))*R2O(ISTAROTHER)
END IF
BCHSPIN(ISTAR) = DHSPIN(ISTAR)/DT + HSPINDT(ISTAR) - HSTF(ISTAR)
END DO
BCHSPIN1 = BCHSPIN(1)
BCHSPIN2 = BCHSPIN(2)
RETURN
END