/******************************************* * * FILE TOTFLUX.C * * This file contains all the functions necessary to calculate * the flux from the binary at a specific orbital phase. * * NOTE: Orbital phase 0.0 is defined to be inferior conjunction * of star2 (= star2 closest to the earth). * Also, internally, the program forces -0.5 < orbphase <= 0.5. * **********************************************/ #include "header.h" double totalFlux( double orbphase ) /************************* * * Calculate the total flux reaching the earth from the binary * at one particular orbital phase. * **************************/ { double x1 [SIZEDIMEN1] [SIZEDIMEN2], y1 [SIZEDIMEN1] [SIZEDIMEN2], mu1 [SIZEDIMEN1] [SIZEDIMEN2], x2 [SIZEDIMEN1] [SIZEDIMEN2], y2 [SIZEDIMEN1] [SIZEDIMEN2], mu2 [SIZEDIMEN1] [SIZEDIMEN2], xdisk [SIZEDIMEN1] [SIZEDIMEN2], ydisk [SIZEDIMEN1] [SIZEDIMEN2], mudisk [SIZEDIMEN1] [SIZEDIMEN2]; double edge1x[LINEPOINTS], edge1y[LINEPOINTS], edge2x[LINEPOINTS], edge2y[LINEPOINTS], outeredgediskx[LINEPOINTS], outeredgedisky[LINEPOINTS], topedgediskx[LINEPOINTS], topedgedisky[LINEPOINTS], inneredgediskx[LINEPOINTS], inneredgedisky[LINEPOINTS]; double inclination, x, y, z, lsky[4], msky[4], nsky[4], originx, coordsign, transmission, flux; long int keepside, itheta, iphi, ir, nedge1, nedge2, max, nouteredgedisk, ntopedgedisk, ninneredgedisk; long int maxthetaindex1, maxphiindex1, maxthetaindex2, maxphiindex2, maxrindex, maxphiindexdisk; double visible1 [SIZEDIMEN1] [SIZEDIMEN2], visible2 [SIZEDIMEN1] [SIZEDIMEN2], visibledisk [SIZEDIMEN1] [SIZEDIMEN2]; /******************************** * * Project all the objects onto the plane of the sky and put the * projected images into the coordinate system centered on star 2. * Then find the edges of all the active objects. Set visible[][] * equal to zero for those points beyond the limbs of the objects. * *********************************/ if( (orbphase < -0.5) || (orbphase >0.5) ) quit("Orbital phase must lie between -0.5 and 0.5 in totFlux."); inclination = starsys.inclination; coordsys( lsky, msky, nsky, inclination, orbphase); nedge1 = 0; nedge2 = 0; ntopedgedisk = 0; ninneredgedisk = 0; nouteredgedisk = 0; if( strcmp( mode.disk, "ON") == 0 ) { maxrindex = disktop.rpoints - 1; maxphiindexdisk = disktop.phipoints - 1; coordsign = -1.0; originx = 1.0; for( ir = 0; ir <= maxrindex; ir++) { for( iphi = 0; iphi <= maxphiindexdisk; iphi++) { x = originx + coordsign * rdisk[ir][iphi] * ldisk[ir][iphi]; y = coordsign * rdisk[ir][iphi] * mdisk[ir][iphi]; z = rdisk[ir][iphi] * ndisk[ir][iphi]; project( &xdisk[ir][iphi], &ydisk[ir][iphi], x, y, z, lsky, msky); mudisk[ir][iphi] = coordsign * normdiskx[ir][iphi] * nsky[1] + coordsign * normdisky[ir][iphi] * nsky[2] + normdiskz[ir][iphi] * nsky[3]; visibledisk[ir][iphi] = 1.0; if( mudisk[ir][iphi] < 0.0 ) visibledisk[ir][iphi] = 0.0; } } nouteredgedisk = diskedges( 1, maxrindex, maxphiindexdisk, xdisk, ydisk, mudisk, outeredgediskx, outeredgedisky); if( (90.0 - starsys.inclination) < disktop.flareangle ) ntopedgedisk = diskedges( 2, maxrindex, maxphiindexdisk, xdisk, ydisk, mudisk, topedgediskx, topedgedisky); if( starsys.inclination > 0.1 ) ninneredgedisk = diskedges( 3, maxrindex, maxphiindexdisk, xdisk, ydisk, mudisk, inneredgediskx, inneredgedisky); } if( strcmp( mode.star1, "ON") == 0 ) { maxthetaindex1 = star1.thetapoints - 1; maxphiindex1 = star1.phipoints - 1; coordsign = -1.0; originx = 1.0; for( itheta = 0; itheta <= maxthetaindex1; itheta++) { for( iphi = 0; iphi <= maxphiindex1; iphi++) { x = originx + coordsign * r1[itheta][iphi] * l1[itheta][iphi]; y = coordsign * r1[itheta][iphi] * m1[itheta][iphi]; z = r1[itheta][iphi] * n1[itheta][iphi]; project( &x1[itheta][iphi], &y1[itheta][iphi], x, y, z, lsky, msky); mu1[itheta][iphi] = coordsign * normalx1[itheta][iphi] * nsky[1] + coordsign * normaly1[itheta][iphi] * nsky[2] + normalz1[itheta][iphi] * nsky[3]; visible1[itheta][iphi] = 1.0; if( mu1[itheta][iphi] < 0.0 ) visible1[itheta][iphi] = 0.0; } } nedge1 = staredge( star1.thetapoints, star1.phipoints, x1, y1, mu1, edge1x, edge1y); cutoffstar1( ninneredgedisk, inneredgediskx, inneredgedisky, nedge1, edge1x, edge1y); cutoffstar1( ntopedgedisk, topedgediskx, topedgedisky, nedge1, edge1x, edge1y); } if( strcmp( mode.star2, "ON") == 0 ) { maxthetaindex2 = star2.thetapoints - 1; maxphiindex2 = star2.phipoints - 1; coordsign = 1.0; originx = 0.0; for( itheta = 0; itheta <= maxthetaindex2; itheta++) { for( iphi = 0; iphi <= maxphiindex2; iphi++) { x = originx + coordsign * r2[itheta][iphi] * l2[itheta][iphi]; y = coordsign * r2[itheta][iphi] * m2[itheta][iphi]; z = r2[itheta][iphi] * n2[itheta][iphi]; project( &x2[itheta][iphi], &y2[itheta][iphi], x, y, z, lsky, msky); mu2[itheta][iphi] = coordsign * normalx2[itheta][iphi] * nsky[1] + coordsign * normaly2[itheta][iphi] * nsky[2] + normalz2[itheta][iphi] * nsky[3]; visible2[itheta][iphi] = 1.0; if( mu2[itheta][iphi] < 0.0 ) visible2[itheta][iphi] = 0.0; } } nedge2 = staredge( star2.thetapoints, star2.phipoints, x2, y2, mu2, edge2x, edge2y); } /**************************** * * Now find which points are occulted by other objects in * the system. * ******************************/ /************************************* * * First occult star 2 by star 1 and the disk. Note that this * version does not handle a disk with a transparent center correctly. * **************************************/ if( strcmp( mode.star2, "ON") == 0 ) { if( fabs( orbphase ) > 0.25 ) { if( strcmp( mode.star1, "ON") == 0 ) { keepside = 1; transmission = 0.0; occult( keepside, transmission, maxthetaindex2, maxphiindex2, x2, y2, visible2, nedge1, edge1x, edge1y); } if( strcmp( mode.disk, "ON") == 0 ) { keepside = 1; transmission = disktop.transmission; occult( keepside, transmission, maxthetaindex2, maxphiindex2, x2, y2, visible2, nouteredgedisk, outeredgediskx, outeredgedisky); } } } /************************** * * Then the disk. * ****************************/ /******************** * * Treat the edge of the disk and the surface of the disk * as if they are separate entities. * The difference between the top of the disk and the edge is * that the top can be occulted by the top edge of the disk * and by star 1. * **********************/ if( strcmp( mode.disk, "ON") == 0 ) { /************************************* * * First do the edge of the disk. The exchange and then * re-exchange of elements is done to avoid writing yet * another special function to handle the edge. * **************************************/ if( strcmp( mode.star2, "ON") == 0 ) { if( fabs( orbphase ) <= 0.25 ) { max = 3; for( ir = 0; ir < max; ir++) { for( iphi = 0; iphi <= maxphiindexdisk; iphi++) { z = xdisk[ir][iphi]; xdisk[ir][iphi] = xdisk[maxrindex - ir][iphi]; xdisk[maxrindex - ir][iphi] = z; z = ydisk[ir][iphi]; ydisk[ir][iphi] = ydisk[maxrindex - ir][iphi]; ydisk[maxrindex - ir][iphi] = z; z = visibledisk[ir][iphi]; visibledisk[ir][iphi] = visibledisk[maxrindex - ir][iphi]; visibledisk[maxrindex - ir][iphi] = z; } } keepside = 1; transmission = 0.0; occult( keepside, transmission, max, maxphiindexdisk, xdisk, ydisk, visibledisk, nedge2, edge2x, edge2y); for( ir = 0; ir < max; ir++) { for( iphi = 0; iphi <= maxphiindexdisk; iphi++) { z = xdisk[ir][iphi]; xdisk[ir][iphi] = xdisk[maxrindex - ir][iphi]; xdisk[maxrindex - ir][iphi] = z; z = ydisk[ir][iphi]; ydisk[ir][iphi] = ydisk[maxrindex - ir][iphi]; ydisk[maxrindex - ir][iphi] = z; z = visibledisk[ir][iphi]; visibledisk[ir][iphi] = visibledisk[maxrindex - ir][iphi]; visibledisk[maxrindex - ir][iphi] = z; } } } } /******************************** * * Now do the top of the disk. * ****************************/ max = maxrindex - 3; if( strcmp( mode.star1, "ON") == 0 ) { keepside = 1; transmission = 0.0; occult( keepside, transmission, max, maxphiindexdisk, xdisk, ydisk, visibledisk, nedge1, edge1x, edge1y); } keepside = 0; transmission = disktop.transmission; occult( keepside, transmission, max, maxphiindexdisk, xdisk, ydisk, visibledisk, ntopedgedisk, topedgediskx, topedgedisky); if( strcmp( mode.star2, "ON") == 0 ) { if( fabs( orbphase ) <= 0.25 ) { keepside = 1; transmission = 0.0; occult( keepside, transmission, max, maxphiindexdisk, xdisk, ydisk, visibledisk, nedge2, edge2x, edge2y); } } } /******************************** * * Finally, do star 1. * *********************************/ if( strcmp( mode.star1, "ON" ) == 0 ) { keepside = 0; transmission = 0.0; occult( keepside, transmission, maxthetaindex1, maxphiindex1, x1, y1, visible1, nedge1, edge1x, edge1y); if( strcmp( mode.star2, "ON") == 0 ) { if( fabs( orbphase ) <= 0.25 ) { keepside = 1; transmission = 0.0; occult( keepside, transmission, maxthetaindex1, maxphiindex1, x1, y1, visible1, nedge2, edge2x, edge2y); } } } /******************************* * * Add all the fluxes from the visible grid points. * ********************************/ flux = 0.0; if( strcmp( mode.star1, "ON") == 0 ) { for( itheta = 0; itheta <= maxthetaindex1; itheta++) { for( iphi = 0; iphi <= maxphiindex1; iphi++) { if( visible1[itheta][iphi] > 0.0 ) { flux = flux + visible1[itheta][iphi] * SpecInt( "star1", star1.fluxsource, localT1[itheta][iphi], mu1[itheta][iphi] ) * mu1[itheta][iphi] * dS1[itheta][iphi]; } } } } if( strcmp( mode.star2, "ON") == 0 ) { for( itheta = 0; itheta <= maxthetaindex2; itheta++) { for( iphi = 0; iphi <= maxphiindex2; iphi++) { if( visible2[itheta][iphi] > 0.0 ) { flux = flux + visible2[itheta][iphi] * SpecInt( "star2", star2.fluxsource, localT2[itheta][iphi], mu2[itheta][iphi] ) * mu2[itheta][iphi] * dS2[itheta][iphi]; } } } } if( strcmp( mode.disk, "ON") == 0 ) { for( iphi = 0; iphi <= maxphiindexdisk; iphi++) { for( ir = 0; ir <= maxrindex; ir++) { if( visibledisk[ir][iphi] > 0.0 ) { if( ir <= (maxrindex-3) ) { flux = flux + visibledisk[ir][iphi] * SpecInt( "disktop", disktop.fluxsource, localTdisk[ir][iphi], mudisk[ir][iphi] ) * mudisk[ir][iphi] * dSdisk[ir][iphi]; } else { flux = flux + visibledisk[ir][iphi] * SpecInt( "diskrim", diskrim.fluxsource, localTdisk[ir][iphi], mudisk[ir][iphi] ) * mudisk[ir][iphi] * dSdisk[ir][iphi]; } } } } } return( flux ); } void coordsys( double *lsky, double *msky, double *nsky, double inclination, double orbphase) /***************** * * This function calculates the unit vectors along the axes * of the sky Cartesian coordinate system in terms of the star * coordinate system. * Input data: * inclination = orbital inclination in degrees. * orbphase = orbital phase. Phase 0.0 is defined to be * inferior conjunction of star 2. * must lie between -0.5 and 0.5. * Output data: * lsky[], * msky[], * nsky[] = unit vectors of the sky coordinate system, where * * lsky[] is in both the plane of the sky and the plane of the orbit * msky[] makes a right-handed, orthonormal coordinate system and * is always positive in the positive z direction. * nsky[] points to the sun * The system rotates counter-clockwise seen from above. * * Note the special handling of orbphase = +/-0.25. * Note also that the orientation at inclination = 0 joins smoothly * onto the orientation at inclination > 0 because the dummy * variable "dummy" is set equal to sin()/cos() instead of ysun/xsun. * *****************/ { double pi, radincl, radphase, xsun, ysun, zsun, dummy, alpha, beta; pi = 3.14159265358979323; radincl = 2.0 * pi * ( inclination / 360.0); radphase = 2.0 * pi * orbphase; xsun = sin( radincl ) * cos( pi - radphase ); ysun = sin( radincl ) * sin( pi - radphase ); zsun = cos( radincl ); nsky[1] = xsun; nsky[2] = ysun; nsky[3] = zsun; if( orbphase == -0.25 ) { beta = 0.0; alpha = 1.0; } else if( orbphase == 0.25) { beta = 0.0; alpha= -1.0; } else { dummy = sin( pi - radphase ) / cos( pi - radphase ); beta = sqrt( 1.0 / (1.0 + dummy*dummy) ); if( (orbphase > -0.25) && (orbphase < 0.25) ) beta = -beta; alpha = -beta * dummy; } lsky[1] = alpha; lsky[2] = beta; lsky[3] = 0.0; msky[1] = -beta * zsun; msky[2] = alpha * zsun; msky[3] = beta*xsun - alpha*ysun; return; } void project( double *pxsky, double *pysky, double xstar, double ystar, double zstar, double *lsky, double *msky) /************************** * * This function projects a position vector in the star coordinate * system (xstar, ystar, zstar) to a position vector in the sky * coordinate system (xsky, ysky, zsky). * ***************************/ { *pxsky = xstar * lsky[1] + ystar * lsky[2] + zstar * lsky[3]; *pysky = xstar * msky[1] + ystar * msky[2] + zstar * msky[3]; return; } void occult( long keepside, double transmission, long maxfirstindex, long maxsecondindex, double x[SIZEDIMEN1][SIZEDIMEN2], double y[SIZEDIMEN1][SIZEDIMEN2], double visible[SIZEDIMEN1][SIZEDIMEN2], long nedge, double *edgex, double *edgey) /**************************** * * This function finds which points in an object are visible * and which have been occulted by an object with an * outline described by the points edgex, edgey. * Input data: * keepside = 0, keep the points on the inside of the curve. * 1, keep the points on the outside of the curve. * x[][], y[][] = the points on the surface of the star projected * to the plane of the sky. * visible[][] = array holding the fractional visibility of each * points. * NOTE: this function will only decrease * visible[][] from 1.0 towards 0.0, never increase * it. Thus, the effect of this function is to * make points cummulatively less visible. * maxfirstindex, * maxsecondindex = the maximum values of the indices of x[][], y[][], * and visible[][]. * edgex, edgey = the points describing the outline of * the occulting object. * nedge = the number of points in the edge. * transmission = the amount of light that is transmitted by the * occulting object. * 0.0 = no light transmitted * 1.0 = all light transmitted * * Output data: * visible[][] = 1.0, if the point is fully visible, * 0.0, if the point is fully invisible, * and between 1.0 and 0.0 for special cases. * ****************************/ { double xlow, xhigh, ylow, yhigh, centerx, centery, xa, ya, rsq, rsqmin, rsqmax, rmin, epsilonr; double twopi, deltar, edger[LINEPOINTS], edgetheta[LINEPOINTS]; long int i, j, maxindex; twopi= 2.0 * 3.14159265358979323846; if( nedge == 0 ) return; /******************************* * * First put the edge into a polar coordinate system with * its center at the center of the occulting object and then * resort for safety. * **************************************/ xlow = edgex[0]; xhigh = edgex[0]; ylow = edgey[0]; yhigh = edgey[0]; for( i = 0; i < nedge; i++) { if( edgex[i] < xlow ) xlow = edgex[i]; if( edgex[i] > xhigh ) xhigh = edgex[i]; if( edgey[i] < ylow ) ylow = edgey[i]; if( edgey[i] > yhigh ) yhigh = edgey[i]; } centerx = (xlow + xhigh) / 2.0; centery = (ylow + yhigh) / 2.0; for( i = 0; i < nedge; i++) { xa = edgex[i] - centerx; ya = edgey[i] - centery; edger[i] = sqrt( xa*xa + ya*ya ); edgetheta[i] = 0.0; if( edger[i] != 0.0 ) edgetheta[i] = acos( xa /edger[i] ); if( ya < 0.0 ) edgetheta[i] = twopi - edgetheta[i]; } sort( nedge, edgetheta, edger ); rsqmin = edger[0]; rsqmax = edger[0]; for( i = 1; i < nedge; i++) { if( edger[i] < rsqmin ) rsqmin = edger[i]; if( edger[i] > rsqmax ) rsqmax = edger[i]; } rsqmin = 0.9 * rsqmin * rsqmin; rsqmax = 1.00001 * rsqmax * rsqmax; /************************** * The factor of 0.04 is based on a few numerical experiments. * Any value between about 0.02 and 0.1 seems to give excellent * results. **************************/ rmin = star1.frontradius; if( star2.frontradius < rmin ) rmin = star2.frontradius; maxindex = maxfirstindex; if( maxsecondindex > maxindex ) maxindex = maxsecondindex; epsilonr = 0.04 * ( 2.0 / maxindex ) * rmin; /******************************* * * Calculate whether each point is on the inside or the outside * of the occulting curve, and multiply visible[][] by "transmission" * if it is on the occulted side of the curve. * If the point is closer to the line than epsilonr, * change visible[][] to 0.5*(1 + transmission) of its previous value * to mimic a half-hidden point. * ***********************************/ for( i = 0; i <= maxfirstindex; i++) { for( j = 0; j <= maxsecondindex; j++) { if( visible[i][j] > 0.0 ) { xa = x[i][j] - centerx; ya = y[i][j] - centery; rsq = xa*xa + ya*ya; if( rsq > rsqmax ) deltar = 1.0; else if( rsq < rsqmin ) deltar = -1.0; else deltar = edgedistance( nedge, edger, edgetheta, xa, ya); if( (keepside == 1) && (deltar < -epsilonr ) ) { visible[i][j] = transmission * visible[i][j]; } if( (keepside == 0) && (deltar > epsilonr ) ) { visible[i][j] = transmission * visible[i][j]; } if( fabs( deltar ) <= epsilonr ) { visible[i][j] = 0.5 * (1.0 + transmission) * visible[i][j]; } } } } return; } double planck( double temperature, double wavelength, char *mode ) /******************** * * Returns the flux emitted by the surface of a black body. * Input data: * temperature in degrees Kelvin * wavelength in centimeters * mode "NU" returns Fnu in units of erg/sec/cm^2/Hz * "LAMBDA" returns Flambda in units of erg/sec/cm^2/A * * The Planck function is normalized such that * * (integral over frequency) = sigma * T**4 * * Thus, it is the monochromatic total flux per unit area * *********************/ { double nu, x, y, Fnu, Flambda; /******************** * * twopi = 2.0 * 3.14159265358979; * h = 6.626076e-27 erg s * c= 2.9979e+10 cm/s * k = 1.38066e-16 erg/K * 2 * pi * h / c^2 = 4.6323525e-47 * h / k = 4.79920911e-11 * *********************/ if( temperature <= 3.0 ) { Fnu = 0.0; } else if( wavelength <= 1.0e-8 ) { Fnu = 0.0; } else { nu = 2.9979e+10 / wavelength; x = 4.6323525e-47 * nu * nu * nu; y = (4.79920911e-11 * nu) / temperature ; if( y > 100.0 ) { Fnu = 0.0; } else { Fnu = x / ( exp( y ) - 1.0 ); } } if( strcmp( mode, "NU" ) == 0 ) { return( Fnu ); } else if( strcmp( mode, "LAMBDA") == 0 ) { /***************** * * Calculate Flambda from Fnu and then * convert from erg/cm^2/sec/cm to erg/cm^2/sec/A * *****************/ Flambda = (Fnu * 2.9979e+10) / ( wavelength * wavelength ); Flambda = Flambda * 1.0e-8; return( Flambda ); } else quit("Unrecognized mode in function planck()."); return( -1.0 ); } double SpecInt( char *object, char *fluxsource, double temperature, double mu) /************************* * * This function returns the specific intensity in the direction * of the observer. It has two parts * 1) The temperature dependent part, which can come either from * adopting a monochromatic black body flux or from tables. * 2) The effect of limb darkening. * ****************************/ { double T[500], Izero[500], a[500], b[500]; double u, norm, inten, wavelength, delta, Io, alpha, beta, x; long int i, ilow, n; /***************************** * * If the emitted radiation is assumed to be Black Body radiation, * the specific intensity is given by * * inten = Planck(T,lambda,mode) * limbdarkening * * where the function Planck() returns the total monochromatic FLUX, * not the monochromatic SPECIFIC INTENSITY. * * The limb darkening is assumed to be linear. * * ld = norm * ( 1 - u + u*mu) * * and the normalization factor is appropriate for conversion of * flux to intensity; that is, * * integral( ld * mu * d(omega)) = 1 * * So, norm = (1 / pi) / (3 / (3 - u)) * ******************************/ if( strcmp( fluxsource, "BB") == 0 ) { wavelength = starsys.wavelength; if( strcmp( object, "star1") == 0 ) u = star1.bblimbdark; else if( strcmp( object, "star2") == 0 ) u = star2.bblimbdark; else if( strcmp( object, "disktop") == 0 ) u = disktop.bblimbdark; else if( strcmp( object, "diskrim") == 0 ) u = diskrim.bblimbdark; else { quit("Unrecognized object in SpecInt."); } norm = 0.954929658 / (3.0 - u); inten = planck( temperature, wavelength, "LAMBDA" ) * norm * (1.0 - u + u*mu); return( inten ); } /***************************** * * If the specific intensity comes from a table look up, * the specific intensity is given by * * inten = Io * ld * * where Io is the specific intensity at mu = 1, and * the functional form of the limb darkening is * * ld = ( 1 - a*(1-mu) - b*(1-mu)^2 ) * * and mu, a, and b have their usual meanings. Note that * the normalization factor on the limb darkening is 1.0 * because one is assumed to be given Io = I(mu = 1). * *****************************/ if( strcmp( object, "star1") == 0 ) { for( i = 0; i < Itabn1; i++ ) { T[i] = ItabT1[i]; Izero[i] = ItabIo1[i]; a[i] = ItabA1[i]; b[i] = ItabB1[i]; } n = Itabn1; } else if( strcmp( object, "star2") == 0 ) { for( i = 0; i < Itabn2; i++ ) { T[i] = ItabT2[i]; Izero[i] = ItabIo2[i]; a[i] = ItabA2[i]; b[i] = ItabB2[i]; } n = Itabn2; } else if( strcmp( object, "disktop") == 0 ) { for( i = 0; i < Itabntop; i++ ) { T[i] = ItabTtop[i]; Izero[i] = ItabIotop[i]; a[i] = ItabAtop[i]; b[i] = ItabBtop[i]; } n = Itabntop; } else if( strcmp( object, "diskrim") == 0 ) { for( i = 0; i < Itabnrim; i++ ) { T[i] = ItabTrim[i]; Izero[i] = ItabIorim[i]; a[i] = ItabArim[i]; b[i] = ItabBrim[i]; } n = Itabnrim; } else quit("Unrecognized object in SpecInt."); /************************************ * * Now interpolate/extrapolate in the tables to get the * correct intensity and limb darkening coefficients. * **************************************/ if( temperature <= T[0] ) { Io = Izero[0]; alpha = a[0]; beta = b[0]; } else if( temperature >= T[n-1] ) { Io = Izero[n-1]; alpha = a[n-1]; beta = b[n-1]; } else { for( i = 0; i < (n-1); i++ ) { if( (temperature >= T[i]) && (temperature < T[i+1]) ) { ilow = i; break; } } delta = ( temperature - T[ilow] ) / ( T[ilow+1] - T[ilow] ); Io = delta * ( Izero[ilow+1] - Izero[ilow] ) + Izero[ilow]; alpha = delta * ( a[ilow+1] - a[ilow] ) + a[ilow]; beta = delta * ( b[ilow+1] - b[ilow] ) + b[ilow]; } x = 1 - mu; inten = Io * ( 1.0 - alpha*x - beta*x*x ); return( inten ); }