In the last post I wrote about my attempt to estimate crater radius and center location using Sartopo data. Lets review.
I use this hypothetical circle to idealize what what a crater would look like. In reality, it isn’t a perfect circle. We can see in topography data that the rims are not equal distances from the center. The point closest to the center is dc, which separates the left and right into two right triangles. Each side can be thought of as some fraction x or y or the total profile length Dp. We can define Dp as follows.
Script based minimizing
I approached this thinking, how can I minimize these two equations? I thought, x, r_left, r_right, and d_c are all solvable if we assume we know where the center is. So I created a script to solve for Dp at each possible point of the center. It seemed promising, where each equation seemed to produce a circle of points that approached the desired Dp/ Subtracting the two from one another, which should produce a 0, would hone in on exactly where they meet, focusing on two points, one near the center found using the radar image and one far outside of view. Unfortunately, it was too far off. I thought it was simply a mistake somewhere in the script, so I proceeded to create a write up starting by reworking through the math I used.
We already have the equation for Dp above, but we also can find the other variables,
such that Dp=Dp. The calculations will be slightly rounded, but the net result approaches Dp because I use circular logic. This was not what I was getting. I was perplexed, and after staring for a while, Kevin lent me a hand. He pointed out that I was multiplying by x instead of dividing. That concluded this failed attempt.
Minimizing with MATLAB function
Beginning again with Dp,
This still leaves us with an function that doesn’t use Dp. I made a mistake in my logic, which I only just realized now, even after writing it up.
I moved x and y to the other side of the equation, but also kept it on the left side. This lead to incorrect and inconsistent results.
Luckily, the steps I used, of which I have more to say (later in this post) are already in place, so I can easily adjust the function and run it again. Where the new function is found to be as follows.
Keep in mind, x can be solved as a function of the other 3 variables are a matter of geometry. The other three numbers have a more direct relation to the center as a function of Dp. To my surprise, the initial results are more promising.
Literature put the center at the red dot in the center, and this estimate puts it at the yellow dot with a red rim. There is a bit more I’d like to do such as play with the constraints and plot a circle to see how well it fits the rim and radar image.
Visualizing and translating the results
This is a necessary thing to discuss because it took several days to take these variables and output values for the center lat and long. I ended up working this out multiple times from different points of view to find the center.
I thought the logic was sound, but I kept getting results that did not make sense. I have x and y values for each point on the profile, but when I found the distance between points using the x and y values, they did not translate into lat and long conversions appropriately. The biggest problem is that the y points don’t change as much as the longitude. X appears to change consistently with longitude, but I found it still to be inaccurate.
I compared changes in longitude with changes in x. I converted each of these to the other and found the two did not translate. The obvious question to ask is how I am converting between the two. Distances in kilometers were converted to lat and longs with,
rLon = dx/ (R *pi()/180 * cos(deg2rad(latitude)));
%where the inputed latitude for the center is used
rLat = dy / (pi()/180 * R);
I wondered if there was something about this equation that wasn’t consistent or that didn’t apply for Titan, so I tried to compare it with another method that goes from changes in lat/lon to changes in kilometers. The method comes from one Andrew Hedges. I chose this because its the recommended method from JPL.
dlon = lon2 – lon1
dlat = lat2 – lat1
a = (sin(dlat/2))^2 + cos(lat1) * cos(lat2) * (sin(dlon/2))^2
c = 2 * atan2( sqrt(a), sqrt(1-a) )
d = R * c
I wondered if I had inputed this wrong, so I changed R to that of earth and tested it with the calculator at the site I link to above and it matched. I used this to calculate the profile length and it was off ~11km from the dbidr. Using just along x, I got a number ~5 times higher than the longitude and nearly as much higher than the dbidr value, meaning its an order of magnitude larger. dbidr ~=77km, lat/long~=66km, dx ~=300km. Therefore, lat long more reliable than the x and y values.
I don’t know if I am wrong to use these lat and lon conversions, but the x and y numbers do not work for solving in the x and y domain. So I used the initial lat and long values to to define the left rim, the dc intersection, and the center in x and y. Changes in lat and long between the left rim and dc, plus the results of the minimization give new x and y for the center. I can then find the change in x and y between dc and the center, or from lr and the center. I did both and got the same resulting lat and lon.
I feel like there must be something about the nature of the sartopo x and y data that I do not understand because the alternative is that it is wrong. Obviously the error is probably with me, but I have gone through this so many times in some many ways to try and hone in on the problem and its lead me here.
Its also worth noting, I wrote the script to use the dx/dy, dlat/dlon depending on the slope of the profile and its general location relative to the center because that will influence whether we are moving to higher or lower lat/lons. Meaning, it isn’t about how I am using the changes in x/y.