With our polar coordinates (r, ν) in hand, there are a few more transformations.
The argument of latitude u is obtained by adding the argument of perigee ω to ν. Argument of perigee (or more generally, periapsis) is taken from the GPS almanac.
u = ν + ω
Next we convert to cartesian coordinates to get position in the orbital plane:
x' = r * cos(u)
y' = r * sin(u)
Now we need to account for rotation of the earth. The right ascension of the ascending node Ω0 is given in the almanac, and it is valid at the beginning of the referenced GPS week, not at the time of applicability. So we add the rotation speed of the earth times the number of seconds since the start of the week.
Ω = Ω0 - Ωe(dot) * t
This will rotate our coordinates so we can get a longitude of the satellite relative to earth longitude. The quantity Ωe(dot) is the rotation rate of the earth, 7.292115 x 10-5 radians per second.
With all these new quantities in hand, we'll use a rotation matrix to move the orbit into earth coordinates: latitude and longitude.
Monday, February 6, 2012
Thursday, February 2, 2012
More on Orbits, Part 3
So far we haven't accounted for the eccentricity of the orbit. We can get this number from the appropriate line in the almanac file. A value of zero means the orbit is circular. As the value moves towards 1, the orbit is more and more enlongated. GPS satellites have very circular orbits, with values less than 0.02. Comets typically have values much closer to 1.
From mean anomaly, we find the eccentric anomaly E, from mean anomaly M that we found in the previous post, and eccentricity ε.
M = E - ε sin E
This equation needs to be solved iteratively. Newton's method can be used to get a very accurate value after a few loops, by finding the root of the equation:
f(E) = E - ε sin E - M = 0
The derivative with respect to E is f '(E) , or: 1 - ε cos E
We can start with a first guess E0 = M. Take E1 = E0 - f(E0)/f '(E0)
So E1 = E0 - (E0 - ε sin E0 - M) / (1 - ε cos E0)
Then find E2 with E1 in the above function, and so on. You can check En in the original equation and keep going until En causes M to converge to the correct value.
E gives us the angle from the center of the ellipse to the satellite. Note that the center is the point between the two foci of the ellipse. The earth will be at one focus, not at the center.
At this point we can get the radius r which represents the distance of the satellite from the center of the earth.
r = a * (1 - ε * cos E)
Finally, to get the angle of the satellite position relative to earth, we convert to true anomaly ν.
cos ν = (cos E - ε) / (1 - ε * cos E)
Now we have the position in polar coordinates (r, ν) of the satellite in the orbital plane. There is still a little more work to get a current latitude and longitude position.
From mean anomaly, we find the eccentric anomaly E, from mean anomaly M that we found in the previous post, and eccentricity ε.
M = E - ε sin E
This equation needs to be solved iteratively. Newton's method can be used to get a very accurate value after a few loops, by finding the root of the equation:
f(E) = E - ε sin E - M = 0
The derivative with respect to E is f '(E) , or: 1 - ε cos E
We can start with a first guess E0 = M. Take E1 = E0 - f(E0)/f '(E0)
So E1 = E0 - (E0 - ε sin E0 - M) / (1 - ε cos E0)
Then find E2 with E1 in the above function, and so on. You can check En in the original equation and keep going until En causes M to converge to the correct value.
E gives us the angle from the center of the ellipse to the satellite. Note that the center is the point between the two foci of the ellipse. The earth will be at one focus, not at the center.
At this point we can get the radius r which represents the distance of the satellite from the center of the earth.
r = a * (1 - ε * cos E)
Finally, to get the angle of the satellite position relative to earth, we convert to true anomaly ν.
cos ν = (cos E - ε) / (1 - ε * cos E)
Now we have the position in polar coordinates (r, ν) of the satellite in the orbital plane. There is still a little more work to get a current latitude and longitude position.
Wednesday, February 1, 2012
More on Orbits, Part 2
I should add a note here about how to identify the satellites. All the satellites transmit on the same frequency, and each satellite has a unique pseudo-random noise (PRN) code that spreads its message data over a wide bandwidth. This is a direct sequence spread spectrum method known as CDMA. There are 32 unique codes so there will be up to 32 active satellites, and they are identified in the almanac by their PRN codes.
The GPS almanac references its time relative to the start of a GPS week. This week starts at midnight UTC time on the first day of the calendar week (Saturday night turning into Sunday morning.) The orbital parameters given in the almanac are valid at the reference time, or epoch, which is given as "Time of Applicability". This time is the number of seconds after the start of the referenced GPS week.
The satellite position is given as mean anomaly at the epoch, on the line in the almanac file that starts with "Mean Anom". Mean anomaly is an angle measurement that increases linearly over time. We get the current mean anomaly using the difference between our time of interest t and the epoch toe:
tk = t - toe
Then the current mean anomaly is
M = M0 + n * tk
There is our mean motion n from the previous post. Now we have an angle measurement that will help us find the position of the satellite in its orbit at time t.
The GPS almanac references its time relative to the start of a GPS week. This week starts at midnight UTC time on the first day of the calendar week (Saturday night turning into Sunday morning.) The orbital parameters given in the almanac are valid at the reference time, or epoch, which is given as "Time of Applicability". This time is the number of seconds after the start of the referenced GPS week.
The satellite position is given as mean anomaly at the epoch, on the line in the almanac file that starts with "Mean Anom". Mean anomaly is an angle measurement that increases linearly over time. We get the current mean anomaly using the difference between our time of interest t and the epoch toe:
tk = t - toe
Then the current mean anomaly is
M = M0 + n * tk
There is our mean motion n from the previous post. Now we have an angle measurement that will help us find the position of the satellite in its orbit at time t.
More on Orbits
Awhile back, I wrote about my orbit simulator. This has been on my mind again, as I'm dealing with some issues at work that have to do with GPS receivers. GPS is useful to synchronize widely dispersed systems such as cellular base stations. It not only provides location, but also very accurate timing. In fact, it is the timing that allows you to triangulate your exact position.
My concern is in predicting which satellites are in view at any given time. There are some websites that track satellite orbits and predict passes over your location, but I'm looking for something more specific to what I want. How to go about this? It turns out that a lot of the same equations I used in my orbit simulator come into play for predicting locations of real world satellites.
The first good reference is the ICD-GPS-200C document which is pretty easy to find on the web. The equations starting in Table 20-IV are what we need.
Since the satellite orbits are perturbed by various influences, the GPS almanac is updated daily. The orbits are generally stable, but will have some small variance from day to day. Almanacs are available here.
The first equation I'll deal with is the mean motion, n.
n = sqrt(GMe / a3)
GMe is 3.986 x 1014 and relates to the gravitational constant and the mass of the earth. The semi-major axis a has to do with the size of the elliptical orbit. The almanac gives the square root of a in the line that starts with SQRT(A).
The mean motion relates to the inverse of the orbital period and gives the number of radians per second that a given satellite is moving around its orbit. This will be important for finding the position at a given time, because the almanac just gives a snapshot of where it was, or will be, at some reference time.
That is all for this post. The next step is finding the position in orbit at a particular time.
My concern is in predicting which satellites are in view at any given time. There are some websites that track satellite orbits and predict passes over your location, but I'm looking for something more specific to what I want. How to go about this? It turns out that a lot of the same equations I used in my orbit simulator come into play for predicting locations of real world satellites.
The first good reference is the ICD-GPS-200C document which is pretty easy to find on the web. The equations starting in Table 20-IV are what we need.
Since the satellite orbits are perturbed by various influences, the GPS almanac is updated daily. The orbits are generally stable, but will have some small variance from day to day. Almanacs are available here.
The first equation I'll deal with is the mean motion, n.
n = sqrt(GMe / a3)
GMe is 3.986 x 1014 and relates to the gravitational constant and the mass of the earth. The semi-major axis a has to do with the size of the elliptical orbit. The almanac gives the square root of a in the line that starts with SQRT(A).
The mean motion relates to the inverse of the orbital period and gives the number of radians per second that a given satellite is moving around its orbit. This will be important for finding the position at a given time, because the almanac just gives a snapshot of where it was, or will be, at some reference time.
That is all for this post. The next step is finding the position in orbit at a particular time.
Wednesday, March 16, 2011
SystemVerilog struct assignment
You may have occasion to initialize a structure in SystemVerilog:
This might seem fine, but the above assignment is actually a concatenation. The simulator will take the set of values and pack them into one big bit vector.
In my experience, Modelsim was okay with positive values, or a mix of positive and negative values. In compilation, it would note that it was promoting an assignment to a concatenation. But if all the values were negative, it would give an error that a packed value was being assigned to an unpacked value.
A quick Google search turned up a book preview that showed an example of how to initialize a structure. Very simple, really. Just add a tick (') before the assignment:
This does the trick, and the notes and errors go away. SystemVerilog uses the '{} construct to differentiate a list of values from a concatenation, {}.
EDIT: Another example, this time with a dynamic array or queue:
typedef struct {
int f1;
int f2;
int f3;
} set_of_values_T;
set_of_values_T set_of_values = {1, 2, -3};
This might seem fine, but the above assignment is actually a concatenation. The simulator will take the set of values and pack them into one big bit vector.
In my experience, Modelsim was okay with positive values, or a mix of positive and negative values. In compilation, it would note that it was promoting an assignment to a concatenation. But if all the values were negative, it would give an error that a packed value was being assigned to an unpacked value.
A quick Google search turned up a book preview that showed an example of how to initialize a structure. Very simple, really. Just add a tick (') before the assignment:
set_of_values_T set_of_values = '{1, 2, 3};
This does the trick, and the notes and errors go away. SystemVerilog uses the '{} construct to differentiate a list of values from a concatenation, {}.
EDIT: Another example, this time with a dynamic array or queue:
bit [0:2] values[$]= '{3,5,6};
Tuesday, February 22, 2011
Engineering Economics
I'm going to deviate from technical discussion and get into money a little. Engineers are often under pressure to reduce costs. Someone asked for advice this week on removing a couple of parts that are not needed for a particular version of a product that we are shipping. Removing those parts has some implications, which include testing to verify that there are no "gotchas", possibly releasing a new program image to account for the parts not being there, further testing of that new image, and so on.
It's good for engineers to have at least some grasp of the tradeoffs involved in their decisions. In this case, we are talking about 14 cents of savings per board. At 5,000 units per year, this comes out to $700. (Network infrastructure does not have the volumes of consumer products.) At 10,000 we are at $1400. Is it worth the time?
How much is engineering time worth? You could try to figure the hourly cost. Divide the salary by 2,000 hours per year. Don't forget to add the cost to the employer of benefits, employer's share of taxes, worker's comp, office space, supplies, all that expensive depreciating equipment in the lab. $100/hour is probably low. What about $500? Maybe.
What else could an R&D engineer be spending his time on? What about new product development?
How much does an engineer bring to the bottom line in new product development? How much of the profits should be allocated to engineering? We grudgingly have to give some credit to marketing, legal, human resources, and other groups that keep the company going. Say engineers make up a third of the overhead. Then maybe it's fair to take credit for a third of the net profit.
How much profit for a new product? We can look at net margins for the business and figure what percentage of overall business profits is due to the new product.
The higher the volumes, the more it makes sense to pay someone to work on the problem, whether it's taking 14 cents of parts off a board, or developing a new product. We engineers are often isolated from these kinds of numbers, especially in a large company. We rely on upper management to set our priorities. But we can at least throw around some rough numbers and figure out if we are pulling our weight.
It's good for engineers to have at least some grasp of the tradeoffs involved in their decisions. In this case, we are talking about 14 cents of savings per board. At 5,000 units per year, this comes out to $700. (Network infrastructure does not have the volumes of consumer products.) At 10,000 we are at $1400. Is it worth the time?
How much is engineering time worth? You could try to figure the hourly cost. Divide the salary by 2,000 hours per year. Don't forget to add the cost to the employer of benefits, employer's share of taxes, worker's comp, office space, supplies, all that expensive depreciating equipment in the lab. $100/hour is probably low. What about $500? Maybe.
What else could an R&D engineer be spending his time on? What about new product development?
How much does an engineer bring to the bottom line in new product development? How much of the profits should be allocated to engineering? We grudgingly have to give some credit to marketing, legal, human resources, and other groups that keep the company going. Say engineers make up a third of the overhead. Then maybe it's fair to take credit for a third of the net profit.
How much profit for a new product? We can look at net margins for the business and figure what percentage of overall business profits is due to the new product.
The higher the volumes, the more it makes sense to pay someone to work on the problem, whether it's taking 14 cents of parts off a board, or developing a new product. We engineers are often isolated from these kinds of numbers, especially in a large company. We rely on upper management to set our priorities. But we can at least throw around some rough numbers and figure out if we are pulling our weight.
Friday, February 11, 2011
Arbitrary voltage source in LTSpice
Here is a post I found when searching for a way to use some data as an input to LTSpice to test a filter.
http://ltspicelabs.blogspot.com/2006/10/using-wav-files-for-io-and-transient.html
This talks about taking a .wav file as input. Which is fine if you are doing audio. But what if you want to generate a .wav file based on, say, a mathematical model?
Python has a nice wave module. So you could write a short script like this:
This will write out your data to a .wav file that can then be used in LTSpice. The example here has a sample width of 2, which would be two bytes, or signed 16 bits. I was modeling a signed 8-bit DAC output, so I scaled it up and offset it to get what I wanted. Also I put in an oversampling factor to get a more step-wise response instead of a linear ramp between sample points.
Don't forget "import wave" in the Python script!
http://ltspicelabs.blogspot.com/2006/10/using-wav-files-for-io-and-transient.html
This talks about taking a .wav file as input. Which is fine if you are doing audio. But what if you want to generate a .wav file based on, say, a mathematical model?
Python has a nice wave module. So you could write a short script like this:
# Send data to output .wave file
# filename = string
# fsys = integer system clock rate
# data = list of integers
def writewave(filename, fsys, data):
oversample_factor = 8
offset_factor = 16384
scale_factor = 60 # 8 bit in to ~15 bit out
w = wave.open(filename, 'wb')
w.setnchannels(1)
w.setsampwidth(2)
w.setframerate(fsys * oversample_factor) # need to oversample to generate the steppy DAC output
w.setnframes(len(data))
d = '' # string to hold byte values in little-endian order
for i in data:
ii = scale_factor * i + offset_factor
rl = int(ii) & 255
rh = (int(ii) >> 8) & 255
for j in xrange(oversample_factor):
d = d + chr(rl) + chr(rh) # little-endian
w.writeframes(d)
w.close()
This will write out your data to a .wav file that can then be used in LTSpice. The example here has a sample width of 2, which would be two bytes, or signed 16 bits. I was modeling a signed 8-bit DAC output, so I scaled it up and offset it to get what I wanted. Also I put in an oversampling factor to get a more step-wise response instead of a linear ramp between sample points.
Don't forget "import wave" in the Python script!
Subscribe to:
Posts (Atom)