.--------------------~- -- - --_.-

NASA/TM-2012-216308

A Probabilistic, Facility-Centric Approach to

Lightning Strike Location·

Author Name: Lisa L. Huddleston

Kennedy Space Center, Florida

William P. Roeder

45/h Weather Squadron, Patrick Air Force Base

Francis J. Merceret

Kennedy Space Center, Florida

January 2012
NASA STI Program ... in Profile

Since its founding, NASA has been dedicated

to the advancement of aeronautics and space

science. The NASA scientific and technical

infonnation (ST!) program plays a key part in

helping NASA maintain this important role.

The NASA STI program operates under the

auspices ofthe Agency Chief Infonnation

Officer. It collects, organizes, provides for

archiving, and disseminates NASA's STI. The

NASA STI program provides access to the

NASA Aeronautics and Space Database and its

public interface, the NASA Technical Report

Server, thus providing one of the largest

collections of aeronautical and space science

STI in the world. Results are published in both

non-NASA channels and by NASA in the

NASA STI Report Series, which includes the

following report types:

• TECHNICAL PUBLICATION. Reports of

completed research or a major significant

phase ofresearch that present the results of

NASA Programs and include extensive data

or theoretical analysis. Includes compilations ofsignificant scientific and technical

data and infonnation deemed to be of

continuing reference value. NASA counterpart of peer-reviewed fonnal professional

papers but has less stringent limitations on

manuscript length and extent of graphic

presentations.

• TECHNICAL MEMORANDUM.

Scientific and technical findings that are

preliminary or ofspecialized interest,

e.g., quick release reports, working

papers, and bibliographies that contain

minimal annotation. Does not contain

extensive analysis.

• CONTRACTOR REPORT. Scientific and

technical findings by NASA-sponsored

contractors and grantees.

• CONFERENCE PUBLICATION.

Collected papers from scientific and

technical conferences, symposia,

seminars, or other meetings sponsored

or co-sponsored by NASA.

• SPECIAL PUBLICATION. Scientific,

technical, or historical infonnation from

NASA programs, projects, and missions,

often concerned with subjects having

substantial public interest.

• TECHNICAL TRANSLATION.

English-language translations of foreign

scientific and technical material pertinent to

NASA's mission.

Specialized services also include organizing

and publishing research results, distributing

specialized research announcements and feeds,

providing help desk and personal search

support, and enabling data exchange services.

For more infonnation about the NASA STI

program, see the following:

• Access the NASA STI program home page

at http://www.sti.nasa.gov

• E-mail your question via the Internet to

help@sti.nasa.gov

• Fax your question to the NASA STI Help

Desk at 443-757-5803

• Phone the NASA STI Help Desk at

443-757-5802

• Write to:

NASA STI Help Desk

NASA Center for AeroSpace Infonnation

7115 Standard Drive

Hanover, MD 21076-1320
NASA/TM-2012-216308

' .. .,...

- : . . :.. ". .:' ...:" . .

,

A Probabilistic, Facility-Centric Approach to

Lightning Strike Location

Author Name: Lisa L. Huddleston

Kennedy Space Center, Florida

William P. Roeder

45th Weather Squadron, PatrickAir Force Base

Francis J Merceret

Kennedy Space Center, Florida

National Aeronautics and

Space Administration

Kennedy Space Center

Kennedy Space Center, FL 32899-0001

January 2012
Acknowledgments

This work would not exist were it not for the generous contributions and suggestions by Dr. Ken Chan of

the Aerospace Corp, who is an expert on spacecraft collision probability, on whose work the probability

of lightning within a radius of interest was based. Dr. Damn Leleux ofthe Johnson Space Center,

Houston, TX, and Dr. Walt Gill of Sandia National Laboratory also provided helpful guidance and

testing. Mr. Jeremy Hinkley and Mr. Pete Hopman ofUnited Space Alliance, the prime contractor for

Space Shuttle operations, provided very valuable efficiency modifications for the visual basic code and

also integrated the probability calculations and closest ellipse point algorithm into the 45th Weather

Squadron Lightning Report Spreadsheet.

The authors appreciate a helpful review of an earlier draft ofthis paper by Mr. John Madura ofthe

Kennedy Space Center. This work was done under the Kennedy Space Center Employee Development

Program.

Available from:

NASA Center for AeroSpace Information

7115 Standard Drive

Hanover, MD 21076-1320

443-757-5802
Table of Contents

List of Figures 2

List of Tables 3

Abstract 4

Nomenclature 5

Introduction 6

Methodology 6

Background 6

The Numerical Integration Technique (Chan, 2011) 7

Evaluation 9

Test Set 1 9

Test Set2 10

Other Applications 15

Conclusion 15

Appendix A - Example Lightning Probability Calculation 16

Appendix B - Excel VBA Macro for 45th Weather Squadron Lightning Spreadsheet that Calculates the

Probability of any Nearby Lightning Stroke being Inside any Radius ofAny Point ofInterest.. 20

Conversions Module 20

Functions Module 21

Module 2 Code 25

Probability Module 28

Strike Class Module 31

References 39
List of Figures

Figure 1. Schematic diagram ofthe angles used in probability calculation for a sample lightning location

error ellipse 7

Figure 2. Sample oflightning strikes where the closest point on the lightning position uncertainty ellipse

was within 0.45 nmi of Launch Complex 39A on 3 August 2009 10

Figure 3. Google Maps visualization ofthe 99% confidence uncertainty ellipse for one ofthe closest

lightning strikes to Complex 39A on 03 August 2009 12

Figure 4. Google Maps visualization ofthe 99% confidence uncertainty ellipse for a lightning strike near

Complex 39A on 03 August 2009 12

Figure 5. Google Maps visualization ofthe 99% confidence uncertainty ellipse for nearby lightning strike

to Complex 39A on 03 August 2009 13

Figure 6. Illustrates a probability of 69.1% of a lightning strike of amplitude -43.0 kA detected by NLDN

occurring 0.26 nmi from the center of Launch Complex 39A on 8/16/2009 14

Figure 7. Illustrates a probability of 74.7% of a lightning strike of amplitude -71.4 kA detected by NLDN

occurring 0.28 nautical miles from the center of Launch.Complex 39A on 10/14/2009 14

Figure 8. Illustrates a probability of 99.9996% of a lightning strike of amplitude -21.7 kA detected by

CGLSS occurring 0.04 nmi from the center of Launch Complex 39B on 6/27/2009 15

2
List of Tables

Table 1. Calculated probability vs. CRC Handbook probability for various inputs 10

Table 2. Input values used for scenarios shown in Figures 3 through 5 11

Table 3. Input values used for scenarios shown in Figures 6 through 8 13

Table 4. Lightning strike probability calculation process 16

3
Abstract

A new probabilistic facility-centric approach to lightning strike location has been developed.

This process uses the bivariate Gaussian distribution of probability density provided by the current

lightning location error ellipse for the most likely location of a lightning stroke and integrates it to

determine the probability that the stroke is inside any specified radius of any location, even ifthat location

is not centered on or even with the location error ellipse. This technique is adapted from a method of

calculating the probability of debris collision~ith spacecraft. Such a technique is important in spaceport

processing activities because it allows engineers to quantify the risk ofinduced current damage to critical

electronics due to nearby lightning strokes. This technique was tested extensively and is now in use by

space launch organizations at Kennedy Space Center and Cape Canaveral Air Force Station. Future

applications could include forensic meteorology.

4
Nomenclature

e angle between the collision plane coordinates and the rotated collision plane coordinates

pxz correlation coefficient of x and z

(JK' standard deviation of x-coordinate ofthe diagonalized covariance ellipse in the rotated

coordinate system

(JH' standard deviation of z-coordinate ofthe diagonalized covariance ellipse in the rotated

coordinate system

(Jx standard deviation of x

(Jz standard deviation of z

(Jx' standard deviation ofx in the rotated coordinate system

(Jz, standard deviation of z in the rotated coordinate system

11K x-coordinate oftarget circle in the (X', Z') coordinate system

flH z-coordinate oftarget circle in the (X', Z') coordinate system

A collision cross-sectional area (nautical miles2

, nmi2)

de the angle between two points on the target ellipse

dH integration step

H intermediate variable in the lightning probability algorithm

P probability

pdf probability distribution function

rA radius of circle of cross-sectional area, A (nautical miles, nmi)

R radius of ellipse (nmi)

Rl distance to first point on target ellipse (nmi)

R2 distance to second point on target ellipse (nmi)

x horizontal rectangular coordinate in the collision plane (nautical miles, nmi)

x' horizontal rectangular coordinate in the rotated collision plane (nmi)

x" transformation variable that circularizes the x-component ofthe probability ellipse (nmi)

Xe nominal distance of closed approach oftwo colliding objects (nmi)

W intermediate variable in the lightning probability algorithm

z vertical rectangular coordinate in the collision plane (nautical miles, nmi)

z' vertical rectangular coordinate in the collision plane (nautical miles, nmi)

z" transformation variable that circularizes the z-component ofthe probability ellipse (nmi)

5
Introduction

The ability to a Introduction accurately estimate the probability that an individual nearby cloudto-ground lightning stroke was within a specified distance of any specified spaceport processing facility at

Kennedy Space Center (KSC) or Cape Canaveral Air Force Station (CCAFS) is important to processing

payloads and space launch vehicles before launch. Such estimates allow engineers to decide ifinspection

of electronics systems aboard satellite payloads, space launch vehicles, and ground support equipment is

warranted due to indu~ed currents from that stroke. Ifinduced current damage has occurred, inspections

ofthe electronics are critical to identify required fixes and avoid degraded performance or failure ofthe

satellite or space launch vehicle. However, inspections are costly both financially and in terms of delayed

processing for space launch activities. As such, it is important these inspections be avoided if not needed.

At KSC/CCAFS, one ofthe main purposes ofthe Four Dimensional Lightning Surveillance System

(4DLSS) (Murphy, 2008, Roeder, 2010) is detection of nearby strokes and determination oftheir peak

current to support decisions to inspect electronics (Flinn, 201 Oa, Flinn, 201 Ob, Roeder, 2005). The high

frequency oflightning occurrence in East Central Florida combined with the large amount of complex

sensitive electronics in satellite payloads, space launch vehicles, and associated facilities makes those

decisions critically important to space launch processing. The 4DLSS provides the data for 50th

percentile location error ellipses for the best location for each stroke, which is then scaled to 95th or 99th

percentile ellipses depending on customer requirements. This error ellipse is necessarily centered on the

best location ofthe lightning stroke. The 4DLSS, however, has not been able to provide the probability

ofthe stroke being within a customer specified distance of a point ofinterest. This paper presents a new

method to convert the 4DLSS 50th percentile location error ellipse for best location of any stroke into the

probability that the stroke was within any radius of any facility at CCAFSIKSC. This technique could be

adapted for use with National Lightning Detection Network (NLDN) data. This new probabilistic facilitycentric technique is a significant improvement over the stroke-centric location error ellipses the 45th

Weather Squadron (45WS) has provided in the past. This technique is adapted from a method of

calculating the probability of debris collision with spacecraft (Chan, 2008, Leleux, 2002, Patera, 2001).

Methodology

Background

In spacecraft collision probability and other applications, at the ins.tant of "nominal" closest

approach, the position uncertainty ofthe collision object relative to the asset being protected is described

by a bivariate Gaussian probability density function (pdf) (Chan, 2008, Patera, 2001, Alfano, 2006,

Alfano 2007), as shown in the following equation

(1)

where O'x and O'z = the standard deviations of x and z, pxz = correlation coefficient of x and z, x and z are

the designations for the rectangular coordinates in the collision plane.

The probability of collision (Eq. 2) is given by the two-dimensional integral, where A is the

collision cross-sectional area which is a circle with radius, rA (Chan, 2008).

6
---------------- ----

p= fff(x,z)dxdz

A

(2)

There is no known analytical solution to the above integral when the two standard deviations ax

and az are not equal. The solution is found by performing a numerical integration ofthe two dimensional

Gaussian pdf (Chan, 2008, Patera, 2001, Alfano, 2006, Alfano 2007).

The geometry used for spaceflight collision probability can also be used for estimation ofthe

probability of an individual nearby lightning stroke contacting the surface within a specified distance of a

specified point ofinterest as shown in Fig. 1. In Fig. 1, a is the heading ofthe semi-major axis ofthe

lightning location uncertainty ellipse from true north and eis the angle between the semi-major axis of

the lightning location uncertainty ellipse and line connecting the center ofthe lightning uncertainty ellipse

and the center ofthe area ofinterest. Two methods ofintegrating the above probability were tested and

gave identical results. The first solution method was based on an algorithm by Patera (2001) and the

second solution method was based on an algorithm by Chan (2011). Chan's algorithm ran much faster

and therefore was selected as the algorithm for the 45WS lightning probability program.

1-0 uncertaintyellipse

around lightningstrike

Area around point of interest

__ @..,LOllPI /

---- a ....t-__ e //

,/\,

unple

The Numerical Integration Technique (Chan, 2011)

This numerical integration technique is one in which the miss distance is given by a non-central

chi distribution with unequal variances (Chan, 2011). The covariance matrix corresponding to the

bivariate Gaussian pdfin Eq. 1 is (Chan, 2008):

7
~-------------

(3)

When the correlation coefficient, Pxz, is not zero, there are undesirable off-diagonal tenus that

overly complicate the calculation. In order to eliminate these tenus, the coordinate system (x,z) is rotated

to a new coordinate system (x',z') such that the major and minor axes ofthe ellipse associated with the

covariance are aligned along the coordinate axes and the new covariance matrix is (Chan, 2008):

(4)

The angle, e, between the two coordinate systems is (Chan, 2008):

(5)

The KSC/CCAFS 4DLSS system does not provide the covariance matrix, but instead provides

the semi-major axis, semi-minor axis, and the orientation ofthe semi-major axis ofthe 50% location error

ellipse relative to north. Therefore the angle, e, in Eq. 5 is found using geometry where eis the angle

between the semi-major axis ofthe lightning location uncertainty ellipse and line connecting the center of

the lightning uncertainty ellipse and the center ofthe area ofinterest.

In the (x', z') system, the Eq. 1 pdf becomes (Chan, 2008)

(6)

and Eq. 2, the collision probability becomes (Chan, 2008)

where

A'= A, rA' =rA' x'p =xe cos (), ZIp =xe sin ()

(7)

(8)

For spacecraft collision, Xe is the nominal distance of closest approach ofthe two colliding objects

and (x'p, z'p) are the coordinates ofthe spacecraft relative to the debris. For lightning strike probability,

Xe (the distance between the position ofthe center ofthe strike location ellipse and the position ofthe

target area) is calculated using the Haversine distance fonuula.

8
The standard deviations in the new rotated coordinate system are calculated by dividing the semimajor and semi-minor axes ofthe 50% lightning positional confidence ellipse by.the scaling constant

used to scale standard error to the 50% confidence level. The scaling constant is:

The probability is given by (Chan, 2011)

k = ~- 2 *In(1- 0.50) (9)

where

(10)

(11)

The parameters ILK and IlH are the coordinates ofthe target circle in the (X', Z') coordinate

system; and (jK and a H are the standard deviations ofthe diagonalized covariance ellipse shown in Eq.

4. The derivation of equations (10) and (11) above is shown in further detail in (Chan, 2011). A detailed

example ofthe calculations using a real-world case is provided in Appendix-A. The Excel Visual Basic

code to implement these calculations is shown in Appendix-B.

Evaluation

The probability that any lightning strike is within any radius of any point ofinterest would be

extremely difficult to estimate intuitively. As a result, given the high impact ofthe decisions on space

launch operations, the tool developed for this application was extensively tested. Tests were conducted

and are discussed in the following sections: 1) known mathematical solutions, and 2) examination of

. real-world events. The new technique passed all ofthe tests. Tests were also conducted to assure the

probabilities calculated using the algorithm ofChan (2011) matched probabilities calculated using the

algorithm ofPatera (2001). The probabilities calculated from the two algorithms were identical and thus

are not shown here.

Test Set 1

The first set oftesting compared the lightning strike probability calculated using the 45WS

lightning strike spreadsheet (which uses an adaptation ofthe numerical integration algorithm by Chan,

2011, to the corresponding circular probability from the CRC Handbook of Tables for Probability and

Statistics (Beyer, 1968). Table 1 shows the probability from the new numerical integration technique for

various inputs and the corresponding correct probability from the CRC Handbook. The values matched to

within a tenth of a percent. These errors in the fmal digit may be due to round-off error.

9
----------------------

Table 1. Calculated probability vs. CRC Handbook probability for various inputs

Semi- Semi- Heading Point Point Strike Strike Radius Calcu- CRC

major mmor ofsemi- Of Of Latitude Long- around lated Handaxis axis major Interest Interest itude Point Of prob- book

(nrni) (nrni) axis from latitude long- Interest ability probtrue itude (nrni) ability

North (Beyer,

1968)

3 3 15 28.6082 -80.6041 28.6995 -80.6041 3 0.095 0.095

3 3 15 28.6082 -80.6041 28.631 -80.6041 3 0.453 0.452

3 3 15 28.6082 -80.6041 28.608 -80.6041 3 0.500 0.499

1 1 15 28.6082 -80.6041 28.608 -80.6041 1 0.500 0.499

1 1 15 28.6082 -80.6041 28.631 -80.6041 1 0.200 0.200

1 1 15 28.6082 -80.6041 28.6995 -80.6041 1 0.000 0.000

1 1 15 28.6082 -80.6041 28.608 -80.6041 2 0.937 0.938

Test Set 2

The second type oftesting analyzed six real-world lightning strikes near Space Launch Complex

39A on 3 August 2009. Figure 2 shows the spreadsheet used to generate the lightning report for those six

strikes. Additional data on three ofthese six strikes are in Table 2. These strikes were selected because the

closest point on the lightning position uncertainty ellipse was within 0.45 nautical miles of Launch

Complex 39A, which is the key radius for assessing the need to inspect electronics for induced current

damage to the Space Shuttle. Figures 3 through 5 are Google Maps depictions ofthree ofthese six

strokes. In Figures 3-5, the black ellipse is the 99% lightning location uncertainty ellipse. The white

circle is a 0.45 nrni radius around the point of interest. The probabilities for a small area around a facility,

even for a nearby stroke, may appear to be surprisingly low. For example, figures 3 and 4 respectively

illustrate a 53.8 % and 7.7% probability that the lightning strike occurred within the area ofinterest

while figure 5 shows that a strike just 0.65 nautical miles away had only a 1.1 % probability of being

within the 0.45 nautical mile radius of Launch Complex 39A. All calculated probabilities are consistent

with these real-world events.

Microsoft Excel

Dati Rf'ric!w \/Itw OlNelOpff Acfd-lns

" x

OUtiln.

1'_ t fit t't [.. Consolkiatt Whlt--lf

.",r;tnfv\lpb<.&teSVa 'I'" ~Is·

DItITools

'i/AcMnud

Sort It Filttf

'I ,.. "'t, f of Rdrnh ."wu )n AJ .. (,. f l.

Get Ellltm.1 Dati ConnectlOl'l1

t§\1.l

',~ HOlM Insert Page layout formulas

Q13495

CG IJghtnlng large RadII of tnterfst correetlon • 2S Oct lO.xlsm

A IS C 0 E f M

"

_ co x

~Mllp I ClO9" ,

.~

5

s~

• .v

19'1 01/05/2009 2l:llt.40.527 '50 0."

'" 0.00 -'UNo 346.07Z5~ 0." 2&.6107 ·80.6114 0.'2 0.1' ""-'" 21.607!J -1O.60S6

1949 0&/03/2000 21:10::40.601 "'" 0." 000 0.12 .",- 2S3.746~ 0.50 21.6U4 ·1O.6U3 056 0.26 500.71) 21.60&6 -10.6019

lJS' 00103/_ lU0:41.240 ,., 0.61 '10 0.20 "9.1 No 412.147405'6 0.41 21.6lU -80.6147 0.42 O.lA 301.50 'UO" -80.607'

2J7I 0&/05/2009 2.1:23:04.062 ... 0.59 '" 0.27 -25.2 No , 7.695324" 0.32 2&.6178 -10.15069 OA' 0.26 ZUOO la.6ll7 -10.600

5160 01/03/2009 21:25:41.710 216 0.65 m 0." -24.6 No 5 1.1lS743" 0.27 21.S99S -10.6113 0.26 O.lA 20.50 21.6036 -10.6090 ,,.. 00103/_ 2U7:1U21 222 0.56 B' 0.32 -7.5 No 4 U6622.OK 0." 28.6014 -lUlU 0.21 014 "90 21.6053 -80.6091

,.." c=::J

Figure 2. Sample oflightning strikes where the closest point on the lightning position uncertainty

ellipse was within 0.45 nrni of Launch Complex 39A on 3 August 2009.

10
_._-~ -------------------------------

The KSC Electromagnetic Environmental Effects (EEE) Panel requested six more real-world

lightning strikes be investigated. These were recently investigated lightning strikes near Launch

Complexes 39A or 39B where there was camera verification ofthe location ofthe strike. The EEE Panel

wanted to compare the results ofthe new facility-centric probabilistic technique to these cases where the

true answers were known unambiguously. The data used for this analysis are in Table 3. Both 4DLSS

and National Lightning Data Network (NLDN) cases were examined, depending upon which sensor

system recorded the stroke. CGLSS strokes were obtained from 45WS 4DLSS. The NLDN usually

provided flash data, so NLDN return stroke data were purchased as special StrikeNet1 reports from

Vaisala Corporation. This was done to match the return strokes routinely provided by 4DLSS. Figures 6

through 8 show the probability results from these cases. In Figures 6-8, the black ellipse is the 99%

lightning location uncertainty ellipse. The white circle is a 0.45 nrni radius around the point ofinterest.

As with the previous real-world tests, all calculated probabilities were consistent with these additional

real-world events.

Table 2. Input values used for scenarios shown in Figures 3 through 5.

Figure Semi- Semi- Confidence Heading Point of Point of Strike Strike Radius

major axis major axis (from interest interest latitude longitude around

of50% of50% true latitude longitude C'N) C'W) point

confidence confidence North) C'N) C'W) of

ellipse ellipse ofsemi- interest

(lan) (lan) major (nmi)

axIS

3 0.4 0.2 0.99 300.7 28.60827 -80.6041 28.6114 -80.6113 0.45

4 0.3 0.2 0.99 293 28.60827 -80.6041 28.6178 -80.6069 0.45

5 0.2 0.1 0.99 20.3 28.60827 -80.6041 28.5995 -80.6113 0.45

I Vaisala, Inc., 2006, "Vaisala StrikeNet Information," URL: http://www.vaisala.com/files/StrikeNet-Brochure.pdf

11
-----------------------------------------------------

Figure 3. Google Maps visualization ofthe 99%

confidence uncertainty ellipse for one ofthe closest

lightning strikes to Complex 39A on 03 August

2009.

Figure 4. Google Maps visualization ofthe 99%

confidence uncertainty ellipse for a lightning strike

near Complex 39A on 03 August 2009.

12
------------~--------------------------------

Figure 5. Google Maps visualization ofthe 99%

confidence uncertainty ellipse for nearby lightning

strike to Complex 39A on 03 August 2009.

Table 3. Input values used for scenarios shown in Figures 6 through 8.

Figure Semi- Semi- Confidence Heading Point of Point of Strike Strike Radius

major axis major axis (from interest interest latitude longitude around

of 50% of 50% true latitude longitude ("N) ("W) point

confidence confidence North) ("N) ("W) of

ellipse ellipse of semi- interest

(kIn) (kIn) major (nmi)

axis

6 0.6 0.4 0.99 82 28.60827 -80.6041 28.6069 -80.6087 0.45

7 0.4 0.4 0.99 95 28.60827 -80.6041 28.6057 -80.6085 0.45

8 0.2 0.1 0.99 72 28.62716 -80.6275 28.6275 -80.6202 0.45

13
14

· .. Figure 6. lllustrates a probability of 69.1% of a

lightning strike of amplitude -43.0 kA detected by

NLDN occurring 0.26 nmi from the center of

Launch Complex 39A on 8/16/2009.

Figure 7. lllustrates a probability of 74.7% of a

lightning strike of amplitude -71.4 kA detected by

NLDN occurring 0.28 nautical miles from the

center of Launch Complex 39A on 10/14/2009.
Figure 8. illustrates a probability of 99.9996% of a

lightning strike of amplitude -21.7 kA detected by

CGLSS occurring 0.04 nmi from the center of

Launch Complex 39B on 6/2712009.

Other Applications

The techniques and methods described in this paper clearly have application reaching far beyond

the space program uses for which they were designed. The list of potential applications is many and

varied and would be of interest to anyone seeking information pertaining to probability oflightning strike

locations, such as the power industry, aviation, or any industry sensitive to electrical overloads. This

methodology is also applicable to forensic meteorology [13], where the question of whether lightning

struck at or near a particular location is an issue in litigation. This technique could be applied to NLDN

or any other cloud-to-ground lightning detection system that provides location error ellipses.

Conclusion

A probabilistic facility-centric approach to lightning location has been developed to calculate the

probability that any nearby cloud-to-ground lightning stroke occurred within any radius of any point of

interest. In practice, this provides the probability that a nearby lightning stroke was within a key distance

of a facility, rather than within the error ellipses centered on the stroke. This process uses the bivariate

Gaussian distribution of probability density provided by the current lightning location error ellipse for the

most likely location of a lightning stroke and integrates it to determine the probability that the stroke is

inside any specified radius. This new facility-centric technique was tested extensively and is much more

useful to the space launch customers.

15
Appendix A - Example Lightning Probability Calculation

This appendix is an example of calculating the probability of any lightning stroke with a known

error ellipse being within a circle of any radius around any point. It is provided to clarify the calculation

process. An example calculation is shown in Table 4.

This example is a real-world event from a lightning strike near-the Space Shuttle launch pad 39A

at 02:35 GMT on 16 August 2009 (ref. Figure 6). Although the lightning data usually are from the cloudto-ground component ofthe Four Dimensional Lightning Surveillance System (CG-4DLSS) (Murphy et

a12008, Roeder 2010) in this example a lightning stroke from the NLDN is used. We sometimes use

StrikeNet reports that provide stroke data from the NLDN to double check the CG-4DLSS report. The

lightning stroke input values are shown in row 1 ofTable 3.

The following assumptions are applicable to the example calculation in Table 4. The location of

Launch Pad 39A is 28.60827486 north latitude (or 0.499309 radians) and 80.60411653 west longitude (or

-1.406807 radians). This is also the center ofthe circle in which the lightning probability will be

calculated. The desired radius for probability oflightning around Launch Pad 39A is 0.45 nautical miles

(nmi). This lightning stroke occurred on 16 August 2009 at 0235 GMT.

Table 4. Lightning strike probability calculation process.

Step Action

Convert semimajor and

semi-minor

axes fromkm

tonmi

Equation and other information

1 nmi = 1.852 km

Calculation and Result

0.6 km = 0.324 nmi semimajor axis

0.4 km = 0.216 nmi semi- . . mmoraXlS

2

3

16

Calculate

distance from

lightning stroke

to center of

circle

Perform

coordinate

system rotation

Haversine Distance Formula:

Distance = Earth Radius * C

• Earth Radius = 3443.920086 nmi

• C = 2 *Atn2[.J1- A, JA)

Atn2 is a two parameter arc tangent

function which returns values in all four

quadrants.

• A = sin(dlat/2)*sin(dlat/2) + cos(

(target lat))* cos«stroke lat))

*sin(dlon/2)*sin(dlon/2)

• dlat = latitude difference from

target to stroke =

28.60827° - 28.6069° =

2.39959x10·5 (radians)

• dlon = longitude difference from

circle to stroke = -80.60411°-

80.6085° = 7.99967X 10.5

(radians)

• X = (Longitude of Target - Longitude

of Stroke) * Cos (Latitude of Strike)

Distance = 3443.920086 *

7.42l7xlO·5

, = 0.2556 nmi

C = 2*Atn2{sqr(1 -

1.377x10·9

),

sqr(1.377x10·9)}

= 7.4217x10·5

A = sin(1.200x10·5/2)*

sin(1.200x1 0.5

/2)*

+

cos(0.4993)*cos(0.4993)*

sin(4.000x1 0.5/2)

*sin(4.000x10·5

/2)

= 1.3770x10·9

X = (-1.4068- (-1.4069) *

Cos (0.4993)
to eliminate the • Z = Latitude of Target - Latitude of = 7.0231 X 10-5

off-diagonal Stroke Z = 0.49931 - 0.49929

term in the • 8 = a - ((nI2) - Atn2(X,Z) = 2.400 X 10-5

covariance • a is the orientation angle ofthe 50%

matrix. lightning positional confidence ellipse ·8 = 1.431-1.571 -

Atn2(7.023 X 10-5

,2.400 • X' = miss distance*Cos(8 (coordinate X 10-5

) = 0.1896

system rotation angle))

• Z' = miss distance*Sin(8) X' = 0.2556*Cos(0.1896)

= 0.2510

Z' = 0.2556*Sin(0.1896)

= 0.0482

4 Calculate the • O'x' = Semi-major axis ofthe 50% k=~-2*ln(1-0.50) =

standard lightning positional confidence ellipse 1.1774

deviations in lelliptical scaling constant used to scale O'x' =0.3240/1.1774 =

the new rotated standard error to the 50% confidence 0.2752

coordinate level O'z' = 0.216011.1774 = system. • O'z' = Semi-minor axis ofthe 50% 0.1834

lightning positional confidence ellipse

lelliptical scaling constant used to scale

standard error to the 50% confidence

level

• Elliptical scaling constant, k, is:

~- 2 *1n(1- probability)

5 Calculate the

probability that

lightning stroke

was within the

target area of

interest by

performing a

numerical

integration

using

Simpson's rule

ofthe lightning

uncertainty

ellipse over the

area ofthe

circle around

the target of

interest.

r;;; -(H - flH)2 12CT1 -(H + flH )2/2CT1

I "W e +e p = ----,,~- J

2,[21iCTH 0

erf(ZI )+ erf(Z2)~H

ZI = [~(W _H2)- Jlk J; v'tO"K

Z2 = [ ~(W _H2)+ JlK J;v'2O"K

• flK and flH are the coordinates ofthe target

circle in the (X', Z') coordinate system

• O'K and O'H are equal to O'x' and O'z' which are

the standard deviations ofthe diagonalized

covariance matrix.

• W = Radius around target2

• N = the number ofiterations to perform

in the integration (for this example, N is

set to 200).

• DH = .JW I N = integration step

• B, C, and D are intermediate variables in

W = Radius around targee

W = 0.452 = 0.2025

DH=.JW IN

DH = .JO.2025 1200 =

0.00225

B = v'2 * O'x'

B = v'2 * 0.2752 = 0.3891

C=X'/B

C = 0.251010.3891 =

0.6451

D = 1 1(2 * ..J21i * O'z')

D = 1/(2 *..J21i * 0.1834) =

1.0874

H = iteration no. * DH

H = 199*0.00225

H = 0.4478

17
18

the algorithm corresponding to various

parts ofthe probability equation shown

above

• B = J2 * ax,

• C=X'/B

• D = 1/(2 * fbi * ad

• A, H, zl, z2, E, F, Erf(zl), Erf(z2), Q and

sum are intennediate variables in the

algorithm corresponding to various parts

ofthe probability equation shown above

• A loop is perfonned for j = 1 to 199.

This example is shown for j = 199.

• Sum=O

• Begin Loop here: H = j* DH

.A=~W-H2

• zl = AIB-C

• z2=AJB+C

• Erf(x) = error function =

fe-t2 dt

",:r 0

• E = (H - Z') 2/(2 * az'2)

• F = (H + Z') 2 1(2 * az'2)

• Q(j) = (e-E+ eF

) * (Erf(Zl) +

Erf(Z2))

• sum = sum + (3 - (-1) j) * Q(j)

• sum = sum + Q(O) + Q(N)

• Probability = D*sum*DH/3

A=~W-H2

A = ~r-O.-'20'-"-2-5--0-.4-4-77-5-2

A = 4.494 X 10-2

zl =A/B-C

zl = 4.494 X 10-2/0.3891

- 0.6451

zl = -0.5296

z2=A/B+C

z2 = 4.494 X 10,2/0.3891

+ 0.6451

z2 = 0.7606

Erf(ZI) =

ErrorFunction(z1)

Erf(-0.5296) = -0.5461

Erf(Z2)

ErrorFunction(z2)

Erf(0.7606) = 0.7179

E = (H - Z') 21 (2 * az'2)

E = (0.4478-

0.0482//(2*0.18342)

E = 2.372

F = (H + Z'/1(2 * az'2)

F = (0.4478 +

0.0482/1(2*0.18342)

F = 3.654

Q(j) = (e,E + e,F) *

(Erf(ZI) + Erf(Z2))

Q(199) = (e,2.372 + e'3.654)*

(-0.5461+ 0.7179) =

2.0467 X 10,2

sum = sum + (3 - (-1) j) *

Q(j)

sum = sum + (3 - (-1) 199)

* 2.047 X 10'2

sum = 844.8952

End Loop

sum = sum + Q(O) + Q(N)

sum = 844.8952 + 2.9361

+ 0 = 847.8317

Probability = D * sum *

DH/3

Probability = 1.0874 *

847.8317 * 0.00225 13=
1

0.6914 E ="(H

- Z') 21 (2*()Z'2)

E= (0.4478-

0.0482il(2*0.1834

2)

E= 2.372

F

= (H+ Z,)21(2*()Z'2)

F= (0.4478+

0.0482)2/(2*0.1834

2)

F= 3.654

Q(j)

= (eE+ eF)* (Erf(Zl) + Erf(Z2))

Q(199)

= (e2

.

372

+ e3

.

654)*

(-0.5461+ 0.7179)=

2.0467X 10-

2

sum= sum+ (3- (-1) j)* Q(j)

sum= sum+ (3- (-1) 199)

* 2.047X 10-

2

sum= 844.8952

End Loop

sum= sum+ Q(O)+ Q(N)

sum= 844.8952+ 2.9361

+0= 847.8317

Probability

=D*sum* DH/3

Probability= 1.0874* 847.8317 * 0.00225 13=

0.6914I

,I

19
Appendix B - Excel VBA Macro for 45th Weather Squadron Lightning

Spreadsheet that Calculates the Probability of any Nearby Lightning Stroke

being Inside any Radius of Any Point ofInterest

Conversions Module

Const pie As Double = 3.14159265358979

Const ecc As Double = 0.081819190842622 'eccentricity (eA )

Const kmPerMile As Double = 1.852

-------------------------------------------,-----------------

Function toRad(Degrees As Double) As Double

'Converts degrees to radians

toRad = Degrees * pie / 180

End Function

Function toDeg(Radians As Double) As Double

'Converts radians to degrees

toDeg = Radians * 180 / pie

End Function

Function applyConfidence(length As Double, newSigma As Double, currentSigma As Double) As

Double

'convert length based on confidence interval

applyConfidence = length * newSigma / currentSigma

End Function

Function revertConfidence(length As Double, newSigma As Double, currentSigma As Double) As

Double

'convert length based on confidence interval

revertConfidence = length * currentSigma / newSigma

End Function

Function kmToNmi(kilometers As Double) As Double

'Converts kilometers to nautical miles

kmToNmi = kilometers / kmPerMile

End Function '

Function geodToECEF(Lat As Double, Lon As Double) As Variant

'convert lat-Ion to EFG

Dim ret(2) As Double

ret(O) = getE(Lat, Lon)

ret(1) = getF(Lat, Lon)

ret(2) = getG(Lat)

geodToECEF = ret

End Function

20
\

Function EceIToGeod(E As Double, F AS,Double, g As Double) As Variant

'convert EFG coordinate to lat-lon

Dim ret(1) As Double 'array of(lat,lon)

ret(O) = Atn(g / ((1 - ecc /\ 2) * Sqr(E /\ 2 + F /\ 2») * 180/ pie

ret(1) = Atn2(E, F) * 180 / pie

EceIToGeod = ret

End Function

Functions Module

Const pie As Double = 3.14159265358979

Const eRad As Double = 3443.920086 'radius of earth in nautical miles

Const ecc As Double = 0.081819190842622 'eccentricity (e/\)

Public Function Atn2(ByVal X As Double, ByVal Y As Double) As Double

Dim at As Double

IfX = 0 Then

IfY < 0 Then

Atn2 = -pie / 2

ElseIf Y > 0 Then

Atn2 = pie/ 2

Else

Atn2 = 0

End If

ElseIf Y = 0 Then

If X > 0 Then

Atn2 = 0

ElseIf X < 0 Then

Atn2 = pie

End If

Else

at = Atn(Y / X)

If X < 0 And Y < 0 Then

Atn2 = at - pie

ElseIf X < 0 And Y > 0 Then

Atn2 = at + pie

Else

Atn2 = at

End If

End If

End Function

Function arcCos(X As Double)

arcCos = Atn((l - X /\ 2) /\ 0.5 / X),

End Function

21
Function addPoints(p10 As Double, p20 As Double) As Variant

'add two EFG coordinates together (

Dim ret(2) As Double

ret(O) = p1(0) + p2(0)

ret(1) = p1(1) + p2(1)

ret(2) = p1(2) + p2(2)

addPoints = ret

End Function

Function getEPointDiffs(eLat As Double, eLon As Double, A As Double, B As Double, phi As Double, t

As Double) As Variant

'get EFG differentials for a point on ellipse at t radians

Dim ret(2) As Double '

Dim xl As Double, y1 As Double

xl = getX(A, B, 0, 90 - phi, t)

y1 = getY(A, B, 0,_ 90 - phi, t)

ret(O) = -xl * Sin(eLon * pie / 180) - yl * Sin(eLat * pie / 180) * Cos(eLon * pie / 180)

ret(1) = X:1 * Cos(eLon * pie / 180) - y1 * Sin(eLon * pie / 180) * Sin(eLat * pie / 180)

ret(2) = y1 * Cos(eLat * pie / 180)

getEPointDiffs = ret

End Function

Function getFocalPointDiffs(eLat As Double, eLon As Double, A As Double, B As Double, phi As

Double) As Variant"

'get EFG differentials for focal point of ellipse

Dim ret(2) As Double

Dim fodDist As Double, xl As Double, yl As Double

fodDist = Sqr(A A 2 - B A 2)

. xl = fodDist * Cos(toRad(90 - phi»

y1 = fodDist * Sin(toRad(90 - phi»

ret(O) = -x1 * Sin(eLon * pie / 180) - Y1 * Sin(eLat * pie / 180) * Cos(eLon * pie / 180)

ret(l) = xl * Cos(eLon * pie /180) - y1 * Sin(eLon * pie /180) *iSin(eLat * pie /180)

ret(2) = y1 * Cos(eLat * pie / 180) .

getFocalPointDiffs = ret

End Function

Function getDistGeod(latl As Double, lon1 As Double, lat2 As Double, lon2 As Double) As Double

'spherical

Dim rLatl As Double, rLat2 As Double, rLonl As Double, rLon2 As Double

rLatl = toRad(latl)

rLat2 = toRad(lat2)

rLonl =toRad(lon1)

rLon2 = toRad(lon2)

22
getDistGeod = getN«(latl + lat2) / 2) * arcCos(Cos(rLatI) * Cos(rLonI) * Cos(rLat2) * Cos(rLon2) +

Cos(rLatl) * Sin(rLonI) * Cos(rLat2) * Sin(rLon2) + Sin(rLatl) * Sin(rLat2)) / 180 * pie

End Function

Function greatCircle(latl As Double, lonI As Double, lat2 As Double, lon2 As Double) As Double

'great circle distance in NM

Dim rLatl As Double, rLat2 As Double, rLonI As Double, rLon2 As Double

rLatl = toRad(latl)

rLat2 = toRad(lat2)

rLonl = toRad(lonl)

rLon2 = toRad(lon2)

greatCircle = (eRad) * arcCos(Cos(rLatI) * Cos(rLat2) * Cos(rLon2 - rLonI) + Sin(rLatl) * Sin(rLat2))

End Function

Function Haversine(latl As Double, lonI As Double, lat2 As Double, lon2 As Double) As Double

'Calculates distance between two latitude and longitude points by the haversine formula in NM

Dim dlat As Double, dlon As Double, C As Double, A As Double

dlat = toRad(latl) - toRad(lat2)

dlon = toRad(lonI) - toRad(lon2)

A = Sin(dlat / 2) * Sin(dlat / 2) + Cos(toRad(latl)) * Cos(toRad(lat2)) * Sin(dlon / 2) * Sin(dlon / 2)

C = 2 * Atn2(Sqr(1 - A), Sqr(A))

Haversine = eRad * C

End Function

Function getX(A As Double, B As Double, H As Double, phi As Double, t As Double) As Double

'get the X cartesian coordinate of a point on the ellipse at t radians

getX = H + A * Cos(t) * Cos(phi * pie / 180) - B * Sin(t) * Sin(phi * pie / 180)

End Function

Function getY(A As Double,B As Double, k As Double, phi As Double, t As Double) As Double

'get the Y cartesian coordinate of a point on the ellipse at t radians

getY = k + B * Sin(t) * Cos(phi * pie / 180) + A * Cos(t) * Sin(phi * pie / 180)

End Function

Function getN(Lat As Double) As Double

'get radius ofthe earth as a given latitude, adjusted for eccentricity

getN = eRad / (1 - «ecc /\ 2) * Sin(Lat * pie /180) /\ 2)) /\ 0.5

End Function

.Function getE(Lat As Double, Lon As Double) As Double

23
/

'get E coordinate from lat-lon

getE = getN(Lat) * Cos(Lat * pie / 180) * Cos(Lon * pie / 180)

End Function

Function getF(Lat As Double, Lon As Double) As Double

'get F coordinate from lat-lon

getF = getN(Lat) * Cos(Lat >Itpie / 180) * Sin(Lon * pie / 180)

End Function

Function getG(Lat As Double) As Double

'get G coordinate from lat

getG = getN(Lat) *(1 - ecc /\ 2) * Sin(Lat * pie /180)

End Function

Function getEPoint(ecLat As Double, ecLon As Double, A As Double, B As Double, phi As Double, t As

Double) As Variant,

'get lat-lon of a point on ellipse at t radians

Dim ret(l) As Double

Dim eCO As Double

Dim pointO As Double

Dim diffsO As Double

eC = geodToECEF(ecLat, ecLon)

diffs = getEPointDiffs(ecLat, ecLon, A, B, phi, t)

point = addPoints(eC, diffs)

getEPoint = EcefToGeod(point(O), point(I), point(2))

End Function

, Function ePointLat(ecLat As Double, ecLon As Double, A As Double, B As Double, phi AsDouble, t As

Double) As Double

'get lat of a point on ellipse at t radians

Dim geodO As Double

geod = getEPoint(ecLat, ecLon, A, B, phi, t)

ePointLat = geod(O)

End Function

Function ePointLon(ecLat As Double, ecLon As Double, A As Double, B As Double, phi As Double, t

As Double) As Double

'get Ion of a point on ellipse at t radians

Dim geodO As Double

geod = getEPoint(ecLat, ecLon, A, B, phi, t)

ePointLon = geod(l)

End Function

/

Function getAzimuth(latl As Double, loni As Double, lat2 As Double, 10n2 As Double) As Double

'get azimuth between two points

Dim angle As Double

24
angle = (Atn2(Cos(toRad(latl)) * Sin(toRad(lat2)) - Sin(toRad(latl)) * Cos(toRad(lat2)) * Cos(toRad(lon2) - toRad(lonl)), Sin(toRad(lon2) - toRad(lonl)) * Cos(toRad(lat2))) * 180/ pie)

If (angle = 0) Then

getAzimuth = 360

Else

getAzimuth = (angle + 360) Mod 360

End If

End Function

Module 2 Code

Dim point As New Strike

Const passWord As String = "PAEcOd3m0nk3y$"

Sub POI_ChangeO

Calculate

End Sub

Sub calculateRowO

Application.EnableCancelKey =xlDisabled

Worksheets("Results").protect Contents:=False, passWord:=passWord

Dim ecefO As Double

Dim interestPoint(1) As Double' (lat,lon)

Dim rowNumber data As Double

Dim rowNumber result As Double

Dim rowData As Variant

Dim AOI(3) As Double

Dim lastRow As Double

IfWorksheets("Results").FilterMode = True Then

Worksheets("Results").ShowAllData

End If

, lastRow = Worksheets("Results").Cells(Rows.Count, "A").End(xlUp)~row

If (lastRow > 6) Then

Worksheets("Results").Range("A7:R" & lastRow).Select

Selection.Delete

End If

AOI(O) = Range("latLower")

AOI(1) = Range("latUpper")

AOI(2) = Range("lonLower")

AOI(3) =Range("lonUpper")

interestPoint(O) = Range("poiLat")

interestPoint(1) = Range("poiLon")

ecef = geodToECEF(interestPoint(O), interestPoint(1))

25
Call point.setStatics(interestPoint, ecef, Range("newSigma"), Range("currentSigma"),

Range("radius"))

Dim xlCalc As XICalculation

xlCalc = Application.Calculation

Application.Calculation = xlCalculationManual

, On Error GoTo CalcBack

lastRow = Worksheets("Data").Cells(Rows.Count, "A").End(xlUp).row

rowNumber result = 7

For rowNumber data = 2 To lastRow

rowData = Worksheets("Data").Range("A" & rowNumber_data & ":1" & rowNumber_data)

If (rowData(l, 3) > AOI(O) And rowData(l, 3) < AOI(I) And rowData(1, 4) > AOI(2) And

rowData(1, 4) < AOI(3)) Then

Call point.init(rowData(l, 3), rowData(1, 4), rowData(l, 6), rowData(l, 7), rowData(l, 8), True)

Worksheets("Results").Cells(rowNumberJesult, 1) = rowData(l, 1)

Worksheets("Results").Cells(rowNumberJesult, 2) = rowData(l, 2)

Worksheets("Results").Cells(rowNumberJesult, 3) = point.centerAzimuth

Worksheets("Results").Cells(rowNumberJesult, 4) = point.centerRange

Worksheets("Results").Cells(rowNumberJesult, 5) = point.criticalAzimuth

Worksheets("Results").Cells(rowNumberJesult, 6) = point.criticalRange

Worksheets("Results").Cells(rowNumberJesult, 7) = rowData(l, 5)

Worksheets("Results").Cells(rowNumberJesult, 8) = point.isInside

Worksheets("Results").Cells(rowNumberJesult, 9) = rowData(l, 9)

Worksheets("Results").Cells(rowNumberJesult, 10) = point.Prob_simpson

Worksheets("Results").Cells(rowNumberJesult, 11) = point.rangeDifference

Worksheets("Results").Cells(rowNumberJesult, 12) = point.Lat

Worksheets("Results").Cells(rowNumberJesult, 13) = point.Lon

Worksheets("Results").Cells(rowNumberJesult, 14) = point.A

Worksheets("Results").Cells(rowNumberJesult, 15) = point.B

Worksheets("Results").Cells(rowNumberJesult, 16) = point.phi

Worksheets("Results").Cells(rowNumberJesult, 17) = point.criticalLat

Worksheets("Results").Cells(rowNumberJesult, 18) = point.criticalLon

'Worksheets("Results").Cells(rowNumberJesult, 19) = point.ProbJician

'Worksheets("Results").Cells(rowNumberJesult, 20) = point.Prob_simpson

rowNumber result =rowNumber result + 1 - - End If

Next rowNumber data

lastRow = Worksheets("Results").Cells(Rows.Count, "A").End(xlUp).row

Worksheets("Results").Range("A7:A" & lastRow).NumberFormat = "MMIDDNYYY"

Worksheets("Results").Range("B7:B" & lastRow).NumberFormat = "HH:mm:ss.OOO"

Worksheets("Results").Range("C7:C" & lastRow).NumberFormat = "000"

Worksheets("Results").Range("D7:D" & lastRow).NumberFormat = "0.00"

Worksheets("Results").Range("E7:E" & lastRow).NumberFormat = "000"

Worksheets("Results").Range("F7:F" & lastRow).NumberFormat = "0.00"

Worksheets("Results").Range("J7:J" & lastRow).NumberFormat = "0.000000%"

Worksheets("Results").Range("K7:K" & lastRow).NumberFormat = "0.00"

Worksheets("Results").Range("L7:M" & lastRow).NumberFormat = "0.0000"

Worksheets("Results").Range("N7:P" & lastRow).NumberFormat = "0.00"

Worksheets("Results").Range("Q7:R" & lastRow).NumberFormat = "0.0000"

26
'Worksheets("Results").Range("S7:T" & lastRow).NumberFormat = "0.000000%"

If (lastRow > 6) Then

With Worksheets("Results").Range("H7:H" & lastRow).Select

Selection.FormatConditions.Add Type:=xlTextString, String:="Yes", _

TextOperator:=xlContains

Selection.FormatConditions(Selection.FormatConditions.Count).SetFirstPriority

With Selection.FormatConditions(l).Font

.Color = -16777024

.TintAndShade = 0

End With

With Selection.FormatConditions(1).Interior

.PattemColorIndex = xlAutomatic

.Color = 5263615

.TintAndShade = 0

End With

End With

End If

Application.Calculation = xlCalc

Worksheets("Results").EnableAutoFilter = True

Worksheets("Results").protect Contents:=True, passWord:=passWord, AllowFiltering:=True,

AllowSorting:=True, userInterfaceOnly:=True

MsgBox ("Data Calculation Complete")

Exit Sub

CalcBack:

Worksheets("Results").protect Contents:=True, passWord:=passWord, AllowFiltering:=True,

AllowSorting:=True, userInterfaceOnly:=True

Application.Calculation = xlCalc

End Sub

Sub openBrowserO

Dim rowNumber As DoubleDim N As Integer

Dim ecefO As Double

Dim interestPoint(l) As Double' (lat,lon)

Dim showPath As Boolean

showPath = ActiveSheet.showPath.value

interestPoint(O) = Range("poiLat")

interestPoint(1) = Range("poiLon")

ecef = geodToECEF(interestPoint(O), interestPoint(1))

27
Call point.setStatics(interestPoint, ecef, Range("newSigma"), Range("currentSigma"),

Range("radius"))

rowNumber = Range("nRow")

N = Range("nPoints")

Call pointinit(Cells(rowNumber, 12), Cells(rowNumber, 13), Cells(rowNumber, 14),

Cells(rowNumber, 15), Cells(rowNumber, 16), False)

Call point.setURL(N, showPath)

'Range("urlString") = point.URL

Set browser = CreateObject("IntemetExplorer.App1ication")

If (Len(point.URL) > 1950) Then

MsgBox ("Your request contains too many characters. Please select a smaller number of perimeter

points or tum on offthe path attribute. ")

Else

browser.Navigate (point.URL)

browser.Visible = True.

End If

End Sub

Probability Module

Const pie As Double = 3.14159265358979

Function RicianIntegral(O As Double, V As Double, m As Integer)

'This subroutine computes the Rician Integral Pm when u, v and m are given.

'It is preferred because both the exponentials are outside the summation loops.

'This formulation gives essentially the same numerical value as the other one.

Dim tm As Double, rm As Double, sm As Double, sumTm As Double, sumTmSm As Double

On Error GoTo overflow

tm= 1

rm= 1

sm= 1

sumTm= 1

sumTmSm= 1

Fori = 1 To m

tm = tm * V / (2 *i)

rm = rm * U / (2 *i)

sm= sm +rm

sumTm = sumTm + tm

sumTmSm = sumTmSm + tm *sm

Next i l

RicianIntegral = Exp(-V / 2) *sumTm - Exp(-(0 + V) / 2) *sumTmSm

Exit Function

overflow:

28
1 _

RicianIntegral = 0

End Function

Function numericIntegral(Xprime As Double, Zprime As Double, SigmaXp As Double, SigmaZp As

Double, Xe As Double, Ra As Double, N As Integer)

Dim xl As Double, zl As Double, r1 As Double, x2 As Double, z2 As Double, r2 As Double

Dim R As Double, sum As Double, dTheta As Double, thisPhi As Double, Q As Double

DPhi = 2 * pie 1N

sum=O

thisPhi = 0

xl = (Xprime + Ra) * SigmaZp 1SigmaXp

zl = Zprime

For i = 1 To N

thisPhi = thisPhi + DPhi

x2 = (Xprime + Ra * Cos(thisPhi)) * SigmaZp 1SigmaXp

z2 = Zprime + Ra * Sin(thisPhi)

r1 = Sqr(xl * xl + zl * zl)

r2 = Sqr(x2 * x2 + z2 * z2)

dTheta = Application.WorksheetFunction.Asin((xl * z2 - x2 *zl) 1(rl * r2))

R = (r1 + r2) 12

sum = sum + Exp(-(R 1SigmaZp) /\ 2 12) * dTheta

xl.= x2

zl = z2

Next i

Q = -sum 1(2 * pie)

If Xe < Ra Then

numericIntegral = 1 + Q

ElseIf Xe = Ra Then

numericIntegral = 1 12 + Q

Else

numericIntegral = Q

End If

End Function

Function simpsonIntegral(Xprime As Double, Zprime As Double, SigmaXprime As Double,

SigmaZprime As Double, Ra As Double, N As Integer)

'This program computes the collision probability between two objects.

'The integration is performed using Simpson's Rule.

Dim Q(1000) As Double, P As Double

Dim A As Double, B As Double, C As Double, D As Double, E As Double, F As Double, W As

Double

Dim R As Double, DR As Double

29
pi = 3.14159265358979

B = Sqr(2) * SigmaXprime

C = Xprime I B

D = 1 1(2 * Sqr(2 * pi) * SigmaZprime)

W=Ra I\2

DH = Sqr(W) IN

'Compute Q(I) array from 0 to (N-l)

For j = 0 To (N - 1)

H=j *DH

A = Sqr(W - H 1\ 2)

zl = AlB - C '

z2=A/B+C

ErfZl = ErrorFunction(zl)

ErfZ2 = ErrorFunction(z2)

E = (H - Zprime) 1\ 2 I (2 * SigmaZprime 1\ 2)

F = (H + Zprime) 1\ 2 I (2 * SigmaZprime 1\ 2)

Q(j) = (Exp(-E) + Exp(-F)) * (ErfZl + ErfZ2)

Nextj

Q(N) = 0

'Compute the integral by Simpson's Rule

sum=O

For j = 1 To (N - 1)

sum = sum + (3 (-1) 1\ j) * Q(j)

Nextj

sum = sum + Q(O) + Q(N)

P = D * sum * DH I 3

simpsonIntegral = P

End Function

Function ErrorFunction(X)

'This subroutine computes the Error Function when X is given.

Dim Z, SqrtPi, tm, sum, Max, m, ErtX, ErfcX

SqrtPi = 1.77245385090552

30
.-------------------------------------------------

Z=Abs(X)

IfZ <= 4 Then

tm=Z

sum=tm

Max = Int(20 * Z) + 1

For m = 1 To Max

tm = -tm * Z * Z *(2 * m - 1) / (m *(2 * m + 1))

sum = sum +tm

Nextm

ErfX = sum * 2 / SqrtPi

E1seIfZ <= 5.736 Then

tm= 1

sum = 1

Form= 1 To 1

tm = -tm *(2 * m - 1) / (2 * Z A 2)

sum = sum + tm

Nextm

ErfcX = sum * Exp(-Z A 2) / (SqrtPi *Z)

ErfX = 1 - ErfcX

Else

ErfX= 1

End If

ErrorFunction = Sgn(X) * ErfX

End Function

Strike Class Module

Const pie As Double = 3.14159265358979

Const eRad As Double = 3443.920086 'radius of earth in nautical miles

Const ecc As Double = 0.081819190842622 'eccentricity (eA )

Private p--poi_geoO As Double

Private p--poi_ecefO As Double

Private p_newSigma As Double

Private p_currentSigma As Double

Private p_strike_geo(1) As Double

Private p_major As Double

Private p_minor As Double

Private p--phi As Double

'lat Ion point ofinterest [lat,lon]

'EFG point ofinterest [E,F,G]

'user defined confidence level

'50% confidence

'lat Ion ofstrike [lat,lon]

'major axis in NM

'minor axis in NM

'stroke orientation in degrees

Private p_a As Double 'scaled major axis in NM

Private p_bAs Double 'scaled minor axis in NM

Private p_strike_ecefO As Double 'EFG strike [E,F,G]

Private p_criticalT As Double 'angle of closest point in radians

Private p_critical_geoO As Double 'lat Ion of closest point [lat,lon]

Private p_critical_ecef(2) As Double 'EFG of closest point [E,F,G]

Private p_aziCenter As Double 'azimuth to strike center

Private pJangeCenter As Double 'distance to strike center

31
Private p_aziCritical As Double 'azimuth to closest point

Private pJangeCritical As Double 'distance to closest point

Private pJangeDifference As Double

Private p_isInside As String

Private p_URL As String

Private pJoci10 As Double

Private pJoci20 As Double

Private pJJrob_rician As Double

Private pJJrob_numeric As Double

Private pJJrob_simpson As Double

Private pJadius As Double 'radius around point ofinterest

Public Sub setStatics(poi_geo As Variant, poi_ecefAs Variant, newSigma As Double, currentSigma As

Double, radius As Double)

pJJoi_geo = poi_geo

pJJoi_ecef= poi_ecef

p_newSigma = Sqr(-2 * Log(1 - newSigma /100))

p_currentSigma = Sqr(-2 * Log(1 - currentSigma /100))

pJadius = radius

End Sub

Public Sub init(Lat, Lon, major, minor, phi, adjustAxis As Boolean)

On Error GoTo initErr

p_strike_geo(O) = Lat

p_strike_geo(1) = Lon

pJJhi = phi

If (adjustAxis) Then .

'input major and minor are in KM, and unscaled

p_major = kmToNmi(CDbl(major))

p_minor = kmToNmi(CDbl(minor))

,

\

Ifp_major = 0 Then

p_major = 0.05

End If

If p_minor = 0 Then

p_minor = 0.05

End If

p_a = applyConfidence(p_major, p_newSigma, p_currentSigma)

p_b = applyConfidence(p_minor, p_newSigma, p_currentSigma)

Else

'input major and minor are in NM and scaled

p_a = CDbl(major)

32
p_b = CDbl(minor)

Ifp_a = 0 Then

p_a = 0.1

End If

Ifp_b = 0 Then

p_b = 0.1

End If

p_major = revertConfidence(p_a, p_newSigma, p_currentSigma)

p_minor = revertConfidence(p_b, p_newSigma, p_currentSigma)

End If

p_strike_ecef = geodToECEF(p_strike_geo(O), p_strike_geo(l))

Call setT

p_critical_geo = getEPoint(p_criticalT)

p_aziCenter = getAzimuth(p-lJoi_geo(O), p-lJoi_geo(l), p_strike_geo(O), p_strike_geo(l))

p_aziCritical = getAzimuth(p-lJoi_geo(O), p-lJoi_geo(l), p_critical_geo(O), p_critical_geo(l))

pJangeCenter = Haversine(P-lJoi_geo(O), P-lJoi_geo(l), p_strike_geo(O), p_strike_geo(l))

pJangeCritifal = Haversine(p-lJoi_geo(O), p-lJoi_geo(l), p_critical_geo(O), p_critical_geo(l))

pJangeDifference = pJangeCenter - pJangeCritical .

Dim diffsO As Double

Dim pointO As Double

diffs = getFocalPointDiffs(p_strike_geo(O), p_strike_geo(l), p_a, p_b, p-lJhi)

point = addPoints(p_strike_ecef, diffs)

p_focil = EcefToGeod(point(O), point(l), point(2))

diffs(O) = -diffs(O)

. diffs(l) = -diffs(l)

diffs(2) = -diffs(2)

point = addPoints(p_strike_ecef, diffs)

pJoci2 = EcefToGeod(point(O), point(l), point(2))

Call setProbability

Exit Sub

initErr:

asd= 0

End Sub

Public Property Get AO

A=p_a

End Property

Public Property Get BO

B=p_b

End Property

Public Property Get criticalTO

criticalT = p_criticalT

End Property J

33
Public Property Get criticalLatO

criticalLat = p_critical_geo(O)

End Property

Public Property Get criticalLonO

criticalLon = p_critical_geo(l)

End Property

Public Property Get criticalAzimuthO

criticalAzimuth = p_aziCritical

End Property

Public Property Get centerAzimuthO

centerAzimuth = p_aziCenter

End Property

Public Property Get criticalRangeO

criticalRange = pJangeCritical

End Property

Public Property Get centerRangeO

centerRange = pJangeCenter

End Property

Public Property Get rangeDifferenceO

rangeDifference = pJangeDifference

End Property

Public Property Get URLO

URL=p_URL

End PropertyPublic Property Get LatO

Lat = p_strike_geo(O)

. End Property

Public Property Get LonO

Lon = p_strike_geo(I)

End Property

34
---------------------------------- --- ------------ -- - ----------------~----------------

Public Property Get phiO

phi = p-phi

End Property

Public Property Get ProbJicianO

ProbJician = p-probJician

End Property

Public Property Get Prob_numericO

Prob_numeric =p-prob_numeric

End Property

Public Property Get Prob_simpsonO

Prob_simpson = p-prob_simpson

End Property

Public Property Get isInsideO

Dim dl As Double, d2 As Double, tOO As Double, d3 As Double, tIO As Double

dl = Haversine(pjocil(O), pjocil(l), p-poi_geo(O), p-poi_geo(1))

d2 = Haversine(pjoci2(O), p_foci2(1), p-poi_geo(O), p-poi_geo(1))

to = getEPoint(O)

tI = getEPoint(pie)

da = Haversine(tl(O), tI(1), to(O), to(I))

If ((dl + d2) < (da)) Then

isInside = "Yes"

Else

isInside = "No"

End If

End Property

Private Sub setTO

Dim thisT As Double, minDist As Double, thisDist As Double, negDist As Double, posDist As

Double, returnT As Double, step As Double

Dim pointLat As Double, pointLon As Double

Dim pointO As Double

returnT = -pie

point = getEPoint(retumT)

minDist = Haversine(p-poi_geo(O), p-poi_geo(l), point(O), point(1))

For thisT = -pie To pie Step pie / 4

point = getEPoint(thisT)

thisDist = Haversine(p-poi_geo(O), p-poi_geo(1), point(O), point(1))

If (thisDist < minDist) Then

minDist = thisDist

35
returnT = thisT

End If

Next thisT

step = pie / 4

While step> pie / (2 /\ 16)

point = getEPoint(returnT - step)

negDist = Haversine(pyoi_geo(O), pyoi_geo(1), point(O), point(1))

point = getEPoint(returnT + step)

posDist = Haversine(pyoi_geo(O), pyoi_geo(1), point(O), point(1))

If (negDist < minDist)Then

returnT = returnT - step

minDist = negDist

ElseIf posDist < minDist Then

returnT = returnT + step

minDist = posDist

End If

step = step / 2

Wend

p_criticalT = returnT

End Sub

Private Function getEPoint(t As Double) As Variant

Dim ret(1) As Double

Dim pointO As Double

Dim diffsO As Double

diffs = getEPointDiffs(p_strike_geo(O), p_strike_geo(1), p_a, p_b, pyhi, t)

point = addPoints(p_strike_ecef, diffs)

getEPoint = EceffoGeod(point(O), point(I), point(2))

End Function

Private Function getCPoint(t As Double) As Variant

Dim ret(1) As Double

Dim pointO As Double

Dim diffsO As Double

diffs = getEPointDiffs(p_strike_geo(O), p_strike_geo(I), PJadius, pJadius, pyhi, t)

point = addPoints(pyoi_ecef, diffs)

getCPoint = EceffoGeod(point(O), point(I), point(2))

End Function

Public Sub setURL(nPoints As Integer, showPath As Boolean)

Dim thisT As Double, tStep As Double

Dim pointO As Double

p_URL = ''http://maps.google.com/staticmap?size=640x640&maptype=satellite&markers="

tStep = 2 * pie / nPoints

For i = 1 To nPoints

thisT = -pie + i *tStep

point = getEPoint(thisT)

36
p_URL = p_URL & Round(point(O), 5) & "," & Round(point(1), 5) & "~,If & "tinyyellowl"

'p_URL = p_URL & point(O) & "," & point(1) & "," & "tinyyellowl"

Next i

tStep = 2 * pie / 10

For i = 1 To 10

thisT = -pie + i *tStep

point = getCPoint(thisT)

p_URL = p_URL & Round(point(O), 5) & "," & Round(point(1), 5) & "," & "tinywhitel"

'p_URL = p_URL & point(O) & "," & point(1) & "," & "tinywhitel"

Next i

p_URL = p_URL & p_critica1_geo(0) & "," & p_critical_geo(1) & "," & "smallpurplel"

'p_URL = p_URL & pjoci1(0) & "," & pjoci1(1) & "," & "smallbluel"

'p_URL = p_URL & pjoci2(0) & "," & pjoci2(1) & "," & "smallbluel"

p_URL = p_URL & p-poi_geo(O) & "," & p-poi_geo(1) & "," & "smallgreenl"

p_URL = p_URL & p_strike_geo(O) & "," & p_strike_geo(1) & "," & "smallred"

If (showPath) Then

p_URL = p_URL & "&path="

tStep = 2 * pie / nPoints

For i = 1 To nPoints + 1

thisT = -pie + i *tStep

point = getEPoint(thisT)

p_URL = p_URL & Round(point(O), 5) & "," & Round(point(1), 5) & "I"

Next i

End If

p_URL = p_URL &

"&sensor=false&key=ABQIAAAAdwjRqRR8FuOdpAOoimTJCBSOxDK05lwxOGB6dIDkgLOxwdqZC

hSForDLWhvadXUn6EEI6WZWYt853w"

End Sub

Public Sub setProbabilityO

On Error GoTo probErr

'This program computes the probability ofof a lightning strike occurring within a

'specified distance of an asset ofinterest.

'The cross-section is a circle.

Dim alpha As Double, LonP As Double, LatP As Double, LonS As Double, LatS As Double

Dim Xprime As Double, Zprime As Double, SigmaXp As Double, SigmaZp As Double, sigma As

Double

Dim X As Double, Z As Double, theta As Double

Dim U As Double, V As Double

Dim m As Integer, AspectRatio As Double

p-prob_rician = 0_

p-prob_numeric = 0

If(p_major = 0 Or p_minor = 0) Then

Exit Sub

End If

LatP = toRad(p-poi_geo(O)) 'Latitude of POI (radians)

37
LonP = toRad(P---Foi.:-geo(I» 'Longitude ofPOI (radians)

LatS = toRad(p_strike_geo(O» 'Latitude of Strike Spot (radians)

LonS =toRad(p_strike_geo(1» 'Longitude of Strike Spot (radians)

alpha = toRad(p---Fhi) 'Heading of Semi~MajorAxis of 50% Confidence Ellipse (radians)

AspectRatio = p_major / p_minor

x = (LonP - LonS) * Cos(LatS)

Z = LatP - LatS

theta = alpha - «pie / 2) - Atn2(X, Z»

Xprime = pJangeCenter * Cos(theta)

Zprime = pJangeCenter * Sin(theta)

SigmaXp =p_major /1.177410023

SigmaZp = p_minor /1.177410023

sigma = Sqr(SigmaXp * SigmaZp)

'Compute Rician Integral

'U = (pJadius / sigma) A 2

'V = Xprime A 2 / SigmaXp A 2 + Zprime A 2 / SigmaZp A 2

'm = Int(Sqr(U * V» + 1

'p---Frob_rician = RicianIntegral(U, V, m)

'skip numeric method if aspect ratio is low

'If (m < 25 And AspectRatio < 10) Then

'Exit Sub

'End If

'Use numerIcal method for computing lightning strike probability:

'p---Frob_numeric = numericIntegral(Xprime, Zprime, SigmaXp, SigmaZp, pJangeCenter, pJadius,

1000)

p---Frob_simpson = simpsonIntegral(Xprime, Zprime, SigmaXp, SigmaZp, pJadius, 200)

Exit Sub

probErr:

End Sub

38
,

-----------------------------

References

Alfano, S., "Satellite Collision Probability Enhancements," J. Guidance, Control, and Dynamics, Vol. 29,

No.3, May-June 2006, pp. 588-592.

Alfano, S., "Review of Conjunction Probability Methods for Short-Term Encounters," AAS/AIAA Space

Flight Mechanics Meeting, Paper AAS 07-148, Sedona, AZ, 2007, 44 pp.

Austin, E., "Forensic Meteorology and Climatology: A Rapidly Emerging Niche," Third Annual CCM

Forum, 9(/' Annual American Meteorological Society Meeting, Paper 1.2, Atlanta, GA, 2010.

Beyer, W. H. (ed.), CRC Handbook of Tables for Probability and Statistics, The Chemical Rubber Co.,

Cleveland,OH, 1968, pp. 151-157.

Chan, F. Kenneth, Spacecraft Collision Probability, The Aerospace Press, El Segundo, CA, 2008, Chaps.

3,4,5,6

Chan, F. Kenneth, "Miss Distance - Generalized Variance Non-Central Chi Distribution," AAS/AIAA

Space Flight Mechanics Meeting, Paper AAS 11-175, New Orleans, LA, 2011, 12 pp.

Flinn, F. C., Roeder, W. P., Buchanan, M.D., and McNamara, T. M., "Lightning Reporting at 45th

Weather Squadron: Recent Improvements," 21 st International Lightning Detection Conference,

Orlando, FL, 2010a, 18 pp.

Flinn, F. C., Roeder, W. P., Pinter, D. F., Holmquist, S. M., Buchanan, M. D., McNamara, T. M.,

McAleenan, M., Winters, K.A., Gemmer, P. S., Fitzpatrick, M. E., and Gonzalez, R. D., "Recent

Improvements in Lightning Reporting at 45th Weather Squadron," 141h Conference on Aviation,

Range, and Aerospace Meteorology, Paper 7.3, Atlanta, GA, 2010b, 14 pp.

Leleux, D., Spencer, R., Zimmerman, P., Propst, C., Heilman, W., Frisbee, J., and Wortham, M.,

"Probability-Based Space Shuttle Collision Avoidance," Space Ops 2002 Conference, Houston, TX,

2002,20 pp.

Murphy, M. J., Cummins, K.L., Demetriades; N. W. S., and Roeder, W. P., "Performance of the New

Four-Dimensional Lightning Surveillance System (4DLSS) at the Kennedy Space Center/Cape

Canaveral Air Force Station Complex," 131h Conference on Aviation, Range, and Aerospace

Meteorology, Paper 7.6, New Orleans, LA, 2008, 18 pp.

( Patera, Russell, P., "General Method for Calculating Satellite Collision Probability," J Guidance,

Control, and Dynamics, Vol. 24, No.4, Jul.-Aug. 2001, pp. 716-722.

Roeder, W. P., Weems, J. W., and Wahner, P. B., "Applications of the Cloud-to-Ground Lightning

Surveillance System Database," ]'1 Conference on the Meteorological Applications of Lightning

Data, San Diego, CA, 2005, 5 pp.

Roeder, W.P., "The Four Dimensional Lightning Surveillance System," 21 st International Lightning

Detection Conference, Orlando, FL, 2010, 15 pp.

39