I recently wanted to know how Maxima (a computer algebra system) implemented something. I searched and searched and couldn't solve my problem. Specifically, I wanted to access a Maxima structure from within Lisp code. However, I didn't know how the structure was implemented in Lisp and therefore didn't know how to access it.
All you need is a single line of inline Lisp:
structures is a global variable which stores the structures that have been defined in the session by using defstruct(). I wasted a few hours looking for this information to come up with 25 characters that would answer my question. Applying the same general idea to an instance of PanelStruct() tells me how it is implemented. By the way, the answer to my questions looks like this:
The near repeat is caused by the print function returning what it has printed and Maxima outputting it. The $ indicates a references to a Maxima variable and the | | symbols indicate a case sensitive reference. Oddly, the case is reversed for RDV (which is documented) but not for PanelStruct which surprises me. (These structures are not built-in but were user-defined in my Maxima session.)
I will still need to figure out how I'm going to use the information, but I have a working concept.
Saturday, July 6, 2013
Wednesday, February 6, 2013
Lighting Distribution
I've written a simple Maxima program which uses linear interpolation and lamp luminous intensity data to determine the illuminance on a flat surface (read floor or work surface) caused by several lamps in different locations. It displays a color contour do show you how the light is distributed. I could probably tweak some things to get a nicer selection of colors for the color contour but as far as a result which is easy to interpret, this works fine. Behold:
Never mind the clutter in the bottom left hand corner. This is just the z-axis numbers written on top of one another because we are looking straight down on the graph. Small price to pay for the cool result and you can edit the picture afterward if you need something for presentation purposes. An additional problem I observe is the aspect ratio. It should look square because my chosen room is square (8' × 8'). Dauntless, I press on.
To implement this program with your own example you need:
Light Selection
You need to select a light and find the manufacturer data. I am using the Cooper RPN3MR-E3MRC. This is a combination of the light fixture and the lamp in the fixture. The data you need will indicate how the candle power changes with vertical angle. The assumed situation is that the light is pointing down. The vertical angle is measured between the vertical line passing through the center of the light fixture and the line drawn through the center of the fixture to a location on the floor (or incident surface). This vertical angle affects the "amount of light" (measured in candelas) that is going in that direction. To input the data for your luminaire-lamp combination you enter it on the line that says "lampData:". The data is in the form of a list of lists where the nested lists (two items each) represent points. The pattern for each point in the master list is [angle in degrees, candle power in candelas]. If your data doesn't go all the way to 90 degrees, then put a final point [90, 0] to indicate there will be no light going in that direction.
Considerations in the Calculation
The vertical angle not only affects the "amount of light" in that direction but the angle of incidence with the floor. We are assuming the light is "facing" the floor. (The floor is horizontal. The light is shining straight down onto it, although, of course it "spreads out" as it goes down.) Remembering that the "angle of incidence" is measured relative to the normal line (perpendicular to the floor), you can see that the angle of incidence is an alternating angle with the vertical angle and therefore equal to it. An angle of incidence affects how the light "spreads out" across a surface because the same amount of light gets spread out more if the surface is at an angle (not "facing head on") to the beam of light. Also, not every point of the floor is the same distance away from each light. As light travels it spreads out in two directions. If you spread the same amount of anything (sound, light, peanut butter) over a larger area, there isn't as much of it in a given area. The basic formula is:
\[fc = \frac{cp}{d^2}\cos{\theta}\],
where
fc = foot-candles (a measure of illuminance: cd/ft2)
cp = candle power (candelas: cd, same unit in metric and imperial, relates to the Watt)
d = distance (in feet)
θ = angle of incidence (measured from the normal).
The formula is the same if you use metric, but the units change. Instead of entering feet and getting foot-candles, you enter meters and get lux. You use the same lamp data which is in candelas. Candelas are the same in metric and imperial. To put it very simply:
lux = cd/m2
fc = cd/ft2
Using the Program
If you want metric, enter all distances in m (you will get lux). If you want imperial, enter all distances in ft (you will get foot-candles).
Note: P2P_Illum() means Point-to-Point Illuminance. Illum() is the sum of the illuminance contributions of each luminaire to a given point.
Possible Improvements to the Program
Never mind the clutter in the bottom left hand corner. This is just the z-axis numbers written on top of one another because we are looking straight down on the graph. Small price to pay for the cool result and you can edit the picture afterward if you need something for presentation purposes. An additional problem I observe is the aspect ratio. It should look square because my chosen room is square (8' × 8'). Dauntless, I press on.
To implement this program with your own example you need:
- Maxima: a free computer algebra system.
- Light data
- A room and light configuration.
Light Selection
You need to select a light and find the manufacturer data. I am using the Cooper RPN3MR-E3MRC. This is a combination of the light fixture and the lamp in the fixture. The data you need will indicate how the candle power changes with vertical angle. The assumed situation is that the light is pointing down. The vertical angle is measured between the vertical line passing through the center of the light fixture and the line drawn through the center of the fixture to a location on the floor (or incident surface). This vertical angle affects the "amount of light" (measured in candelas) that is going in that direction. To input the data for your luminaire-lamp combination you enter it on the line that says "lampData:". The data is in the form of a list of lists where the nested lists (two items each) represent points. The pattern for each point in the master list is [angle in degrees, candle power in candelas]. If your data doesn't go all the way to 90 degrees, then put a final point [90, 0] to indicate there will be no light going in that direction.
Considerations in the Calculation
The vertical angle not only affects the "amount of light" in that direction but the angle of incidence with the floor. We are assuming the light is "facing" the floor. (The floor is horizontal. The light is shining straight down onto it, although, of course it "spreads out" as it goes down.) Remembering that the "angle of incidence" is measured relative to the normal line (perpendicular to the floor), you can see that the angle of incidence is an alternating angle with the vertical angle and therefore equal to it. An angle of incidence affects how the light "spreads out" across a surface because the same amount of light gets spread out more if the surface is at an angle (not "facing head on") to the beam of light. Also, not every point of the floor is the same distance away from each light. As light travels it spreads out in two directions. If you spread the same amount of anything (sound, light, peanut butter) over a larger area, there isn't as much of it in a given area. The basic formula is:
\[fc = \frac{cp}{d^2}\cos{\theta}\],
where
fc = foot-candles (a measure of illuminance: cd/ft2)
cp = candle power (candelas: cd, same unit in metric and imperial, relates to the Watt)
d = distance (in feet)
θ = angle of incidence (measured from the normal).
The formula is the same if you use metric, but the units change. Instead of entering feet and getting foot-candles, you enter meters and get lux. You use the same lamp data which is in candelas. Candelas are the same in metric and imperial. To put it very simply:
lux = cd/m2
fc = cd/ft2
Using the Program
If you want metric, enter all distances in m (you will get lux). If you want imperial, enter all distances in ft (you will get foot-candles).
- As always, choose a coordinate system and stick to it. Pick a corner of your room to be (0,0) and one of the sides to be the positive x-axis.
- Enter the length and width of your room. I have called these x_length and y_length to avoid misunderstandings based on the usual understanding of width < length. Just pick a dimension to be the x direction length and one for the y direction length. Doesn't matter which is set to which. Just be consistent with step 1. Sorry, only rectangles. You're welcome to program something more interesting if you get bored some day. Or, maybe your room can pretend to be a rectangle?
- Enter the locations (positions) of your luminaires. These are in terms of [x, y, z]. I have entered the height of the luminaires for the the z value (12.5 feet). You could instead enter their height above a work surface (the results would only be correct for the work surface). Alternatively, you can change the Illum(x,y) function definition to use a z-value of the height of the work surface instead of 0 (same caveat as previous suggestion).
Note: P2P_Illum() means Point-to-Point Illuminance. Illum() is the sum of the illuminance contributions of each luminaire to a given point.
Possible Improvements to the Program
- Change the color selections for the plot. Perhaps yellow for higher values. Might be a fiddly thing to get right.
- Include a more explicit boundary (wall) definition mechanism. A list of vertices would satisfy the wall definition needs. But then some programming would need to be done to determine which lights have a line of sight to a given x,y point. There would also be a need to either crop the end result or indicate some neutral value for places in the plot area but not within the boundaries of the room.
- Allow the incident surface shape to be specified.
- This might be as simple as allowing inputs of areas at different heights. This would be the simplest and most widely useful approach. It would be helpful to indicate the locations of the different surfaces on the plot as well. It might be necessary to directly use the draw package to do this. [First, load(draw). Then learn how to use draw. I haven't looked much into it yet.]
- More complex surface entry would require more complex angle of incidence calculations. The angle of incidence would not necessarily be equal to the vertical angle. If the incident surface (or part thereof) was represented as an analytic function of two variables, you could use the partial derivatives to determine the tangent plane and thus the direction of the normal. Given the direction of the normal and the direction toward a given luminaire, you take the dot product of these (normalized/unit) vectors to get the cosine of the angle of incidence (which is what you need; see formula above).
Thursday, December 6, 2012
Degree of Saturation versus Relative Humidity
These two quantities are similar and as I worked through ASHRAE Handbook Fundamentals (2009) IP, chapter 1, I was initially baffled as to what the difference was. First, the definitions based on ASHRAE Fundamentals:
Relative humidity: The mole fraction of water vapour in a sample divided by the mole fraction of water vapour in saturated air at the same temperature and pressure.
Degree of saturation: The humidity ratio of water vapour in a sample divided by the humidity ratio of water vapour in saturated air at the same temperature and pressure.
Humidity ratio: The humidity ratio is the mass of water vapour divided by the mass of dry air in a sample. This value can also be expressed in terms of molar fractions.
Discussion of Differences
In relative humidity, imagine taking two samples, one of the air which you want to know the relative humidity of and one of saturated air at the same temperature and pressure. For example, we might take samples of 1 mole each. In the non-saturated air, we will have less water vapour than in the saturated air. On the other hand, we will have more dry air in the non-saturated sample than in the saturated sample. The same general relationship will be true of mass, but the masses and the mole fractions will relate differently due to differing molecular masses. Here's the bottom line: we compare the samples based on moles (which amounts to the number of molecules) not on a given mass or volume. When you're dealing with gases, moles is the way to go. In other words, our non-saturated and saturated samples relate by the equation:
where the items on the left hand side are the mole fractions of water vapour and dry air in the non-saturated sample and the items on the right hand side are for saturated air.
Understanding that we need to compare these quantities with a fixed number of total moles as in the previous equation is what makes ASHRAE's equation (14) work out. From the last given equation we have
Therefore,
which verifies equation (14) as given in ASHRAE Fundamentals (2009).
It is worth noting that we only had to worry about this molar equality in developing this equation. When it comes to computing the degree of saturation from the relative humidity and vice versa, we don't have to think about that. The formula will do that for us. Also, we can still compute the degree of saturation directly from values off of a psychrometric chart (for example) by reading across to the right to get \(W\)and finding \(W_S\) by finding where the dry bulb temperature intersects the saturation curve and reading across to the right. The definition of degree of saturation doesn't require us to account for the number of moles, only the relative humidity (and yes, if you're using charts you can read it off there as well). If you are familiar with reading psychrometric charts, you will know that reading W values is generally much easier than reading RH values.
So if you're getting different values for these quantities in your work, it's not (necessarily) a mistake or a misreading of the chart—they really are different.
Relative humidity: The mole fraction of water vapour in a sample divided by the mole fraction of water vapour in saturated air at the same temperature and pressure.
Discussion of Differences
In relative humidity, imagine taking two samples, one of the air which you want to know the relative humidity of and one of saturated air at the same temperature and pressure. For example, we might take samples of 1 mole each. In the non-saturated air, we will have less water vapour than in the saturated air. On the other hand, we will have more dry air in the non-saturated sample than in the saturated sample. The same general relationship will be true of mass, but the masses and the mole fractions will relate differently due to differing molecular masses. Here's the bottom line: we compare the samples based on moles (which amounts to the number of molecules) not on a given mass or volume. When you're dealing with gases, moles is the way to go. In other words, our non-saturated and saturated samples relate by the equation:
Understanding that we need to compare these quantities with a fixed number of total moles as in the previous equation is what makes ASHRAE's equation (14) work out. From the last given equation we have
It is worth noting that we only had to worry about this molar equality in developing this equation. When it comes to computing the degree of saturation from the relative humidity and vice versa, we don't have to think about that. The formula will do that for us. Also, we can still compute the degree of saturation directly from values off of a psychrometric chart (for example) by reading across to the right to get \(W\)and finding \(W_S\) by finding where the dry bulb temperature intersects the saturation curve and reading across to the right. The definition of degree of saturation doesn't require us to account for the number of moles, only the relative humidity (and yes, if you're using charts you can read it off there as well). If you are familiar with reading psychrometric charts, you will know that reading W values is generally much easier than reading RH values.
So if you're getting different values for these quantities in your work, it's not (necessarily) a mistake or a misreading of the chart—they really are different.
Convective Heat Transfer on a Building Envelope (Wind Chill?)
I'd like to make a brief mathematical investigation into the concept of a wind chill factor and how one particular scheme common in popular culture (and even affecting some in industry) does not relate well to modelling of a building envelope. Wind chill is the temperature which a human being perceives accounting for the effects of wind and temperature together. The preceding sentence is awkward because it is trying to say too much at once, so let me try again:
hNW AS (TS – TWC) = q = hactual AS (TS – Tactual),
which simplifies to
TWC = TS – (hactual / hNW)(TS – Tactual).
A Little Problem
This leaves me with a problem. I want to have some type of reasonable value for the convective heat transfer coefficient ratio (hactual / hNW). How do I get that? I will engage in a little pragmatism and pretend that the wind chill equation will still give me a way to get an (hactual / hNW) ratio. This is perhaps a dubious step, but may suffice for illustrative purposes. From the wind chill equation we get a 15 mph wind at -31°F being equivalent to no wind at -59°F. Here's what we get for the ratio for human skin (which I have to use since that's what the wind chill equation is based on):
(hactual / hNW) = (TS – TWC) / (TS – Tactual)
= (88 – (-59)) / (88 – (-31))
= 1.235
(I've used 88°F for human skin, but that temperature will be different depending on the ambient temperature. Oh, well, what do you do when you live in a shoe?)
Theoretically, the convective heat transfer coefficient is largely independent of temperature and area, but it will change based on the interaction of the fluid and the surface which will include things like geometry, orientation, surface roughness, moisture, etc. Unfortunately, temperature differences may affect the actual coefficient. They will affect the result if you consider radiation losses as well, since radiation transport is proportional to absolute temperature to the fourth power.
Wind Chill Effect on Surfaces at Different Temperatures Than Human Skin?
So, let's consider a surface which is maintained at a temperature of 10°F and is surrounded by an ambient temperature of -31°F and the wind speed is 15 mph. (This is not the temperature in the boundary layer. The point of the convective heat transfer equation is to deal with this phenomenon in a simplified way—essentially to circumvent it.) So, what's my "equivalent temperature in the absence of wind"?
TWC = TS – (hactual / hNW)(TS – Tactual)
= 10 – 1.235 (10 – (-31))
= -40.6 °F
So, the rate of heat loss from Surface X under the actual conditions is the same as the rate of heat loss for Surface X with no wind and a temperature of -40.6°F.
We obtained a different value than we did for the wind chill index applicable to human skin. The reason is straightforward: the rate of heat loss is proportional to the difference between the ambient air temperature and the surface temperature. Now, if it is -31°F outside, what temperature is the outside surface of a building? Perhaps -21°F? This makes quite a difference:
TWC = TS – (hactual / hNW)(TS – Tactual)
= -21 – 1.235 (-21 – (-31))
= -33.4 °F
So the effect is not zero, but it isn't much to talk about. If my h ratio is actually much higher, we will get a greater difference.
Considerations on Building Outside Surface Temperatures
The thickness and properties of the building envelope materials may change the relationship between ambient air temperature and surface temperature. Such changes will change the rate of heat loss in proportion to the temperature difference. Is there a higher difference at glazing surfaces than concrete surfaces?
There's also a sticky point regarding solar heat gains. Such gains on the surface will increase the rate of heat loss at the surface. But watch out! This is mainly behaviour occurring at the surface. We're dealing with a cold climate situation. Radiant heat from the sun heats up the surface only to be removed by increased rates of convection and radiation loss. This does not help us analyze the heat losses from the interior of the building! Plain "horse sense" tells me that adding thermal energy to the exterior does not cause a net increase to the rate of heat loss from the building interior. At worst, the increased rate of convective (and radiation) loss will remove all of the radiant energy added. To deal thoroughly with this issue, we would either need to model both together (allowing them to interact numerically with each other—for transient analysis) or use a "net direct solar gain" approach that deals with the interaction of convection and radiation (gains and losses) at the surface in an integrated way. The obvious outcome is that the rate of heat loss from the interior is reduced, but determining how much could be fun times.
Note that in warm climate conditions (cooling conditions), heat on the outside surface migrates toward the interior and we must therefore consider it. It is no longer starting to go in only to double back* (as it were), but actually conducting through the envelope and so directly affecting interior conditions. (*Analogical language enjoys talking about physics.)
CAUTION: This is NOT a Final Answer
In closing, this article is not trying to present the final answer to this problem. Rather I am attempting to expose the fallacy that the wind chill index (and corresponding formula) can be naively applied to energy modelling problems in building science. Furthermore, my criticism is not directed against nor does it address considerations of the effects of wind on air leakage analysis.
Note that ASHRAE Fundamentals chapter 26 (2009) gives effective R-values which can be used. The point of these values is to bypass all of the stuff I'm writing about in this article and account for convection and radiation losses at the surface under a 15 mph wind design condition. You don't actually need all that wind chill stuff—chapter 26 R-values already give a reasonable estimate of the wind effects, i.e., R-0.17 at 15 mph on the exterior instead of R-0.68 on the interior for still air at a wall.
A wind chill index (expressed in units equivalent(ish) to temperature measurements) is answering a question: What temperature TWC (with wind speed = 0) will be perceived by a human being as equivalent to temperature Tactual (with wind speed = Vactual)? (We are talking here more particularly about the exposed skin of a human being. Also, it isn't purely a matter of perception.)
Human beings, not buildings. Consider the following differences:
Human beings, not buildings. Consider the following differences:
- buildings (in cold climates) have a much lower surface temperature than humans
- the amount (and temperature) of moisture on their surfaces is generally different and has a different internal transport mechanism (brick has basically no moisture and that moisture is already cold, human skin normally has moisture—perspiration—which evaporates more rapidly in windy conditions than calm conditions causing heat loss, because evaporation requires energy)
- unlike human beings, buildings do not have psychological comfort issues affecting their perception of temperature (and, admittedly, cannot reasonably be considered to perceive anything)
The equation used today for wind chill index has some theoretical basis but in the end is empirically derived. That means they did experiments and found something that "works". It is a very pragmatic way to deal with complex issues, especially ones that include psychological/cognitive components. (I'm not "disrespecting" empirical equations, they're useful.)
So, where do I start? I would hard pressed with my current knowledge/understanding to deal with the evaporative cooling effects on human skin as compared with buildings, though such investigation might be meritorious. I will not attempt to compare the psychology of buildings with that of humans—I hope you understand. I will however seek to show how the concept of an "equivalent temperature in the absence of wind" yields different results depending on the surface temperature of the material in question. As such, although we will not compare human skin with face brick, we will compare warm Surface X with cool Surface X.
Convective Heat Transfer
The convective heat transfer equation is simple enough:
q = hc AS (TS – Tair)
where q is the heat energy (per unit time), hc is the heat transfer coefficient, AS is the surface area, TS is the temperature of the surface and Tair is the ambient air temperature.
The heat transfer coefficient varies with wind speed. Suppose we want to find an air temperature TWC (no wind, hNW) which causes the same rate of heat loss as Tactual (wind speed = Wactual, hactual). In that case, we need to relate:
Convective Heat Transfer
The convective heat transfer equation is simple enough:
q = hc AS (TS – Tair)
where q is the heat energy (per unit time), hc is the heat transfer coefficient, AS is the surface area, TS is the temperature of the surface and Tair is the ambient air temperature.
The heat transfer coefficient varies with wind speed. Suppose we want to find an air temperature TWC (no wind, hNW) which causes the same rate of heat loss as Tactual (wind speed = Wactual, hactual). In that case, we need to relate:
hNW AS (TS – TWC) = q = hactual AS (TS – Tactual),
which simplifies to
TWC = TS – (hactual / hNW)(TS – Tactual).
A Little Problem
This leaves me with a problem. I want to have some type of reasonable value for the convective heat transfer coefficient ratio (hactual / hNW). How do I get that? I will engage in a little pragmatism and pretend that the wind chill equation will still give me a way to get an (hactual / hNW) ratio. This is perhaps a dubious step, but may suffice for illustrative purposes. From the wind chill equation we get a 15 mph wind at -31°F being equivalent to no wind at -59°F. Here's what we get for the ratio for human skin (which I have to use since that's what the wind chill equation is based on):
(hactual / hNW) = (TS – TWC) / (TS – Tactual)
= (88 – (-59)) / (88 – (-31))
= 1.235
(I've used 88°F for human skin, but that temperature will be different depending on the ambient temperature. Oh, well, what do you do when you live in a shoe?)
Theoretically, the convective heat transfer coefficient is largely independent of temperature and area, but it will change based on the interaction of the fluid and the surface which will include things like geometry, orientation, surface roughness, moisture, etc. Unfortunately, temperature differences may affect the actual coefficient. They will affect the result if you consider radiation losses as well, since radiation transport is proportional to absolute temperature to the fourth power.
Wind Chill Effect on Surfaces at Different Temperatures Than Human Skin?
So, let's consider a surface which is maintained at a temperature of 10°F and is surrounded by an ambient temperature of -31°F and the wind speed is 15 mph. (This is not the temperature in the boundary layer. The point of the convective heat transfer equation is to deal with this phenomenon in a simplified way—essentially to circumvent it.) So, what's my "equivalent temperature in the absence of wind"?
TWC = TS – (hactual / hNW)(TS – Tactual)
= 10 – 1.235 (10 – (-31))
= -40.6 °F
So, the rate of heat loss from Surface X under the actual conditions is the same as the rate of heat loss for Surface X with no wind and a temperature of -40.6°F.
We obtained a different value than we did for the wind chill index applicable to human skin. The reason is straightforward: the rate of heat loss is proportional to the difference between the ambient air temperature and the surface temperature. Now, if it is -31°F outside, what temperature is the outside surface of a building? Perhaps -21°F? This makes quite a difference:
TWC = TS – (hactual / hNW)(TS – Tactual)
= -21 – 1.235 (-21 – (-31))
= -33.4 °F
So the effect is not zero, but it isn't much to talk about. If my h ratio is actually much higher, we will get a greater difference.
Considerations on Building Outside Surface Temperatures
The thickness and properties of the building envelope materials may change the relationship between ambient air temperature and surface temperature. Such changes will change the rate of heat loss in proportion to the temperature difference. Is there a higher difference at glazing surfaces than concrete surfaces?
There's also a sticky point regarding solar heat gains. Such gains on the surface will increase the rate of heat loss at the surface. But watch out! This is mainly behaviour occurring at the surface. We're dealing with a cold climate situation. Radiant heat from the sun heats up the surface only to be removed by increased rates of convection and radiation loss. This does not help us analyze the heat losses from the interior of the building! Plain "horse sense" tells me that adding thermal energy to the exterior does not cause a net increase to the rate of heat loss from the building interior. At worst, the increased rate of convective (and radiation) loss will remove all of the radiant energy added. To deal thoroughly with this issue, we would either need to model both together (allowing them to interact numerically with each other—for transient analysis) or use a "net direct solar gain" approach that deals with the interaction of convection and radiation (gains and losses) at the surface in an integrated way. The obvious outcome is that the rate of heat loss from the interior is reduced, but determining how much could be fun times.
Note that in warm climate conditions (cooling conditions), heat on the outside surface migrates toward the interior and we must therefore consider it. It is no longer starting to go in only to double back* (as it were), but actually conducting through the envelope and so directly affecting interior conditions. (*Analogical language enjoys talking about physics.)
CAUTION: This is NOT a Final Answer
In closing, this article is not trying to present the final answer to this problem. Rather I am attempting to expose the fallacy that the wind chill index (and corresponding formula) can be naively applied to energy modelling problems in building science. Furthermore, my criticism is not directed against nor does it address considerations of the effects of wind on air leakage analysis.
Note that ASHRAE Fundamentals chapter 26 (2009) gives effective R-values which can be used. The point of these values is to bypass all of the stuff I'm writing about in this article and account for convection and radiation losses at the surface under a 15 mph wind design condition. You don't actually need all that wind chill stuff—chapter 26 R-values already give a reasonable estimate of the wind effects, i.e., R-0.17 at 15 mph on the exterior instead of R-0.68 on the interior for still air at a wall.
Tuesday, November 6, 2012
Linear Interpolation Made Convenient in Maxima
(I usually put the load statement in a separate cell in wxMaxima so that I am not re-executing the load every time I correct my main code and rerun the cell. Minor point.)
Maxima can turn a set of points into a function for linear interpolation. For that matter, it can make a cubic spline out of the same points. I'm more interested in the linear interpolated function today—the other is for your reference. There are three main points of interest in the above code:
- linearinterpol() takes points, sorts them by the first ordinate, and produces an expression.
- linearinterpol() produces expressions which use the characteristic function, rather than taking cases to make it transparent to other Maxima routines.
- The use of the double single quote operator ('') and the single quote operator (') helps me to define the function f(x) such that x is a formal variable and not embedded as an expression.
- linearinterpol() produces an expression which extrapolates beyond the specified points (from negative infinity to positive infinity) and it is left to the user's discretion to determine the limits to which this extrapolation is reasonable or applicable.
In a Nutshell
All you need is a list of 2D points. There is more than one way to specify these, but I find it easiest to use a list of lists. The main list elements are all lists of 2 representing x and y ordinates. The function sorts these points according to the x-coordinate before proceeding. You can specify the variable which will be used in the expression. Provided that you precede that variable with a single quote operator, you can use a variable that is elsewhere assigned a value. More on that below.
Transparent Output
All you need is a list of 2D points. There is more than one way to specify these, but I find it easiest to use a list of lists. The main list elements are all lists of 2 representing x and y ordinates. The function sorts these points according to the x-coordinate before proceeding. You can specify the variable which will be used in the expression. Provided that you precede that variable with a single quote operator, you can use a variable that is elsewhere assigned a value. More on that below.
Transparent Output
The characteristic function, which is often assigned the Greek letter chi ("χ"), returns either 0 or 1. It is a function of two variables, a set and a value. If the value is an element of the set, it returns 1, otherwise it returns 0. In the context of our linear interpolation expression, we have a linear expression for each interval in our set of points and each of these is multiplied by a characteristic function with a different set. Since the sets used for the characteristic functions form a partition of the real numbers, we won't have more than one of the terms of our overall expression from linearinterpol() evaluate to a non-zero value at a time. This expression can be integrated, differentiated, or combined with other expressions without having to take special cases. Take the monkey off your back, let Maxima do it.
Cool, Useful, but Initially Abstruse Operators
The single quote operator (beside the key) is a way to tell Maxima, "Don't evaluate this thing." If I write x and x is assigned value, then Maxima will do what is normally very helpful to me, namely, make the substitution for me. If I don't want Maxima to do that, because I actually want the expression x, not the value that x stands for (because it has been assigned a value), I need to write 'x instead. In contrast, the double single quote operator ('') tells Maxima to replace the given expression with the value of the expression. This way, when I use the delayed assignment operator (:=) it does not delay the execution of the linearinterpol() function, only the evaluation of the expression as a function of x is delayed. You may recall from my previous post that parameters in function declarations are formal. Even if x is elsewhere assigned a value, f(x) := some_expression_including_x + x, will still work intuitively.
Here is some sample input followed by the output. This is the best way to clarify in your mind the use of these operators:
x:2;
y:2*x;
z: 2*'x;
d(x) := z;
d(2);
d(3);
h(x) := ''(z);
h(2);
h(3);
(%o46) 2
(%o47) 4
(%o48) 2*x
(%o49) d(x):=z
(%o50) 2*x
(%o51) 2*x
(%o52) h(x):=2*x
(%o53) 4
(%o54) 6
Extrapolation
The expression that linearinterpol() produces does not protect you from yourself. If will produce a result from negative to positive infinity without telling you that you have gone beyond the boundaries of your initial data set. Actually, this is as it should be, in my opinion. The program gives you the freedom to ask all the what-ifs your little heart desires, but it is up to you to discern what is really true to life and applicable to the particular problem you are trying to solve—whether the extrapolation is truly justified by the data is rightly left as a matter for user judgement.
Cool, Useful, but Initially Abstruse Operators
The single quote operator (beside the
Here is some sample input followed by the output. This is the best way to clarify in your mind the use of these operators:
x:2;
y:2*x;
z: 2*'x;
d(x) := z;
d(2);
d(3);
h(x) := ''(z);
h(2);
h(3);
(%o46) 2
(%o47) 4
(%o48) 2*x
(%o49) d(x):=z
(%o50) 2*x
(%o51) 2*x
(%o52) h(x):=2*x
(%o53) 4
(%o54) 6
Extrapolation
The expression that linearinterpol() produces does not protect you from yourself. If will produce a result from negative to positive infinity without telling you that you have gone beyond the boundaries of your initial data set. Actually, this is as it should be, in my opinion. The program gives you the freedom to ask all the what-ifs your little heart desires, but it is up to you to discern what is really true to life and applicable to the particular problem you are trying to solve—whether the extrapolation is truly justified by the data is rightly left as a matter for user judgement.
Wednesday, October 17, 2012
Making a Double Out of Parts in VB.NET
Recently I picked up on a piece of code I had started several months ago and didn't have time to finish. I was a bit stumped actually, though I've got it pretty well beat into shape now (maybe more on that later). I was making a VB.NET type Fraction using BigIntegers as the base type. As such, the class would theoretically be able to represent any rational number. But what if I wanted to convert that result to a double? So, I started into the code necessary to do the job. I began with the easy part, the cases where the numerator and denominator can both be converted to doubles. In that case, .NET already knows how to do the needed conversions. But if the numerator or denominator is too large to be converted to a double, the result of division (mathematically) may still be within the range of a double. For example,
4.445 × 10400 / 4.222 × 10399 = 10.528185693983894
An important finishing piece to solving this programming problem is the following function:
First, for a good summary of the double format, see Double-precision floating-point format (Wikipedia).
Code Walkthrough
The double type has 64 bits as does the Int64 type. We begin with initializing our result variable to 0. Depending on whether the number is indicated as positive or not, we set the sign bit. The sign bit (bit 63, the most significant bit) needs to be set (1) to indicate a negative number. For simplicity of use, I wanted to let the caller work in terms of the exponent (base 2) rather than with a biased result. The 1023 "excess" or bias is one of those values you don't want to have to remember. This function obviates the need to be overly conscious of that value.
After making our biased exponent (biasedExp), we check that it is in range. We have exactly 11 bits to work with. If our result is greater than 11 bits (0x7FF in hex, 11111111111 in binary), then we complain to the caller with an OverflowException. If the result is negative, then it still does not fit into the 11 bits and we complain with an OverflowException in this case as well. The mantissa is required to be no more than 52 bits. That's 13 "nibbles" (4 bits, a single hex digit) or 6 1/2 bytes. If the mantissa has any bits set in the upper 12 bits, we complain. By calling the variable mantissa I have already declared my draconian intent with this function. I will not coddle the caller by shifting things that are in the wrong place to the right place. If the caller sends a mantissa which is too large, they may have done the shifting wrong or they may have incorrectly included the implicit 1 at the beginning of the "binary scientific notation." Either way, the caller has made an error and we cannot know which one it is. If we could know which error they made, it might make sense to correct it silently and move on. However, a wrong assumption here could be somebody's program's undoing. Therefore, we just complain and let the caller figure out what to do about it.
The biased exponent is now shifted into place—this, like the 1023 excess, falls into the category of "painful detail we want to encapsulate so we can forget about it." (This is not ambiguous like the too-large-mantissa.) In our final steps we "Or" the biased exponent and mantissa into place and use a call to BitConverter.Int64BitsToDouble() to ensure we have the proper format.
Room For Improvement
To make this function more robust, an exception class could be defined, either several different ones (one for each error) or a class with properties and enumerations to allow the caller to efficiently test for the particular condition that caused the "overflow". I note also that I have used OverflowException where an underflow has occurred. In program debugging, an accurate error message ("over" versus "under") can provide a clue to the person debugging the program.
4.445 × 10400 / 4.222 × 10399 = 10.528185693983894
An important finishing piece to solving this programming problem is the following function:
First, for a good summary of the double format, see Double-precision floating-point format (Wikipedia).
Code Walkthrough
The double type has 64 bits as does the Int64 type. We begin with initializing our result variable to 0. Depending on whether the number is indicated as positive or not, we set the sign bit. The sign bit (bit 63, the most significant bit) needs to be set (1) to indicate a negative number. For simplicity of use, I wanted to let the caller work in terms of the exponent (base 2) rather than with a biased result. The 1023 "excess" or bias is one of those values you don't want to have to remember. This function obviates the need to be overly conscious of that value.
After making our biased exponent (biasedExp), we check that it is in range. We have exactly 11 bits to work with. If our result is greater than 11 bits (0x7FF in hex, 11111111111 in binary), then we complain to the caller with an OverflowException. If the result is negative, then it still does not fit into the 11 bits and we complain with an OverflowException in this case as well. The mantissa is required to be no more than 52 bits. That's 13 "nibbles" (4 bits, a single hex digit) or 6 1/2 bytes. If the mantissa has any bits set in the upper 12 bits, we complain. By calling the variable mantissa I have already declared my draconian intent with this function. I will not coddle the caller by shifting things that are in the wrong place to the right place. If the caller sends a mantissa which is too large, they may have done the shifting wrong or they may have incorrectly included the implicit 1 at the beginning of the "binary scientific notation." Either way, the caller has made an error and we cannot know which one it is. If we could know which error they made, it might make sense to correct it silently and move on. However, a wrong assumption here could be somebody's program's undoing. Therefore, we just complain and let the caller figure out what to do about it.
The biased exponent is now shifted into place—this, like the 1023 excess, falls into the category of "painful detail we want to encapsulate so we can forget about it." (This is not ambiguous like the too-large-mantissa.) In our final steps we "Or" the biased exponent and mantissa into place and use a call to BitConverter.Int64BitsToDouble() to ensure we have the proper format.
Room For Improvement
To make this function more robust, an exception class could be defined, either several different ones (one for each error) or a class with properties and enumerations to allow the caller to efficiently test for the particular condition that caused the "overflow". I note also that I have used OverflowException where an underflow has occurred. In program debugging, an accurate error message ("over" versus "under") can provide a clue to the person debugging the program.
Tuesday, October 9, 2012
How to make a .NET (Excel) Add-In Using Visual Basic Express 2010
Automation add-ins are easier than what follows and are the way to go if you just want a few custom functions. Actually, for efficiency's sake, you might want to just use VBA functions added to Module1 (must be Module1—no renaming to a more suitable name) and then save the workbook as an Excel Add-In (*.xlam). But, if it is user interaction or in some way making Excel "do stuff" (rather than just return a single result or an array of results), an add-in which implements the
IDTExtensibility2 and IRibbonExtensibility interfaces may be a good strategy. (The main alternative is to use VBA.)
The functions RegisterFunction and UnregisterFunction are invoked by utilities like regasm.exe which comes with the .NET Framework. I have created a batch file which I run in order to register the assembly. You will need to modify it to suit your assembly.
Note that once the Add-in is registered and added, you will not need to reregister it after every Build (while developing).
Making an Excel Add-in in VB.NET involves the following (the ordering of the steps is not absolute):
- Ensure you have Office installed on the development machine and have the Primary Interop Assemblies (PIAs) installed in the Global Assembly Cache (GAC). You can download them from Microsoft. (The Microsoft installer will put the PIAs in the GAC.)
- Start a Class Library project in VB.NET.
- It's a good idea to generate your own strong name using the sn.exe tool (part of the .NET framework). You only need one for multiple projects. In a nutshell, you use the sn.exe tool to create a *.snk file which you will tell VB.NET about.
- let VB.NET know where your *.snk file is by going to Project/Properties..., under Signing check "Sign the assembly" and choose a strong name key file. You will be able to browse to where your *.snk file is
- if your add-in is intended to be widely distributed, and especially if you want to sell your add-in, you would be well-advised to learn more about this topic to ensure your end-users don't get warning messages that cause them to be (perhaps unduly) alarmed.
- Apparently, step 4 is unnecessary, but I did step 4, so here it is: In the project properties, go to Application and click the button Assembly Information.... In this screen, check the Make assembly COM-Visible box.
- Create a GUID. There are several ways but a program called guidgen.exe may already be on your computer which will do the job. Just search your computer for it. It is located under C:\Program Files\Microsoft SDKs\Windows\... on my computer.
- Copy and Paste the GUID into the Assembly Information... dialog box referred to in step 4.
- Encapsulate your project in a namespace either by typing it explicitly in the files or by specifying a root namespace in the project properties.
- Add references to needed PIAs, and extensibility (.NET tab)
- It is recommended that you add these by referencing the Interop Assembly on the .NET tab, not the COM tab. You may need to browse for this: Look under C:\Windows\assembly\gac_msil\Microsoft.Office.Interop.Excel\...
- Create a class named Connect (seems that the name has to be Connect) which implements the IDTExtensibility2 and IRibbonExtensibility interfaces.
- Implements Extensibility.IDTExtensibility2, Microsoft.Office.Core.IRibbonExtensibility
- Implement the methods of the interfaces.
- Decorate the class with COM visibility specifications and include functions for creating and deleting custom registry settings.
- Register your assembly with regasm.exe.
Steps 8, 9, and 10 are implemented in the following code snippet, which has been parameterized somewhat.
The functions RegisterFunction and UnregisterFunction are invoked by utilities like regasm.exe which comes with the .NET Framework. I have created a batch file which I run in order to register the assembly. You will need to modify it to suit your assembly.
Note that once the Add-in is registered and added, you will not need to reregister it after every Build (while developing).
Subscribe to:
Posts (Atom)