Adventures Among the Alarmists

I've been spending far too much time reading 'climate change' websites and papers lately. Mostly I wanted to see if my skepticism about lack of model validation in climate science was well founded (it is, claims of impending doom are not supported by any sort of validated model predictions). I came across a good blog that I added to my reading list, written by a scientist in Japan working on climate modeling, and his frustrations with alarmism and appreciation of Bayes / Jaynes reflect my own (nice to see I'm not a complete crank), however his ensemble evaluation method which claims to 'do away with the truth-centered paradigm' is a bit troubling from a validation perspective (it seems like giving up on correcting any bias, might have a post on this after I read it again more carefully). So just for finding that site, the descent into the internet climate-swamp was worth it.

That's the extent of the good, here's some of the bad (or at least amusing and misinformed). This site is run by an NGO employee with a biology background. It has the sort of alarmist cheer-leading and appeal to authority that you would expect from a motivated, but uncritical social worker. One of the interesting things is the challenge page, where the challenge is to debunk 'anthropocentric global warming'. Now, I've never heard that particular turn of phrase (I've always heard it as 'anthropogenic'), and I thought the challenge was a little vague, so I asked them to clarify. The response I got was surprising, apparently 'anthropocentric' involves a philisophical / moral connotation as well as a scientific attribution one.  Here's the exchange:

jstults

that climate change is anthropocentric

Do you mean that the anthro forcing is the largest? That it is one of the biggest? That it is significant (an equal among peers)?
Would results showing that it was measurable but not practically significant meet your desiderata?

• on February 5, 2010 at 11:19 pm | Reply Martha

  The meaning of anthropocentrism is more about how we dominate other species and generally overlook our duty towards others. Sadly, it has led to viewing the environment purely in terms of its value to humans; discussing climate change as if we are the only species impacted; and perceiving some members of our species as more important than others. As such, one can see that ‘anthropocentric’ has deeply philosophical and sociological meaning.
  I think you know that the ‘A’ in ‘AGW’ is usually taken to mean ‘anthropogenic’, as in ‘human caused’ — as in the overwhelming scientific evidence that the main driver of the current warming trend is C02 emissions from human activity.
  Did you have evidence to the contrary?
  “glaciers are melting… and not one entity that we trust with our money to look out for our interests as a species has any freaking clue as to why”
  There are entities that have a clue. They’re called climate scientists. They are not the entities that deal with your money.
  cheers
  
  
  • on February 5, 2010 at 11:33 pm | Reply jstults

    Martha said:
    

    The meaning of anthropocentrism is more about how we dominate other species and generally overlook our duty towards others

    I did not realize the philosophical connotations of that term. My personal ethical philosophy leans towards deontology, so I’m a bit familiar with the ethical concepts of ‘duty’. Would you say we have a ‘perfect duty’ towards the survival of other species?
    

    …as in the overwhelming scientific evidence that the main driver of the current warming trend is C02 emissions from human activity.

    I think the primary literature (and especially the recent findings about stratospheric water content) show that human CO2 emissions have effects with magnitudes similar to other forcings (at least on short time-scales, that’s one of the reasons it’s challenging to find the ’signal’ among the ‘noise’ in attribution studies). Can you point me at an attribution study that shows CO2 emissions are in-fact the primary driver?
    Thanks for addressing my questions and not flaming me.
    Cheers to you!
    
    
  • on February 6, 2010 at 12:08 am | Reply jstults

    “glaciers are melting… and not one entity that we trust with our money to look out for our interests as a species has any freaking clue as to why”

    Nice work on finding that! I wasn’t sure what that was from until I stuck it in to Google. You realize I wasn’t the one who said that right? So what about the claims about the AGW-disaster connections? Any thoughts on that?
    
    

on February 6, 2010 at 6:43 pm | Reply Martha

Hi,
1) Re. C02 Science. What primary literature are you speaking of? The overwhelming majority of published climate science says C02 is the main driver of the current warming trend. This science is linked and discussed all over this site and many other credible science sites. You aren’t the first visitor I’ve seen ignore it and I’m sure you won’t be the last.
I see you depend alot on Pielke. I encourage you to know that Pielke has been extensively debunked by competent scientists.
2) Re recent ‘findings’ regarding stratospheric water vapour. Is it Susan Solomon’s research you’re talking about?
http://www.realclimate.org/index.php/archives/2010/01/the-wisdom-of-solomon/
Romps and Kuang?
http://www.sciencedaily.com/releases/2009/04/090420121421.htm
These scientists are not challenging the core science of AGW.
Solomon, for example, wonders if changes in stratospheric water vapour might have helped keep the warming from being even worse in the past decade.
You know, no one is suggesting that all feedbacks are known or understood. Global climate and the carbon uptake cycle is a complex system. Climate knowledge is an evolving science. However, we don’t have to know everything to know alot, and the overwhelming majority of climate scientists say the evidence continues to agree with AGW.
The question is why you think they are wrong.
If you have evidence that they are wrong, please post it. This site has a thread for Challenging the Core Science.

 Of course, they wouldn't publish my reply!  So here's what you would have seen if they weren't such wieners:

C02 is the main driver of the current warming trend

I guess I'd still appreciate a pointer to that study; what I've been able to find so far is that radiative 'forcing' from CO2 has an effect with a size comparable to lots of other stuff (hence the problem of 'natural variability' sometimes masking our ability to measure the effect), and estimating the climate sensitivity to CO2 is still an open research problem.

...you seem to rely on Pielke a lot...

Well, when it has to do with his area of expertise (policy and disasters), and how his research was misrepresented, I do (you brought it up btw). I'm aware of 'name calling' on the internet directed at Pielke, but not any published rebuttals of his work (again, pointers to the lit appreciated).

Solomon, for example, wonders if changes in stratospheric water vapour might have helped keep the warming from being even worse in the past decade.

That was kind of my point, that effect was one of similar magnitude and of opposite direction to the CO2 forcing, so they 'canceled' (at least in the short term up to present).

the overwhelming majority of climate scientists say the evidence continues to agree with AGW

If by AGW you mean 'humans are affecting climate', then most 'skeptics' would agree too. If by AGW you mean some sort of moral prescription (you didn't answer my duty question btw) or particular policy package, then I'm not sure what science has to say about that.

...extensively debunked by competent scientists...overwhelming majority of climate scientists say...

Appeal to authority much? I wouldn't mind if you link to a search or two on scholar.google.com, and anyone that comes along latter would benefit from it too.

The question is why you think they are wrong.

You are assuming an awful lot about what I think. My question was about clarifying the meaning of 'anthropocentric' (I hadn't heard that particular phrasing before), and what exactly would 'count' in answering the challenge on this post. I wasn't trying to tilt at some windmill or pick a fight ('Challenging the Core Science' as you so reverently put it).

I thought, surely this was an easy slow one over the plate for someone to answer.  Link to the attribution study that shows CO2 is the 'main driver' (the other questions about duty and disasters were just for poking fun).  So I googled for 'quantifying anthropogenic forcing', here's the papers Google found.

The most recent is from 2006, and it is far from a slam dunk (I was surprised, I expected this to be a very straight-forward result based on what I'd seen on alarmist websites and the IPCC). The attribution method is simplistic (multiple linear regression optimal fingerprinting), and they scale the model outputs so that they fit the observations! Are you calibrating or attributing!? The other thing I thought was interesting was this snippet:

It is misleading and ultimately fruitless to suggest there can be some kind of ‘‘global’’ detection and attribution analysis capable of summarising the climate-related information-content of a dataset in a single estimation procedure.

I am a skeptic (which I was able to figure out by consulting this resource), and I am more optimistic about the fruitfulness of newer / better methods than that. I think the only way forward is a fully integrated Bayesian approach where you tackle the whole big uncertainty quantification / parameter fitting problem in a coherent manner. [Update: I re-read that paragraph and thought it needed a little more explanation. Attribution studies and model validation are closely linked to uncertainty quantification and data assimilation, the whole process is an inverse problem much like tomography, unfortunately those sorts of problems can be quite sensitive to the regularization (assumptions), so out-of-sample validation becomes even more critical. This is the sort of thing that leaves them open to Lindzen's 'prosecutor's fallacy' criticism:

Most importantly, Hasselmann (1997) noted that, for any formal detection and attribution procedure to get started, it is necessary to confine attention a priori to a relatively small number of competing explanations for observed climate change. If the number of allowed model-simulated signals, m, is too large, then it becomes increasingly likely that at least one signal will closely resemble a linear combination of the others, leading to a so-called ‘‘degenerate’’ estimation problem in which the data are insufficient to constrain the bi. We will consider cases up to m = 4 in this paper, which allows us to cover the main known drivers of recent near-surface temperature change: greenhouse gases, anthropogenic aerosols, solar variability and volcanic activity. We will find that, even with m = 3 or m = 4, many results become ill-constrained by the kind of large-scale data considered in this paper.

See, it's an inverse problem (with noise): Bayes leads the way! I won't touch their assumption about linear superposition of forcing effects (but, but it seems to work over this parameter range), this update is already too long.]

Any way, that's my trip report. If you have any interesting attribution studies that I missed link 'em in the comments. Thanks.

Joint Targeting Zen

Sometimes the most important part of the targeting cycle is deciding what targets not to engage to achieve the effects we want.

"It's possible (a strike) could be used to play to nationalist tendencies," Petraeus, head of the U.S. Central Command region, which includes Iran, said in an interview this week. "There is certainly a history, in other countries, of fairly autocratic regimes almost creating incidents that inflame nationalist sentiment. So that could be among the many different, second, third, or even fourth order effects (of a strike)." --Patraeus Says Strike on Iran Could Provoke Nationalism

The response of the British populace to The Blitz provides a good historical example of the kind of thing Gen Patraeus is describing.

Thirty spokes
Round one hub.
Employ the nothing inside
And you can use a cart.
Knead the clay to make a pot.
Employ the nothing inside
And you can use a pot.
Cut out doors and windows.
Employ the nothing inside
And you can use a room.
What is achieved is something,
By employing nothing it can be used.
--Tao Te Ching, 11

The Dayton Aerospace Cluster

This post is about the competitive advantage Dayton may or may not enjoy as an Aerospace Innovation Hub. This month’s Aerospace America has an interesting article titled “Changing aerospace cluster dynamics”, by Philip Butterworth-Hayes. A cluster is a group of similar or competitive businesses that locate in the same region, also known (in Ohio) as an innovation hub. The reasons given in the article for the increase of attractiveness of aerospace clusters are

• development of research institutions
• local and national government subsidies
• specialized work-force availability
• just-in-time production techniques

For these reasons, locating in a cluster gives a business a competitive advantage (see the Cluster Meta-study for more details). Like everything else in the world economy, the dynamics of clusters is changing and businesses in clusters are developing links with clusters around the world. Butterworth-Hayes claims that aerospace clusters in “North America and Europe will do less manufacturing and concentrate more on developing high-end integration and high-technology skills, outsourcing the labor-intensive manufacturing work to clusters in low-wage economies.” A familiar refrain to the folks of Dayton who have been struggling with the effects of the off-shoring (or down-southing) of manufacturing for decades. He continues, “This is exactly what seems to be happening. The key to transforming manufacturing clusters in North America and Europe to centers of aerospace excellence is access to innovation skills – corporate and academic. There is a new imperative for traditional manufacturing clusters to attract companies skilled in design and engineering, rapid prototyping, software development, part manufacturing, testing and research capabilities.”

Ohio’s governor has been in the news lately after designating Dayton an Aerospace Innovation Hub. The main focus of the Governor’s plan seems to follow the “transformation formula” described by Butterworth-Hayes. Grants to local research institutions and trade deals with foreign countries where the manufacturing will take place. So how good an idea is the Governor’s strategy? First it’s a good idea to find out what we’re starting with. Figures 1 through 5 come from the Cluster Mapping Project. They illustrate the relative size / employment / wages in traded clusters in the state of Ohio.

Figure 1 shows that Ohio has a significant and growing share of the turbine engine manufacturing employment in the US. The other good news from Figure 1 is the “Production Technology” cluster, relatively large and growing. The bad news shouldn’t be any surprise, “Automotive” and “Metal Manufacturing” are both large and shrinking.

---

Figure 1: State of Ohio: Specialization by Traded Cluster, 1998 – 2007

---

Figure 2 shows the part of the reason Ohio (and Dayton) is feeling some pain. The 2^nd and 3^rd largest employers are the shrinking ones. Production Technology is 9^th on the list and Aerospace Engines is 24^th.

---

Figure 2: State of Ohio: Employment by Traded Cluster, 2007

---

---

Figure 3: State of Ohio: National Employment Rank by Traded Cluster, 2007

---

Figure 4 shows that our Ohio workers in Aerospace Engines are well-paid compared to the national average. Our Aerospace Vehicles and Defense folks (2^nd on the list) are a paid a little less than the national average.

---

Figure 4: State of Ohio: Wages by Traded Cluster, 2007

---

The key political graph is Figure 5. Which sectors are creating the most jobs? The surprising (at least to me) big winner is “Transportation and Logistics”. Aerospace Engines added a tiny number of jobs over the time-period. Production Technology lost jobs over this period (Figure 1 showed the relative nationwide change).

---

Figure 5: State of Ohio: Job Creation by Traded Cluster, 1998 – 2007

---

For the Dayton region, the contribution to the Aerospace Engines cluster is basically GE in Cincinnati, that’s close enough to Dayton to be part of a virtual cluster (though their immediate outlook probably depends significantly on the fate of the F-136). The Aerospace Vehicles and Defense is basically WPAFB. So the cluster map seems to indicate that designating Dayton an ’aerospace innovation hub’ is at least somewhat reasonable. There’s a significant number of people in the region already doing the job. What about the rest of the “cluster checklist”?

• development of research institutions: DAGSI, AFRL
• local and national government subsidies: well $250k is a start…
• specialized work-force availability: GE’s big facility in Cinci, various small local shops in Dayton
• just-in-time production techniques: this isn’t really a characteristic of the region, just what efficient manufacturers are doing these days

This actually looks pretty good for Dayton. According to the formula, the well-educated, white-color professionals at the “top” of the value chain will be “innovating” and inking international trade deals while coordinating production at low-wage manufacturing clusters in Mexico and Asia. So what’s the problem? The problem for Dayton is Figure 5. Tiny growth in aero-manufacturing and defense / research jobs is a poor dressing for the sucking chest wound of automotive and metal manufacturing job losses. If we want Dayton to succeed, if we want all of Dayton to succeed, then Dayton needs a different approach than the commonly accepted “transformation formula” for aerospace clusters.

Why not coordinate manufacturing in a (relatively) low-wage cluster in Dayton’s hollow core? Wages couldn’t be as low as developing nations, but there is significant value in locating the engineers and the assembly lines in close proximity (see Gulfstream in Savanah or Cessna in Wichita), and coordinating the far-flung manufacturing is not without it’s significant costs/risks (see Boeing’s 787 experience). Don’t get me wrong, this won’t be a return to the glory days of gold-plated labor contracts. Labor in Dayton has to understand their competitive position and the value proposition that they bring. The availability of inexpensive foreign alternatives must temper their demands, but it doesn’t have to be a race to the bottom. The value of having design and manufacturing collocated means that wages can be higher than in the foreign low-wage alternatives, while still offering a competitive return on capital.

Dayton is different, generic transformation recipes don’t fit. Real innovation is finding a way to go after the top of the value chain (research / innovation / design) and provide opportunities for people in the middle and the bottom of the value chain at the same time. Not everyone can be a Chief, we need ways for Indians to contribute constructively too, or as Esrati put it: What do elevator attendants do in the information age?

Closing the Post-Shuttle-Gap Faster

Photo credit: Justin Ray/Spaceflight Now

SpaceX has a hangar full of rocket at the Cape. The latest rumors are that the new NASA budget from the President [Update: not rumors any more] puts the focus for near-term human space-flight (to LEO / ISS) on COTS-like providers. Based on SpaceX's claims (and their demonstrated commitment to develop, build and fly hardware), if this is true, the gap in US capability post-shuttle will actually be less (or eliminated) than if we continued to fund the Constellation boondoggle.

The latest AIAA news roundup on the topic follows.

Coverage continued on the new NASA plans expected to be announced Monday. Media coverage focused on the shift towards commercial users, but stressed that a fight is expected in Congress. The AP (1/31, Borenstein, Chang) reported that with the Obama Administration's expected proposal Monday that NASA buy launch services from private companies, there is "some concern...from former NASA officials worried about safety and from congressional leaders worried about lost jobs." Commercial companies see this very beneficial to the US, and executives "dismiss" safety concerns raised by various panels, citing the airline industry. According to the article, instead of the "established aerospace giants," it is the "newer space guard that brings some excitement to the field."

        The Washington Post (2/1, A8, Achenbach) notes, "NASA's grand plan to return to the moon, built on President George W. Bush's vision of an ambitious new chapter in space exploration, is about to vanish with hardly a whimper." Even though the upcoming budget will "effectively...kill" NASA's previous plans, "it remains to be seen whether Congress will accede to Obama's change in direction. Industry insiders expect a brutal fight in Congress." Meanwhile, White House spokesperson Nick Shapiro said Sunday the budget will show President Obama's "dedication to NASA. NASA is vital not only to spaceflight, but also for critical scientific and technological advancements."

        The Houston Chronicle (1/31, Powell, Berger) also reported that a "titanic political fight" is expected "as lawmakers from Texas move to protect jobs, astronauts and the Johnson Space Center." Congressmen from states with NASA centers "are drawing their battle lines to try to protect aspects of NASA's existing manned space program that account for jobs, payroll and subcontracts in their states." Scott Pace, head of the Space Policy Institute, said people "should not be surprised that the existing moon program has come to an end. ... Bush's budget decisions stressed the program; Obama's cuts drove it off the cliff."

        Meanwhile, Florida Today (1/30, Halvorson) reported Sen. Bill Nelson "says Obama shouldn't pump all of a $6 billion NASA budget increase into the development of commercial crew transportation services. Some of the increase should be put toward building a new heavy-lift launch vehicle." Nelson is concerned NASA will fall further behind in its development otherwise. "Investing in a super-sized rocket would leverage the $9 billion NASA spent on Project Constellation over the last six years" and "could increase the number of test flights at KSC, a move that would further make up for shuttle job losses

        According to BBC News (2/1, Amos), even though organizations like the Commercial Spaceflight Federation endorse the change, "there are major bi-partisan concerns in Congress about the proposed new path, and the impact it could have on jobs tied up with current Nasa programmes." AFP (1/31, Santini) reported than a White House advisor on Friday said, "Constellation is dead." However, "While Constellation is dead, it does not mean human space exploration is also dead," according to John Logsdon, former director of the Space Policy Institute, who points to plans to increase international cooperation.

        Bolden: NASA Will Deal With Non-Traditional Partners. Space News (1/30, subscription required) reported that NASA Administrator Charles Bolden, speaking in Israel, "said U.S. President Barack Obama's vision for NASA will elevate the importance of international cooperation." Bolden said NASA will be "vigorously engaged with Israel and other nontraditional partners as we move forward with the new vision."

        New Heavy Launch Vehicle Could Be Based On DIRECT. NASASpaceFlight.com (1/30, Bergin) reported, "The Michoud Assembly Facility (MAF) have confirmed they have almost enough External Tank resources to allow for one ET-sized 'In Line' Shuttle Derived Heavy Launch Vehicle (SD HLV) test flight and up to three Block I SD HLVs. The news comes as NASA managers insist the workforce should wait for official news, and not to be distracted by reports on Ares' demise." According to the article, a HLV "will be contracted out...moving to a multi-company effort led by Boeing," with partners like Lockheed Martin, Northrop Grumman, and United Space Alliance. A plan is progressing "at a healthy pace" towards a test flight "of what is now heavily confirmed as based on the DIRECT team's Jupiter-241 Stretched Heavy Launch Vehicle." The article noted that this plan, contrasted with the "sluggish" Ares V rocket, is reportedly "receiving significant support" from Congress members.

Chemistry Drives Convection

There's a short article with a neat set of flow-viz pictures and a movie on the APS site. It's kind of funny how research results get written up in the popular press.

As a demonstration of the chemistry "tail" wagging the fluid dynamics "dog," researchers have...

Although this chemically-induced convection was not a surprise, the convection patterns in the simulations were surprising. In some cases the convective plumes went only upward from the boundary, rather than creating a symmetric pattern that sent fingers both up and down from the mid-line.

It turns out that even without a chemical reaction, convection would be expected if the reactants had different diffusion rates.

The new simulations show that convection makes the reactions progress faster, just as stirring would.

I'm not sure how surprising any of that would be. Turbulent mixing is a well known way of increasing reaction rates, for instance in flame-holders and high-g combustors. It's also kind of funny how everything has to do with climate change:

This could mean that researchers have underestimated reaction and mixing rates in natural phenomena, like the Earth's mantle and supernova explosions, as well as for carbon sequestration, in which carbon dioxide is stored underground. Modelers have so far not considered how convection induced by the interaction of carbon dioxide with ground water may affect the long-term containment of this greenhouse gas...

Now that's the funding tail wagging the science dog.

FFT-based Abel Inversion Tutorial

I was putting together a little data post-processing tutorial for some guys in my lab, and it turned into a tour de force for that scientific computing swiss army knife, the FFT. So I figured I'd share it here.

The explicit recognition of its [the FFT’s] importance is one of the few significant advances in numerical analysis in the past few decades. Much has been learned about various aspects of computation, but no other recent discovery has so profoundly affected so many different fields as has the fast Fourier transform.
–R.W. Hamming [1]

This tutorial will demonstrate Gaussian convolution / deconvolution and Abel inversion of something resembling microwave interferometry data. The tutorial uses Scipy [2], but the concepts (as well as most of the function names and even the underlying FFT libraries) transfer directly to other environments (Matlab, Octave, etc). When in doubt, Google or Hutchinson [3] are your friends.

1 Gaussian Convolution / Deconvolution

The Microwave Interferometer beam is assumed to have a Gaussian distribution of intensity. The equation for a Gaussian distribution is

(1)

This means that the measurement at each spatial position is a “mixture” of effects from neighboring areas along the beam’s axis. To understand this, we’ll do some convolutions with Gaussians. First off, we need an array of points from a Gaussian, but we need them in a particular order. Here’s a Python function which does what we need:

def periodic_gaussian(x, sigma): 
    nx = len(x) 
    # put the peak in the middle of the array: 
    mu = x[nx/2] 
    g = gaussian(x, mu, sigma) 
    # need to normalize, the discrete approximation will be a bit off: 
    g = g / sum(g) 
    # reorder to split the peak across the boundaries: 
    return(sp.append(g[nx/2:nx], g[0:nx/2]))

The reason for reordering is because we’ll be using FFTs to do the convolutions, and the FFT assumes periodic data. Figure 1 shows the results for several values of the standard deviation.

---

Figure 1: Reordered Gaussians

---

Once we have the proper distributions defined, the convolution just requires an FFT of the theoretical curve, and the Gaussians, point-wise multiplication in the transformed-domain, followed by an inverse transform. The Python to do this is shown below:

dx = x[1] - x[0] 
sigma = [25.0, 50.0, 100.0, 150.0] 
g = [] 
cyl_conv = [] 
fcyl = fft(cyl) 
for i,s in enumerate(sigma): 
    g.append(periodic_gaussian(x, s*dx)) 
    cyl_conv.append(ifft(fcyl * fft(g[i])))

Figure 2 shows a theoretical phase shift curve for a calibration cylinder along with convolutions of that same curve with Gaussians of various standard deviations.

---

Figure 2: Convolutions of Calibration Cylinder Theoretical Phase Shift with Gaussians

---

Convolution is easy enough, but what we have in the case of an interferometry measurement is the convolved signal, and we’d like the ’true’ measurement. We can get this by deconvolution. This is simply division in the frequency space as opposed to multiplication:

# now try a deconvolution: 
sigma_deconv = [94.0, 100.0, 106.0] 
cyl_deconv = [] 
# use one of the previous convolutions as our ’measurement’: 
fmeas = fft(cyl_conv[2]) 
for i,s in enumerate(sigma_deconv): 
    cyl_deconv.append(ifft(fmeas / fft(periodic_gaussian(x, s*dx))))

Figure 3 shows the results for a couple test σ’s. Unfortunately, deconvolution is numerically ill-posed and very sensitive to noise, it tends to amplify any noise present. Unless we are very lucky and guess the exact σ and our measurement has no noise, we end up with ’wiggles’.

---

Figure 3: Deconvolution of σ = 6 ’measurement’

---

We can fix things by first ’smoothing’ with a small Gaussian, and then performing some test deconvolutions.

# now try a deconvolution with some ’regularization’: 
sigma_reg = 50.0 # about half our ’true’ sigma 
cyl_deconv_smooth = [] 
# use one of the previous convolutions as our ’measurement’, smooth it: 
fmeas = fft(cyl_conv[2]) * fft(periodic_gaussian(x, sigma_reg*dx)) 
for i,s in enumerate(sigma_deconv): 
    cyl_deconv_smooth.append(ifft(fmeas / fft(periodic_gaussian(x, s*dx))))

---

Figure 4: Deconvolution of σ = 6 ’measurement’, with smoothing

---

The deconvolution that uses σ > σ_true still displays overshoots even with the smoothing, it also becomes too ’sharp’. The calibration measurements can be used to find the optimal smoothing and sharpening σ for the particular interferometer / frequencies being used.

2 Abel Inversion

Many of the same tools presented in Section 1 can be extended to the calculation of the Abel inversion. In microwave interferometry we have a ’chord-averaged’ measurement, but we would like to understand the radial distribution of the quantity (phase shift / plasma density).

(2)

where F is our chord-averaged measurement, and f(r) is the radial distribution that we would like to recover. To calculate the inverse transform we first need an estimate of the derivative ∂F∕∂y. Like deconvolution, numerical differentiation tends to amplify noise (numerical or otherwise). To avoid this we can convolve the derivative of a Gaussian with the function we wish to differentiate and then perform an inverse transform to recover a ’smoothed’ numerical derivative. The derivative of the Gaussian is

(3)

A Python function that gives a properly re-ordered first derivative is

def periodic_gaussian_deriv(x, sigma): 
    nx = len(x) 
    # put the peak in the middle of the array: 
    mu = x[nx/2] 
    g = dgaussiandx(x, mu, sigma) 
    # need to normalize, the discrete approximation will be a bit off: 
    g = g / (-sp.sum(x * g)) 
    # reorder to split the peak across the boundaries: 
    return(sp.append(g[nx/2:nx], g[0:nx/2]))

In the same way that we used the Gaussian for smoothing we can use its derivative for differentiation. Here’s the Python to do that:

# calculate some derivatives: 
sigma_deriv = 0.5 # this can be small because we’ve already smoothed a lot 
cyl_deriv = [] 
for i,c in enumerate(cyl_deconv_smooth): 
    cyl_deriv.append( 
        ifft(fft(c) * fft(periodic_gaussian_deriv(x, sigma_deriv*dx))).real)

Figure 5 shows the results. The wiggles that we don’t manage to smooth in the initial sharpening stage become amplified by differentiation.

---

Figure 5: Derivative of De-convolved Functions

---

Once we have a numerical approximation for the derivative, calculating the integral is straight-forward:

def abel(dfdx, x): 
    nx = len(x) 
    # do each half of the signal separately (they should be the same 
    # up to noise) 
    integral = sp.zeros((2,nx/2), dtype=float) 
    for i in xrange(nx/2, nx-1): 
        divisor = sp.sqrt(x[i:nx]**2 - x[i]**2) 
        integrand = dfdx[i:nx] / divisor 
        integrand[0] = integrand[1] # deal with the singularity at x=r 
        integral[0][i-nx/2] = - sp.trapz(integrand, x[i:nx]) / sp.pi 
    for i in xrange(nx/2, 1, -1): 
        divisor = sp.sqrt(x[i:0:-1]**2 - x[i]**2) 
        integrand = dfdx[i:0:-1] / divisor 
        integrand[0] = integrand[1] # deal with the singularity at x=r 
        integral[1][-i+nx/2] = - sp.trapz(integrand, x[i:0:-1]) / sp.pi 
    return(integral)

The one minor ’gotcha’ is dealing with the integrable singularity in 2. We just copy the value of the nearest neighbor over the Inf, the integral will approach the correct value in the limit of small Δx (lots of points).

---

Figure 6: Inverse Abel Transform Results

---

The analytical inverse Abel transform of the theoretical phase shift curve can be easily calculated (Maxima [4] script shown in section A.2). The theoretical phase shift has the form

(4)

which gives the derivative as

(5)

and finally, the inverse Abel transform as

(6)

2.1 Centering

Often, despite our best efforts, the data are not centered. There is a method to estimate the optimal shift required to center the data using the FFT [5]. Figure 7 shows one of the “measurements” shifted by a small amount.

---

Figure 7: Shifted Measurement

---

The method of Jaffe et al. [5] is to maximize

(7)

Where L_nh is the Fourier transform of the data, and m_kΔ is a shift. Since the data is real, the real part of its transform is the even part of the signal. The shift which maximizes the magnitude of the even part is the one which makes the data most even (centered). A graphical analysis of the data can provide a range of shifts to try. The Python to perform this maximization is shown below (this relies on being able to index arrays with arrays of integers, which can be done in Matlab as well).

# a shifted measurement 
shift = nx / 10 # a ten percent shift 
cyl_shifted = cyl_conv[2][scipy.append(range(shift,nx), range(0,shift))] 
# find the shift needed to center the measurement (should match shift above): 
indices = [] 
# get range of shifts to try from ’eye-balling’ the measurement 
try_shift = range(nx - 2*shift, nx) 
for i in try_shift: 
    indices.append(scipy.append(range(i,nx), range(0,i)).astype(int)) 
even_norm = sp.zeros(len(indices), dtype=float) 
for i, cyclic_perm in enumerate(indices): 
    even_norm[i] = sp.sum(fft(cyl_shifted[cyclic_perm]).real**2) 
max_even_shift = try_shift[even_norm.argmax()]

Figure 8 shows that the objective function (equation 7) reaches a local maximum at the shift which was originally applied (shown in Figure 7).

---

Figure 8: Centering Objective Function

---

3 Calibration Model Fitting

The techniques presented in section 1 can be combined to fit a calibration model to interferometry measurements of a calibration cylinder (or any other object for which you have the theoretical phase-shift curve). The purpose of a calibration model is to provide estimates of parameters like the standard deviation of the Gaussian beam for use in deconvolution and inversion of subsequent test measurements.

There steps in fitting the model are

1. Center the data (2.1)
2. Find an error minimizing set of deconvolution and smoothing σs (these are the parameters of our model).
   1. Define a function of the data and the σs which returns an error norm between the theoretical phase shift curve and the smoothed, deconvolved data.
   2. Apply a minimization algorithm to this function. Choose an initial guess for the deconvolution σ which is likely to be smaller than the true σ.

The Python functions to accomplish this are shown below. One error function uses a sum of squared errors and the other uses a maximum error.

def smooth_deconv(s, x, data): 
    fsmooth = fft(data) * fft(periodic_gaussian(x, s[0])) 
    fdeconv = ifft(fsmooth / fft(periodic_gaussian(x, s[1]))).real 
    return(fdeconv) 
 
def cal_model_ls_err(s, x, data, theoretical): 
    fdeconv = smooth_deconv(s, x, data) 
    return(sp.sqrt(sp.sum((fdeconv - theoretical)**2))) 
 
def cal_model_minmax_err(s, x, data, theoretical): 
    fdeconv = smooth_deconv(s, x, data) 
    return((abs(fdeconv - theoretical)).max())

The error in the fitted model following this approach is shown in Figure 9. Both the least squares error and the min-max error methods recover the correct deconvolution σ.

---

Figure 9: Error in Calibration Fit

---

4 Further Improvements / Variations

Jaffe et al. [5] also suggest “symmetrizing” before calculating the Abel transform by removing the imaginary component of the transformed data. They also estimate the variance of the inverse Abel transform based on both the portion of the data removed by filtering (smoothing) and symmetrizing.

Rather than calculating the derivative necessary for the inverse Abel transform with a convolution, a pseudo-spectral method could be used. The trade-off being that very accurate derivatives could be attained at the expense of greatly amplifying any noise. A more practical variation, that is not susceptible to noise, would be to calculate the integral with a pseudo-spectral approach rather than the composite trapezoidal rule. This would give a more accurate integral, and bring the numerical and analytical curves in Figure 6 into closer agreement.

A Complete Scripts

A.1 Python

The complete Python script to do all of the calculations and generate all of the plots in this tutorial is available.

A.2 Maxima

The Maxima batch file used to calculate the inverse Abel transform of the theoretical phase-shift curve is available.

References

[1]   Hamming, R.W., “Numerical Methods for Scientists and Engineers,” 2nd. ed., Dover Publications, 1986.

[2]   Jones, E., Oliphant, T., Peterson, P., “Scipy: Open source scientific tools for Python,” 2001–, http://www.scipy.org

[3]   Hutchinson, I.H., “Principles of Plasma Diagnostics,” 2nd ed., Cambridge University Press, 2002.

[4]   Schelter, W., et al., “Maxima: a computer algebra system,” http://maxima.sourceforge.net

[5]   Jaffe, S.M., Larjo, J., Henberg, R., “Abel Inversion Using the Fast Fourier Transform,” 10th International Symposium on Plasma Chemistry, Bochum, (1991).