Development and Evaluation of New Methods for Estimating Albedo-Areas for Three-Axis Stabilized Geosynchronous Satellites

Although direct measurements of the projected areas of various Geosynchronous Earth Orbit (GEO) satellite facets are impossible without high-resolution imaging, estimates of the albedo-Area (aA) product lead to the possibility of inferring the geometric size of an RSO based on the illuminated area. Such size estimates are an integral part of its identity because satellite mission capabilities can be inferred. We are engaged in parallel development of two methods for calculating aA for the body/communication antennas structures and one method for the solar panels. We have previously reported on the Two Facet Model (2FM) method for body aA, and here we discuss a method based on differences between new observations and a baseline catalog that has been constructed from the GEO Observations with Longitudinal Diversity Simultaneously (GOLDS) data. We report on evaluations of the 2FM and Differential Method (DM) algorithm results. We also discuss a new method of estimating solar panel aA by fitting new data that include specular glints. All of these measurement methods are compared to models and simulations that serve as a proxy for ground truth. Because of the partially directional nature of the composite Bi-directional Reflectivity Distribution Function (BRDF) of all bus-mounted appendages, variance of body aA results is expected to be significant. Short-term and long-term variance of the derived aAs will also be discussed.


Introduction
We have been conducting both observational and analysis programs for geostationary satellites. The primary databases that we have created and used for algorithmic development are the Geo Color Photometry Catalog (GCPC, Gregory and Payne) and the Geo Observations with Longitudinal Diversity Simultaneously (GOLDS, Gregory, et al.). The observations were carried out with a variety of commercial, offthe-shelf (COTS), Raven-like sensors that are primarily located in continental US (CONUS) and Maui. The photometric observations were obtained through the Johnson (or its derivatives) filter system, with most of the data being in the B and R filters, thus providing both brightness and color information. Several analyses of the data have been conducted and show, among other important findings, that there are important diurnal, seasonal, and geographic dependences of both brightness and color. The research presented here is aimed at deriving effective estimates of the sizes of these objects by calculating the albedo-Area (aA) products of the bus/payload/com antennas and solar panels based on "new data," using the GCPC and GOLDS catalogs to develop our algorithms. In general, these new methods can be used with any new ground-or spacebased sensors and can be generalized to different filter/detector combinations. Sizes are an important intrinsic characteristic of a satellite because inferrences to satellite power or payloads can be made.

Satellite-Based Coordinates
One of our goals has been to remove as much of the sensor location and seasonal variations as possible by carrying out our calculations in a coordinate system that is based on the Resident Space Object (RSO) itself. The manner in which these groundbased geographic and seasonal effects contribute to observed brightness is determined from the calculation of several body frame angles.
In defining the RSO reference frame, the Geo object is in a nadir-aligned stabilization (i.e., the normal of the nadir-facing facet is anti-parallel to the vector that points from Earth center to the satellite). Figure 1 shows the geometry of the RSO coordinate frame. The surface represents the nadir-pointing facet of the object, and we identify the surface normal n with the body axis z B . For a satellite in orbit, the direction of its velocity vector, |v|, lies in the orbital plane, as does n, and so the cross product n x |v| is normal to n and points to orbital South. We identify this vector as y B . If the object is in a perfectly circular orbit (e ≡ 0), the velocity vector v would complete the triad, and the body frame would be complete; however, for e ≠ 0, the third reference frame axis, x B , is defined by y B x z B ≡ (n x |v|) x n. Figure 2 shows the four body angles that will be used in the work described herein; these are the Sun zenith and azimuthal angles (θ s , Φ s , respectively) and observer zenith and azimuthal angles (θ o , Φ o , respectively). Zenith angles are the separations between the surface normal and the vectors that point to the Sun or observer. Azimuthal angles are measured between one of the axes (x B , in the figure) and the projection of the Sun, or observer, vector onto the facet surface. The sign of these angles (i.e., the direction of rotation from the x B axis) must be taken into account. For our purposes, one only needs to calculate |ΔΦ o | = |Φ s -Φ o |, the absolute angular difference between the projections of the Sun and observer vectors. The vectors shown in Figs 1 and 2 are not fixed in inertial space but instead evolve over time.

The Two Facet Model (2FM) Method for Body aA
The measured brightness of a three-axis stabilized geosynchronous satellite can be decomposed into the contribution due to the RSO solar panel and the contribution due to the RSO bus. Following the assumptions that the solar panel tracks the Sun (possibly with an offset) and that the bus tracks nadir, the dynamics of a RSO can be approximated with a two-facet model. The panel's tracking of the Sun is due to its need for power regulation and this tracking is assumed to not be perfect. That is, we allow for an East-West offset angle in our calculations. One facet in the model represents the RSO solar panel and tracks the Sun. This facet is assumed to be approximately planar and is assumed to have both specular and Lambertian reflective properties. As the solar panel is assumed to be planar, it is decomposed into a singular albedo-Area product. The second facet of the model represents the RSO bus. The bus is assumed to be nadir pointing and is a complex threedimensional shape. We only model the nadir facing side of the bus since the respective normal vectors of the other faces, assuming a rectangular shaped bus,  are oriented 90 degrees from the nadir vector. This would cause the reflections from those faces to be minial relative to the observer near the nadir direction.
The bus is further assumed to exhibit primarily Lambertian reflection based on the following argument. Being a complex three dimensional object, there are a multitude of normal vectors from each of the many facets that point in a variety of directions. Each facet has its own BRDF that governs how much light is reflected off-normal. Since the body is spatially unresolved, the aggregate reflected light off the body is a superposition of the specular and diffuse reflections from all of these facets. To first order, the resulting bundle of light is Lambertian mainly due to the mixture of BRDFs and mixture of normal vector angles.
The reflective properties are governed by the geometry of the Sun and the observer and therefore pose-dependent. Finally, it is important to note that the algorithm's accuracy is dependent on proper identification of the specular regions of the signature. If the specular reflection from a panel, or dish, are incorrectly attributed to the body, then the specular reflectance will introduce errors into the Body aA calculation. Due to the pose-dependence, the RSO bus is decomposed into a number of different albedo-Area products to accommodate differential albedo-Areas under different observation conditions.
At present the method for solving for the various albedo-Area products of the bus is reviewed. A method for integrating the various albedo-Area products into a single average albedo-Area product for the RSO bus is additionally discussed. The method for solving for the albedo-Area products of the RSO bus leverages the assumed Lambertian reflectance of the bus to represent the albedo-Area products in terms of a system of linear equations which can be efficiently solved. In contrast to the bus, the specular nature of the solar panel's reflectance requires the reflectance be modeled in terms of a function that has high reflective properties near the point of ideal specular reflection and low reflective properties further away. The specular nature of the Solar Panel's reflectance requires the usage of a non-linear shape parameter. Due to this introduced non-linearity, the current method cannot be applied to the solar panel. An alternative method for computing the albedo-Area of the solar panel is later discussed in "The Solar Panel (SP) Method" section. As such, the current discussion concerns the decomposition of the body contribution to a signature into a number of posedependent albedo-Areas.
The RSO bus is modeled with a Lambertian basis function. That is, the projected albedo-Area of the RSO bus is proportional to cos(θ s ) cos(θ o ) where θ s is the angle between the RSO-nadir vectors and the Sun, and θ o is the angle between the RSO-nadir vector and the observer. As the observed brightness of the bus is assumed to be pose dependent, the bus is modeled with multiple basis functions, one per orbit angle bin, multiplied by participation factors w i (γ) which are functions of the orbit angle. The orbit angle (OA) is the angle between the projection of RSO-Sun vector into the orbital plane and the RSO-nadir vector as seen in Fig. 3. The participation factors are weights modulating the contribution each observation provides to its respective bins. The eleven bins are 30 degree regions spaced 15 degrees apart throughout the range [−75, 75] degrees orbit angle (the first and last bin are cut off and thus are only 15 degrees wide). This means the bins overlap across ten 15-degree intervals, and every point an interval will fall within exactly two bins. The particular choices of participation factors w i are somewhat arbitrary; however, here they are assumed to be piecewise continuous linear functions. For the bus, then, the total projected albedo-Area g aA BUS is a sum of the projected albedo-Areas across all orbit angle ranges: aA i is the in-orbit angle albedo-Area of the RSO bus, and k is the total number of orbit angle bins. The above equation is linear in terms of unknown coefficients aA i . Thus the parameters can be solved for by constructing a linear system of equations: X is a design matrix of basis functions per orbit angle bin: ð Þcos θ on ð Þw 1 γ n ð Þ cos θ sn ð Þcos θ on ð Þw n γ n ð Þ ⋯ cos θ sn ð Þcos θ on ð Þw k γ n ð Þ β is a column vector of the unknown aA i , and Y is a column vector of observed projected albedo-Areas. Projected albedo-Areas can be directly computed from photometric measurements as per [1]. The ordinary least squares solution can be computed as β = pinv(X)Y where pinv represents the Moore-Penrose pseudoinverse pinv(X) = ( X T X) −1 X T . However, using the pseudoinverse does not guarantee non-negativity of the β vector. To force non-negative values in β we have implemented an iterative approach which solves for a strictly non-negative least squares solution. This solving method is a variation on the work by Seung-Jean Kim, et al. [2]. Now that we have shown how to compute a vector of pose dependent albedo-Areas, we will discuss how to convert this vector into a singular albedo-Area value to represent the satellite as a whole. We wish to obtain a single aA value because it is easier for an analyst to work with one value rather than a vector of values. It is also easier to perform change detection using one aA value rather than the entire vector.
The method we use is as follows: We only use the lambertian observations with γ(orbit angle) values in [−75, 75] to calculate the body diffuse aA. There are 11 orbit angle bins, each spanning 30 (BinSize) degrees measured in γ. From the previous steps, we obtain β which holds 11 aA values. In what follows, the symbols ⌈ x⌉ and ⌊x ⌋ represent the ceiling and floor functions respectively: The subscript L below corresponds to the value for the bin to the left of the observation and the subscript R below corresponds to the value for the bin to the right of the observation.
For example, if γ i = − 53 then index iL = 1, w iL = 0.5333, index iR = 2, and w iR = 0.4667. β i is the value of the i th entry of β for i = 0, 1, …, 10. Each corresponding to one of the bins.
The Differential Method (DM) for Body aA Many of our methods treat a signature as the primary analysis source. In our usage, "signature" refers to a plot of brightness or color index vs. the longitudinal component of the total sensor/RSO/Sun angle (solar phase angle) (LPA). Figure 4 shows the vectors to the sensor (Observer), RO, and the Sun, RS, originating from the RSO. These two vectors are projected into the Earth's equatorial plane (dotted lines). The angle between the two projections is the LPA. LPA = 0 o is defined when the object lies directly opposite the Sun in its diurnal motion. Hence negative LPA values occur earlier in the night, and positive LPA values occur later. An experienced analyst can qualitatively extract a number of important features from a signature, and we are continuing development to automate these kinds of analyses. For example, a signature can provide information about the offset of the solar panels from their nominal direction of pointing at the Sun. It can also often provide the number and placement of large communication (com) antennas. Since different manufacturer/bus types tend to have self-consistent offsets and large structures, signature types are correlated with manufacturer/bus types.
As examples, we present Figs. 5 and 6, which illustrate simple and complex signature types, respectively. In both plots, we use the range corrected (36,000 km) Johnson R band magnitudes (i.e., R abs , to indicate brightness). We note that all of the DM analysis presented here will use R abs . The reason for concentrating on the R band observations is that we are trying to analyze the portions of the RSO that do not include  Fig. 5 shows only one brightness peak that is located near LPA = 0 o and has monotonically decreasing brightness on either side. In contrast, Fig. 6 presents the signature of SES-1 (bus type STAR 2), and one can see shoulders, which are located between 30 o and 40 o on either side of the central peak, which, in this case, is offset by about 8 o from the Sun (i.e., the solar panel normal is offset by about 4 o from the direction to the Sun). At some times during the year, as seasonal illumination conditions vary, these shoulders turn into broad peaks. Our interpretation of these shoulders and side peaks is that they arise from diffuse reflections from large com antennas, with the peak at negative LPA arising from an antenna located on the eastern or leading side of the satellite.
We have pointed out the degree of complexity of signatures because this is an important consideration in applying the DM method. We use the brightness information in the range abs(LPA) > 40 o in the body aA calculation. In the RSO-based coordinate system discussed in "Satellite-Based Coordinates" section, this is also approximately the same range in ΔΦ at these angles that lie far from 0 o . To emphasize these points, note that the brightness in the LPA > 40 o range is higher for the complex-signature satellites than that of the simple-signature satellites, and we allow for this natural difference in our DM method.
The Differential Method is based on determination of the mean difference in magnitude between the new observations of a satellite and a baseline catalog for that satellite's signature type. These catalogs were built up from observations taken during the GOLDS project, in which multi-site observations were made in 2013. They feature many exposures per night in B and R. We have constructed four baseline catalogs-that of the Single Peak, the Double Peak (Left), the Double Peak (Right), and the Triple Peak objects. Double Peak (Left) signatures have a single com antenna whose corresponding brightness peak is usually somewhat lower than the central peak, which is generally caused by specular or near specular reflections from the solar panels and are located typically in the range -10 o < LPA < +10 o . The lower, broad side peaks tend to be located in the range 30 o < abs(LPA) < 40 o . The Double Peak (Right) signatures are essentially mirror reflections of the Double Peak (Left) signatures. The Triple Peak signatures show evidence for peaks at both negative and positive LPAs. For the four baseline catalogs, we need a baseline body aA value. These body aAs were estimated using satellite models. Figure 7 presents the magnitude and body-referenced angle information found in the baseline catalogs. For illustrative purposes, we have added all of the catalogs together. There are two main groupings in the ΔΦ coordinate, and the noticeable tracks that lie at a nearly constant value of ΔΦ but move from front left to back right in θ s are the trajectories that the data takes when observations are taken many times throughout a single night. The two groupings (roughly half positive and half negative) in ΔΦ represent the positive and negative portions of the Phase Angle distribution, i.e., those observations made prior to the moment of 0 o LPA (the positive ΔΦ values) and those made after 0 o LPA.
Although the DM method is based on catalogs built up from moderate-cadence observations taken during a series of whole nights, any number of new observations can be used to estimate aA DM . Conceptually, one simply plots the new θ s , ΔΦ, and R abs observations onto a plot similar to that of Fig. 7 and then finds the mean distance in magnitude between the new data point and the mean R abs of the nearest catalog neighbors. This magnitude difference is converted into a flux ratio, which is proportional to the aA ratio between the newly observed satellite and the mean aA of the satellites used to construct the catalog. The equation that we use is aA = aA mean,i * (0.4*10^(-ΔR abs )). In this equation, the subscripts "mean" and "i" indicate the mean aA value adopted for the i th of the four catalogs, and ΔR abs is the average difference between the new data and the catalog mean. An example of the process is shown in Fig. 8. It illustrates most of the features of the DM method. The red plus sign markers indicate the distribution of the new data in the body-referenced coordinate system for the satellite AMC-18 on the night of Feb. 5, 2013 as seen from the Kirtland Raven. The blue circle markers indicate the mean value of the N = 25 nearest neighbors that are found in the catalog (neighbors in the θ s and ΔΦ coordinates). The new estimate of aA DM is directly calculated from the average value of the R abs difference between the catalog means and the new data.

Comparisons Between 2FM and DM Methods for Body aA
One problem in developing algorithms that calculate aA is there generally does not exist ground truth for this quantity. One can examine photographs and artists' conceptions in order to estimate the projected area of components when they are in orbit and then guess the albedo of the materials that appear to be important contributors of the total reflection. However, this is not very satisfactory, especially given the rather wide range of albedos for tested materials.
Given the lack of ground truth, we primarily rely on comparisons between our two methods in order to perform a low level of verification and validation. In addition to the fact that the two methods generally agree well, we add in the facts that the calculated values agree well with comparisons to web illustrations and also that the disagreements, when they do exist, agree with an analysis of the physical nature of the satellites' major components. Table 1 presents our comparisons for the eight primary target satellites in the GOLDS1 observing campaign, which took place in the late winter and spring of 2013. Column 1 lists the common satellite name (or obvious abbreviation). Column 2 provides comments about the complexity and structure of the photometric signatures for each satellite. The single peak and triple peak objects have been discussed above. However, the two satellites, AMC-15 and AMC-18, have signatures that are of a higher complexity. During the portion of the observing campaign that took place prior to the Vernal Equinox, they each showed Double (Left) signatures, while after the equinox, they had Double (Right) signatures. This behavior is quite unusual in our experience and suggests that the com antennas on these two satellites might be mounted in planes that are tilted with respect to the Celestial Equator. Column 3 lists our best understanding of the manufacturer/bus type. The A2100 buses and derivatives are made by Lockheed Martin. The SS/L objects are made by Space Systems Loral. The BSS satellites are made by Boeing Space Systems, and the STAR 2 is made by Orbital Sciences Corp. Columns 4 and 5 list the estimated aA values resulting from the DM and 2FM methods, respectively. The quoted uncertainties are standard error of the mean. All quantities are in units of m 2 .
We note the four Single Peak satellites show a high level of agreement between the two body aA methods of calculation. However, the four satellites with complex signatures show much less agreement, and we see the 2FM method produces higher values of aA than the DM method does for each of these four RSOs.

The Solar Panel (SP) Method
Here we discuss and derive our method for the calculation of the Solar Panel albedo-Area (SPaA). This section is broken up into a preliminary and an algorithm section. In the preliminary section, we discuss the geometry of the problem and give details on various sub-algorithms needed for calculation for SPaA. In the algorithm section, we discuss the step-by-step process from taking a light signature and returning a value for SPaA.

The Geometry of The Solar Panel and Its Albedo
We must consider the geometry involved with the solar panel and the incident light relative to the observer. Refer to Fig. 9 for the geometry discussed. We denote the solar panel's normal vector asN. From this vector, the angle towards the Sun's incident light is labeled θ i , and the angle towards the observer's location is denoted θ r . The angles ϕ i and ϕ r are defined as the angles, in the panel's plane, going from the x − axis to the projection of incident and reflected light into the panel's plane. We assume the panel's Bidirectional Reflectance Distribution Function (BRDF) to have a Gaussian shape as a function of off-reflection angle, as suggested in [3], under this assumption the panel's reflected light will have a rotational symmetry. This is the same geometric convention used by [4]. Note that we use θ i here in reference to the incident light to the panel's surface whereas in the body-frame we use θ s . The BRDF is defined to be the function f(θ i , θ r , ϕ i , ϕ r ) which describes the scattering of radiance, reflected off of the panel's surface, in all directions above the hemisphere. This is written as: where L r is the spectral radiance of the reflected light, and E is the spectral irradiance. The albedo is the total (unitless) ratio of reflected to incident light off of a given surface throughout the entire hemisphere of possible reflections. This quantity, denoted a, is defined as: This albedo formula (6), as discussed in [4,5], is derived in detail in Appendix 1.
Another angle used to simplify the notation is the glint angle θ * or off-reflection angle (see Fig. 9). This is the angle from the observer to the reflection of the incident light across the panel's normal. This angle can be written as: The derivation of the off-reflection angle formula (7) is given in detail in Appendix 2. Calculating off-reflection angle from a solar panel offset angle Another formulation of the off-reflection angle which is needed comes from simple geometry. We assume that our signature data will be collected in terms of longitudinal phase angle, and if we also know the solar panel's correct offset angle, then we can convert our signature into brightness versus off-reflection angle (from phase angle). To see this, we define the off-reflection angle as the angle between the panel's specular reflection vector, P ! spec , and the panel to observer vector, RO ! (RSO to observer vector).
Assuming these vectors are normalized, then, an equivalent form of (7) is: RO ! is the difference between the position vectors of the observer and satellite with respect to the Earth-centered frame, the panel's specular reflection vector can be calculated by: The RS ! vector is the difference between the position vectors of the Sun and satellite, and the panel's normal vector can be found by fixing some solar panel offset angle, SPO, and calculating: where d Orb is the orbital north vector defined as d RSO ÂV. d RSO andV are the normalized vectors of satellite position and velocity, respectively. Hence, if we have the solar panel's offset angle, then Equations (10), (9), and (8), in that order, give us a conversion from phase angle to off-reflection angle.

Derivation of Solar Panel aA
The solar panel specular albedo-Area is defined as the product of its albedo and its area, thus: with a being the albedo and A the area of the panel in m 2 . Let g(θ * ) be the shape of the specular BRDF. We assume the panel's BRDF to have a Gaussian shape as a function of off-reflection angle, as suggested in [3], given by: where σ is the width parameter of the curve and is determined by the linear fitting method (LFM) described in "The Linear Fitting Method (LFM)" section . The use of LFM is done in the log space where the function g(θ * ) becomes linear after a trigonometric change of variable. To be used as the BRDF's shape, g(θ * ) must be normalized as to allow the BRDF to give the distribution of albedo over the hemisphere of possible observer locations. To achieve this, we define a normalizing factor G norm . This constant will be the value by which the volume under g(θ * ) is normalized to 1 over the hemisphere: This allows for us to define our Gaussian BRDF for the specular case by: Equation (14) denotes the scaled albedo apparent to the observer at all off-reflection angles in the viewable hemisphere. Next we consider the model we are using for the panel's reflection, which consists of both specular and diffuse parts as per [5]. We estimate that the solar panel radiant intensity, measured in W sr , would be composed of a primarily specular pattern (at least 80%), while the remaining intensity would be due to a diffuse contribution (less than 20%). This resulted in our model: where E is the in-band Solar irradiance (E≈642 W m 2 for open band). Since we are using inverse modeling, we attempt to capture the specular component of the solar panel in its entirety. This method results in the weighting constants being modified so that the panel has an entirely specular nature (C diff is neglected). That is, our model in Eq. (14) is simplified to become: Using a light curve, we apply the LFM to the observed data to get: for some scaling constant P, which represents the peak of the specular contribution whose shape is assumed to be Gaussian. Equation (17) is extracted from the light curve data, then we compare the fitted values to our theoretical model. By equating the fit to the model, S fit (θ * ) = S model (θ * ), we get: which by cancellation says that which is then easily solved for the solar panel albedo-Area, units of (aA) panel are m 2 : The Linear Fitting Method (LFM) This algorithm was developed in order to better estimate the specular region of the light curve in longitudinal phase angle. We approximate the specular shape by the function f θ * À Á ¼ Pg θ * À Á ¼ e a*cos θ * ð Þþb where the light curve is radiant intensity (RI), measured in w sr , as a function of off-reflection angle, θ * . Using this model has a few advantages. Note then that the natural log of radiant intensity is linear with respect to the cosine of off-reflection angle. We will exploit this relationship and call the algorithm the linear fitting method (LFM). For the LFM, we assume a Gaussian BRDF as in Equation (12) which can be written as f θ * À Á ¼ Pg θ * À Á ¼ e a*cos θ * ð Þþb . While the exponent containing the cosine term is equivalent to the sine squared version through a trigonometric identity shown in Equation (20), we opt for the cosine form. Also note, this model is symmetric, so the sign of the off-reflection angle is irrelevant (the solar panel specular contribution to a signature is symmetric with respect to θ * ). This means it is easy to tell where the signature is symmetric; it is symmetric where two halves of the signature, one with positive off-reflection angles and the other with negative offreflection angles, fall on top of one another. Where in the light curve the Gaussian assumption is true, to the plot of log(RI) vs cos(θ * ) is linear. This region is where the reflectance is specular. Note that, although the specular region is generally contained between −20 and 20 degrees phase angle, the LFM uses the entire range of phase angles to ensure a robust detection of the entire specular region. For example, this allows for inclusion of specular reflectance outside of this range, which the authors have seen in real data.
Using this new model, we perform an iterative search for where the two halves of the data are symmetric and linear. The specular region is centered around the off-reflection angle θ * = 0, so at each iteration we remove the same fraction, say 1/4th, of the data with the lowest off-reflection angle and process the data through LFM again. We stop once we find a subset of observations that is sufficiently linear; checking if the data is linear incorporates a check for if the data is symmetric because if two halves do not fall on top of each other in the log(RI) vs cos(θ * ) space, then the data cannot be linear. In order to check if the data is linear, we fit the data with the line produced by the observation points with the largest and smallest off-reflection angle. If any of the remaining points are further from the line than some predetermined distance δ, then the data is not sufficiently linear, and the next iteration of the algorithm is executed. The value of δ was determined by examining many candidate values and choosing the one which gave the best results for typical light curves. The value of δ used is 0.001. Once the data remaining in the LFM's active region satisfies our maximum distance bound, δ, we use the line through the observations with the largest and smallest off-reflection angles as the specular fit for the data. The method used to define the line here is done to ensure that the point at the peak will be included in the fit which is not guaranteed by a least squares procedure. In experimentation it was found that the results were very sensitive to the inclusion of the peak observation. A flow chart of the algorithm can be seen in Fig. 10. Figure 11 shows an example run of the algorithm. The light curve is plotted by taking the natural logarithm of brightness as a function of cos(θ * ) and is approximately linear in the specular region. Each red line represents one iteration of LFM. This process occurs from left to right, so the first iteration is represented by the left most line, the second by the line directly to the right of that, and so on. For a given iteration, all of the data with smaller off-reflection angles (larger cosine values) than the line is used for the current iteration. The data to the right of the rightmost line represents the observation that were found to be linear.  Fig. 10 Flowchart of the linear fitting algorithm Next consider the natural log This is the equation of a line in the log-cosine space and it follows that m ¼ 1 σ 2 and b ¼ ln P ð Þ− 1 σ 2 . These relationships provide us with σ 2 ¼ 1 m and P = e m + b .

The Solar Panel Method Algorithm
The process of generating the solar panel albedo-Area takes three stages to complete. First, a solar panel offset angle must be correctly identified. This is the angle between the solar panel's normal vector and the Sun. Second, using the correct offset angle allows for the light curve to be converted from longitudinal phase angle to off-reflection angle. Finally, in the off-reflection angle space, the LFM allows for Gaussian parameters to be extracted from the data. These Gaussian parameters are directly applied to calculate the SPaA. The overall process occurs in the following manner with regard to a single light curve which we assume to be data collected in radiant intensity, RI, in units of W sr , vs longitudinal phase angle:

Solar Panel aA Result Summary for Usno Dataset
Here we present a summary of our results on the solar panel albedo-Areas found by the method presented. The photometric observation data was collected by sensors over the course of two years starting in 2013. The mean aA values shown in Table 2  are averages over all light curves returning a value as some nights do not contain enough data to produce a meaningful value. As part of our verification process, the aA values were compared with satellite models with generally good results. Over the spectral range of panchromatic data, the "true" albedo of solar panel materials varies by a large amount; it also varies significantly between the material types of Silicon and Gallium Arsenide. Averaging over the wavelength, the panchromatic albedo is about 10%. Therefore our aA results will almost always be significantly smaller than the actual size in m 2 . The range of values for SPaA were also large, since the reflectivity of the solar panels varies a great deal with season. The observation angles during the Equinoxes are much more closely aligned with the glint angle of the panels than during the Solstice periods of the year. Therefore, our generally good results were found from Equinocial light curves with poor results for light curves collected during the Summer and Winter months.

Discussion and Conclusion
We have presented results on one new algorithm (DM) and one algorithm under continuing development (2FM) for extracting the albedo-Area of the satellite body alone. For those satellites with simple signatures (i.e., one central brightness peak located near LPA = 0 o ), the agreement between the two methods is quite good. However, for those objects with more complex signatures (i.e., those with broad shoulders or secondary peaks that are located in the general range of 20 o < abs(LPA) < 40 o ), we find that the two methods disagree and that the 2FM method computes consistently larger body aAs than the DM method does. There is an explanation for the disagreements. First, the DM method excludes brightness data that lies in the range abs(LPA) < 40 o , whereas the 2FM method often includes brightness data that lies in the range 20 o < abs(LPA) < 40 o . The inclusion of data that is closer to the central LPA regime does not affect the single peak signature objects because the bright near-specular solar panel reflections do not encroach on the LPA range of either method in such signatures. However, for those satellites that show side peaks or shoulders, the "excess" brightness does often fall in the LPA range of the 2FM method (cf. Figs. 5 and 6). These differences do not strike us as flaws in either method. Indeed, we believe that it might be possible to exploit the differences in a manner that allows for a richer set of features that can be gleaned and thereby provide for better SSA analyses. The differences between the DM and 2FM methods, to the best of our current interpretation of the data, lie in whether or not a substantial portion of the reflected light of the com antennas are included in the aA calculation. Hence the DM estimates of aA appear to be more directly related to characterize the size of the main body/payload of the satellite, while the 2FM estimates of aA appear to be more directly related to the body/payload plus the com antennas. We are currently evaluating whether or not these differences can be fully exploited.
Another new algorithm has been developed and tested for estimating the albedo-Area of the solar panels based on the specular contribution of the overall reflected intensity from the satellite. The diffuse contribution from solar panels, by their design and function, is negligibly low and so is ignored in this approach. The sizes of the solar panels vary with satellite design and type of solar panel material. The areas can range from 12 to 70 m 2 for the satellites we have studied. Since albedo is less than one, reflectance aAs are only proxies for the true physical size. Further research might indicate if the aAs can be studied over time to determine effects of space weather on the reflectivity of the materials.
The accuracy of these approaches for estimating the component aAs is dependent on sampling the various parts of the satellite well. Research is ongoing to understand the effect of sparse data on the resulting aAs. Another challenge in validating these approaches is the availability of truth data. Our research group has been using satellite model dimensions as a proxy for truth data, which is the most accurate comparison since specific albedos are unknown. Even measured albedos in databases can be inaccurate since the albedos can change with different manufacturers and with space aging.
One future area of work is to investigate the use of splines for improving the estimate of the specular regions of the signatures. This would improve both the solar panel offset calculation and the solar panel aA. Another area of work will be to investigate how to fuse the DM and 2FM methods for improved body aA estimation with sparse data and possibly to quantify the contribution of the antennas.