Best SNR

I feel like the discussion of the maximum of SNR (signal-to-noise ratio) will benefit my thesis, because my boring mathematical derivation can easily occupy some space to lengthen the number of pages.

To make things easy, I assume that the signal of the observed target can be described by a Gaussian distribution, namely,

I(\rho) = I_{0} \exp \left(-\frac{\rho^2}{2 \sigma^2} \right) ,

where I_{0} is some referenced intensity, \rho is the radius to the centre of the source, and \sigma describes its FWHM. The integrated intensity of the source within an aperture r is then calculated from

F(\rho) = \int_{0}^{\rho} 2 \pi r I(r) \mathrm{d}r.

After doing some simple algebra the result is simply

F(\rho) = 2\pi \sigma^2 I_{0} \left[1 - \exp \left(-\frac{\rho^2}{2\sigma^2} \right) \right] .

The SNR of the source is given by the following equation:

\mathrm{SNR} = \frac{F}{\sqrt{F + \pi \rho^2 f}}.

For simplicity, we assume the sky background remains constant across the targeted signal source, so f can be regarded as a combination constant value related to the sky background intensity and the readout noise of the CCD. To obtain the maximum SNR, we ought to find out when \mathrm{d~SNR} / \mathrm{d} \rho = 0 , which can be transformed to

 F' F + 2 f \pi \rho^2 F' - 2 f \pi \rho F = 0,

where F' = 2 \pi \rho I_{0} \exp \left(-\frac{\rho^2}{2 \sigma^2} \right).

With some algebra, this equation can be changed to

\left(k e^{-u} - 1 \right) \left(1 - e^{-u} \right) + 2u e^{-u} = 0 ,

where k = I_0 / f, and u = r^2 / 2\sigma^2 . Unfortunately, this is a transcendental equation so I cannot find out its solution analytically. Instead, let me use the numerical method.

Assuming the targeted source is very faint, which means that k \to 0, I obtain \rho \approx 1.6 \sigma. Conversely, if the source is extremely bright, i.e., k \to \infty, then we have \rho \approx \sqrt{2 \ln k} \sigma .

Anyway, the best-SNR radius as a function of k is plotted in the following figure. Note that I have converted the radius in FWHM by using \mathrm{FWHM} = \sqrt{2 \ln 2} \sigma \approx 2.355 \sigma .


In conclusion, using \rho \approx 1 FWHM looks to be a good choice.

Best SNR

Observation Run at KPNO (Mar 2017)

My colleague Ariel and I went to KPNO near Tucson for our observation run just right before the spring break.

Day 1 (Mar 21)

The flight from LAX to Tucson was ~ 8 am, so I had to wake up at ~ 5:15 am, and then quickly packed up my backpack. After eating some leftover bánh trôi in the fridge, I set off to catch the 6am Flyaway At the Westwood. Within half an hour, I arrived at LAX, where I met Ariel an hour later. The flight was smooth. On board I recognizes your destination – KPNO in distance, thanks to the super obvious signature just next to KPNO,  Baboquivari Peak. Of course all the domes appeared as small what dots merely.

We took a Uber to the KPNO in downtown Tucson, on campus of University of Arizona. Since the shuttle wouldn’t leave until 11 am, Ariel and I walked around the campus very briefly, but soon we decided not to continue because the temperature was scorching. Actually I’ve visited the campus three years ago during the DPS. However, it was at night so I didn’t see things clearly. It turned out to be that Ariel and I were the only two persons going to KP, and Ariel had to drive the shuttle.

The bicentennial Moon tree in UA’s campus. Last time I saw it in dark.

The drive was ~1.5 hrs. The weather atop KP was visually perfect, although the wind was pretty strong. Starving we were, we devoured quite a lot at the cafeteria. We then went to the 36″ dome, waiting for Flynn, who showed up ~2:30 pm and gave us some very brief instructions on operating the dome and the telescope. Basically we still remember every procedure because our last visit was three months ago only. Flynn mentioned that recently there have been several reports of bear and mountain lion witness around KP and asked us to be cautious.

Ariel and I were both exhausted, so we went back to the dorm and took a nap. After dinner, we met at the dome and started working. Ariel arrived earlier. When I got there, she was taking flats already. Having adjusted the focus of the telescope by taking images of θ Gem, we started observations. Our targets were all comets, SPCs and LPCs. We did plan to observe some unusual asteroids with TJ < 2, including several Damocloids, but eventually didn’t have enough time. During the early run, we unconsciously stuck the telescope because the pointing exceeded the minimum altitude limit twice. As a result, the telescope had to be manually slewed by me through a control panel upstairs. We have now been experts on this!

We observed a PCCP object YF982D7 through BVRI filters. I quickly calibrated all the images, and measured astrometry. No cometary feature can be seen from our coadded images. It finally was designated as 2017 EF12 by the MPC on Day 3, so the cometary feature wasn’t confirmed at all. I’m wondering if some observer has been fooled by bad seeing. Otherwise it’d be quite interesting and exciting, were it a comet; it has a low-eccentricity orbit, and the semimajor axis is small too.

Before daybreak, we attempted to observe several newly discovered comets, such as C/2017 E4 (Lovejoy), only to realise that the influence from the Moon was awfully serious, so had to give up. The observing run ended at ~5:40 am, when the twilight became too strong.

After taking a shower, I went to bed and felt asleep soon.

Day 2 (Mar 22)

I woke up several times around noon, and forced me to sleep. Finally I got up at 4 pm and went to have dinner (breakfast in reality). The weather had been overcast, but the cloud dispersed and revealed a crystal clear sky before sunset. I went to the dome and met Ariel there. Again she was taking flats.

The scene around Baboquivari in the south viewed from the dome.

As the twilight faded, we began the observing run with the synoptic observation about a variable star in the region of M42, which happened to be the only observed target tonight. As soon as we finished taking BVRI images of it, we saw the instruction of shutting down the dome, due to the strong wind and rising humidity. So we did. The humidity went up quickly and reached ~100% later. Ariel left the dome around 2 am, and I left around 4 am.

On the way back I felt the ambience a bit horrible. The smog was pretty strong; I could hardly see things beyond ~5 m even with a torchlight. The gust kept producing weird noise. Anyway, I successfully made my way to the cafeteria, and ate some hot soups there.

The spicy chicken soup.

Before I went to bed, it started raining. What a night!

Day 3 (Mar 23)

Again I woke up several times before the time I planned to get up. This is equivalent to suffering from jet-lag. Before I went to have my dinner, there were considerable amounts of clouds in the sky. But fortunate enough, they dispersed towards dusk. The temperature was also much colder.

We still had plenty of time after having finished taking flats. To kill the time I was busy calibrating images from Day 1 and performing astrometry. Of course I couldn’t finish the all before the sky got dark enough. So we adjusted the focus. I should point out that seeing was extremely awful. However we changed the focus, images of field stars looked fat. Anyway we began the observing run with 128P, which somehow appeared fainter than expected, perhaps due to the fact that the twilight was still there a bit.

Around 9:30 pm, a staff (I forgot his name) brought some visitors into our dome, showing them “what astronomers are doing”. It didn’t last for long, and they left within ~15 min. We then concentrate back on imaging comets. During the run we did attempt to observe Damocloid 2015 RK245, but we couldn’t spot it in individual images.

Before the end, we observed comet C/2017 E4 (Lovejoy), which appeared super bright. Judging from the morphology, it’s a gassy one. It moved pretty fast too. We obtained 10-ish images, which did give us enough SNR of it. We finally challenged the telescope a bit by observing another bright comet C/2017 E1 (Borisov) — then the twilight had been quite strong already. We did get it, but way less impressive than C/2017 E4, maybe because of the bright sky background. The session ended at ~5:45 am. The slender crescent moon looked amazing in the eastern sky.

Day 4 (Mar 24)

Maybe I’m really aging; I again woke up first at 12pm, and then woke up every hour. I dreamed a lot, probably due to which I didn’t sleep well. Anyway, the weather for this observing run appeared to be the best since we arrived at KP; the humidity was low enough, <~ 40%, basically windless, and the temperature was pretty steady, ~9 ℃. As a result, seeing was way better than on Day 3.

We began the observing run from comet C/2010 U3 (Boattini), in Perseus, followed by the synoptic observation of the M42 field. Comet C/2017 C1 (NEOWISE) was observed successfully. Originally we expected it to be faint, since no one had observed it since Feb, but it turned out to be a good one. It was readily visible in individual images, moving fast. The synoptic observation also required us to observe a field containing ρ Oph, so we wasted some time on this.

My night lunch for Day 4. The highlight is that spoon. Because I couldn’t find a spoon, I had to make one from foil. It looked fancy, even like from the antiquity, but actually functioned awfully…

The session was ended with C/2017 E1 (Borisov). Compared to Day 3, we were able to get two more images before the twilight was too bright.

Day 5 (Mar 25)

I slept a bit better today, at least my first wakeup was close to 1pm. I could clearly remember what I dreamed during the sleep — my family, and my calligraphy teacher 張老師, and his master pieces. He suffered from a stroke last Sept, resulting in weak disability of his right arm. Since then I had had quite a lot of dreams, in which he had convalesced and could write things as before. But a good thing is that I’ve heard that he’s getting better and better.

Pano pic of KP, taken by my phone.

The late afternoon witnessed a considerable amount of clouds. But thanks to them, I enjoyed a spectacular sunset. The clouds were like burnt, against a violet background of the sky. As time evolved they changed their shapes. I was simply standing outside and looking at the sky, doing nothing else. As one of my favourite things in my childhood which I did a lot was looking at the sky during sunsets, the scene certainly brought me back, just like yesterday once more. I thought I were only a little kid at a moment.

The awesome sunset view, taken by my phone.

I went back into the dome after the show faded. The clouds were great decorations during the sunset, but form a bad thing for our observing run. The dome was shut down because of them. We had to wait till ~11 pm until the sky turned clear. Our first target was 315P, which displayed a fan-shaped broad tail. A super bright comet we observed was C/2015 ER61 (PANSTARRS), which was initially discovered as an asteroid by Pan-STARRS, but later a cometary feature was found.

This day’s observing run ended again by observing C/2017 E1. Earlier, our advisor Dave emailed us asking whether we’d observed it; it might be a potentially disintegrating candidate due to the rotational instability. Frankly I didn’t think about this point, but just thought that we shall observe it to improve the orbit; the current solution was still pretty nasty, because we saw a significant discrepancy between the observed and predicted positions. Anyway, this time we put more efforts to this comet, and we were able to obtain seven R-band images on it.

Before going back to the dorm, we saw the very elegantly thin crescent Moon rising in the east.

Day 6 (Mar 26)

I finally managed to adapt to the timezone for astronomers, which is good. But unfortunately it was the last day of our observing trip, which is bad.

We planned to end the whole KP observing run by observing C/2017 E1 again. However, for some reason an attempt image showed that the image quality couldn’t be even more nasty, and thus we decided to observe C/2013 C2 (Tenagra), which was predicted to be bright enough, provided that the prediction didn’t go wrong. It hadn’t been observed for almost a year, based on the MPC. So we did. But we couldn’t see anything obvious, so at least it wasn’t that bright. I haven’t yet reduced the data. Maybe something will show up, hopefully.

It’d been pretty bright at the time we shut down the dome. I knew that Venus should be visible, so on the way back we searched for it for awhile. Indeed we managed to spot it by naked eye! It just passed the inferior conjunction with the Sun a day ago! This definitely becomes my record of seeing Venus closest to the Sun during non-transits by naked eye. Alas my camera was left in my office, a bit pity.

Before I left the dome, I drew a super ugly graffiti of a man’s head. It used to be a maze but then I thought it was kinda boring. Some observers might feel asleep during observing runs, so putting something on like this horrible may help them stay awake. I’m considerate, ain’t I?

I got up at noon the next day. Before we left the mountain I grabbed something for lunch. After checking out, Ariel and I drove back to downtown Tucson at 2 pm, and then took a Uber to the airport. Our flight was ~5:15 pm and landed in LAX an hour later. So the whole trip came to an end, but I still have to reduce and measure all the images from Day 6.

Observation Run at KPNO (Mar 2017)


Today saw my first success of establishing an orbital linkage for a SOHO comet, SOHO-3165 , with observations from three consecutive apparitions. I have never done this work on my own — all the previous similar works I’ve done were either verifying Rainer’s solutions, or simply lucky because of lacking influential nongravs. SOHO-3165 doesn’t seem to be an easy one, as neither Bill nor Karl managed to get a solution. Bill managed to get a linkage with the first two apparitions, and third one was left unsuccessful with huge residuals. Therefore I was curious if I could solve this problem.

I haven’t done this job for long, perhaps over three years. As a result it took me a while before I could get onto the right track. A few tweaks of starting conditions enabled Aldo’s Exorb to get a decent orbit for the first apparition. In this step, a parabolic orbit had to be assumed. I then forced the program to solve an orbit for semimajor axis a = 3.05~\mathrm{AU}. With observations from the second apparitions the orbital solution was quickly refined. The trick was to fit for the first observation from the second apparition initially. Once the RMS decreased, more subsequent data could be fed in.

If there were no nongravs, the third apparition could be linked without much difficulty. However, this is not the case this time, which exactly complicated the computation work. In fact there is a huge offset in positions in the first observation of the third apparition (which was the only observation from the third apparition at this step). So there must be a nonzero transverse nongrav parameter. I tried A_2 = 10^{-8} \mathrm{AU ~ day}^{-2}, no success, and even worsened the solution to the previous two apparitions. Yet the program turned the optimised A_2 into a negative number, from which I immediately realised that the third apparition was postponed. Then -10^{-10}~ \mathrm{AU~day}^{-2}. Success! The program quickly shrunk the RMS of the fit. My last step was adding all the remaining observations and let the program find the best-fit solution.


Yet seems like there is something wrong with the code about completing iterations. The code somehow continuously looked for but never succeeds in yielding a final convergence. I therefore manually paused the program. Luckily because the code was generally looping around small RMS. Although the above solution may well not have the smallest RMS, it has one close to already.

The motivation of this writing is simply to provide me as a reminder in the future about how I can obtain an orbital linkage when nonzero nongravs are presented.


Further to My Previous Blog, rgd RVSF

I got in touch with Dr. Nalin Samarasinha at PSI, who wrote the codes for all of the special cometary processing filters at PSI website. He has beeing patiently answering every of my questions so I now can understand what his RVSF code is basically doing, though I know very little of FORTRAN. The discrepancy I spot between the source code and the explanation file is indeed a typo.

I added a new keyword to my IDL routine so that the user could have the option whether they like to switch on or off the function of sub-pixel sampling of the image prior to the function of the filter or not. Through some simple tests I realized that the differences between sub-pixelization and non-sub-pixelization of the image would not be very obvious in terms of visual inspection. However, the speed matters considerably, especially when the input image size or the part of regions selected is quite large — the processing time, in this case, could be quite consuming. So it’s a good idea to preview the enhanced image without sub-pixel sampling, not only would you save a great amount of time, but also this would give you a basic knowledge about if your kernel parameters are appropriate, and finally, to process the image by sub-pixelization with reasonable kernel parameters.

I couldn’t discern any difference in the enhanced images with or without sub-pixelization, however, a subtraction between the two clearly reveals what is previously hidden behind. See the following image.

The difference between non-sub-pixelization and sub-pixelization sampling. Linear stretch.
The difference between non-sub-pixelization and sub-pixelization sampling. Linear stretch.

Additionally, I’m very delighted to find that the resulting enhanced image processed by means of sub-pixelization shares great similarity with the one processed through Nalin’s code or PSI’s online tool. So I think my routine is quite successful.

At last, before the end of this update, I need to confess that I had mistakenly regarded the author of the source codes to be Padma Yanamandra-Fisher. I had been confused by their names, all look rather long for me… ><

Further to My Previous Blog, rgd RVSF

Testing My Radially Variable Spatial Filter Code

I think I need to update my blog; it has been so long since my last update… Also should I prove that I’m still alive here.

I saw a message in the comets-ml posted by Martino Nicolini that the PSI website has released a web tool for processing cometary images, including azimuthal median/average/renormalization filters, 1/rho coma model division, and, quite unfamiliar to me, radially variable spatial filter ( I attempted to process an image taken by HST regarding comet C/2012 S1 (ISON) in May 2013 with the online tool, however, I had difficulty in retrieving the enhanced or processed data. Weirdly the size of the to-be-downloaded file was always 0 KB, obviously problematic. With consideration that the network is not permanently available for me, it sounds absurd and waste-of-time for me to wait for accessible network before performing specially processed cometary images, and therefore I decided to write codes for my own purpose.

So I did by taking reference to the explanation file of the enhancement techniques, understanding the global idea of how the radially variable spatial filter works. Realizing the algorithm of the filter was proved easy, and it didn’t take me long to accomplish the IDL routine.

I did some tests with the image CometCIEF_test.fits provided in this page. The following was generated with kernel A = 4.0, B = 4.0, N = 0.4. It looks correct anyway, quite similar to the appearance in the tutorial file.

RVSF test imageHowever, I found my result would look somewhat different from those presented in the tutorial file if the kernel size was smaller — probably the scaling plays a role there, yet anyway my outcome would look less detailed. The corresponding FORTRAN source code of the filter in the PSI page seems to has a typo, which I have already reported to Padma Yanamandra-Fisher. I’m still comparing my codes against the PSI’s…

Testing My Radially Variable Spatial Filter Code



最近想徹底解決掉遺漏許久的彗星塵埃粒子在力場中的問題,如今所寫出的IDL代碼已經可以完美地解決了\beta = \frac{F_{rad}}{F_{grav}}< 1的情況。如下圖所示,模擬C/2011 L4 (PanSTARRS)彗星在今年4月1日的部分塵埃尾巴情況。

Synchrone and syndyne grids for the comet
Synchrone and syndyne grids for the comet

解決這個問題的本質其實就是不斷地計算不同塵埃粒子的Kepler’s Equation。由於粒子大小不同,釋放時間也不同,其mean anomaly和偏心率也會不同。計算方法參考的是Probstein和Finson於1968年所發表的一片論文《A theory of dust comets. I. Model and equations》。不過論文僅針對了軌道是拋物線運行的彗星的情況。但只要按照同樣的思路,橢圓軌道、雙曲線軌道的彗星情況同樣可以迎刃而解。然而論文中只討論了\beta < 1,對於大於及等於1之情況卻未有照顧。不等式取等號時對應塵埃粒子將以勻速直線運動方式運動,而大於時則是受到中心斥力。於是我很自然地聯想到,在中心斥力力場作用下,塵埃粒子的運動方程是否仍然有類似於中心吸引力的類似形式。


\ddot{r} - r\dot{\theta}^2 = \frac{\mu}{r^2}——(1)


其中\mu = \left ( \beta - 1 \right )GM_\odot 。首先先推導軌道類型。令u = \frac{1}{r}代入(1)式,並結合(2)式,可以得到如下一個二次非齊次線性微分方程:

\frac{\mathrm{d}^2 u}{\mathrm {d}\theta^2}+u=-\frac{\mu}{h^2}

可得其通解爲u=\frac{\mu}{h^2}\left [ e\cos \left ( \theta - \omega \right ) - 1\right ],其中e和\omega爲待定常數。將其表達回r,得到:

r=\frac{h^2/\mu}{ e\cos \left ( \theta - \omega \right ) - 1}——(3)

由解析幾何知,這是一條以外焦點爲原點的雙曲線。可記得Rutherford Scattering,被散射粒子正是沿着雙曲線軌道運行的!

我們定義q = r_{min} = \left | a \right |(1+e),a爲對應雙曲線之半長徑,爲負值。剩下的工作是找出r和\theta關於時間t的含數。


\ddot{r} - \frac{h^2}{r^3}=\frac{\mu}{r^2}——(4)

由於\ddot{r} = \dot{r}\frac{\mathrm{d} \dot{r}}{\mathrm{d} r},故(4)式可以化爲:

\dot{r}\mathrm{d} \dot{r}=(\frac{h^2}{r^3}+\frac{\mu}{r^2})\mathrm{d}r




對於r=q時,有\dot{r}=0,代入(5)可以解得K_1 = \frac{\mu}{\left | a \right |},即

\dot{r}^2=-\frac{2\mu}{r}-\frac{h^2}{r^2}+\frac{\mu}{\left | a \right |}——(6)

注意到v^2 = \dot{r}^2 + \dot{\theta}^2 r^2 ,對(6)變換下形式,於是立即獲得了一個與活力公式十分相似的一條描述粒子速率的一條公式:

v^2 = \mu (\frac{1}{\left | a \right |}-\frac{2}{r})——(7)



\mathrm{d}t = \frac{r\mathrm{d}r}{\sqrt{\frac{\mu}{\left | a \right |}}r^2 - 2\mu r - \mu \left | a \right |(e^2 - 1)}


\sqrt{\frac{\mu}{\left | a \right |}}\mathrm{d}t=\frac{r\mathrm{d}r}{\left | a \right |\sqrt{(r - \left | a \right |)^2 - a^2 e^2}}

我們令r - \left | a \right |=\left | a \right |e\cosh{F},則有:

\sqrt{\frac{\mu}{\left | a \right |^3}}\mathrm{d}t = (e\cosh{F}+1)\mathrm{d}F


\sqrt{\frac{\mu}{\left | a \right |^3}}t = e\sinh{F}+F+K_2

廣義Kepler’s Equation雛形初現!!上式中K_2爲一積分常數,用同樣的方式,考慮r = q時,選擇t=0爲過近日點之時,則可得到K_2 = 0 。於是廣義Kepler’s Equation爲

\sqrt{\frac{\mu}{\left | a \right |^3}}t = e\sinh{F}+F——(8)

可見,在斥力下的方程與引力下的方程形式上非常接近。仿照真正意義上的Kepler’s Equation,我們且將F亦稱爲eccentric anomaly(偏近點角)。這是個超越方程,只要將F解出,根據

r = \left | a \right |(1 + e\cosh{F})——(9)

則可求得r。r既得,\theta亦可根據(3)式求出。不過此處我們想看看\theta與F之間的關係爲何。爲簡單見,令 f = \theta - \omega

將(9)式代入(3)式,我們可以得到\cos{f} = \frac{e + \cosh{F}}{1 + e \cosh{F}}。利用半角公式2 \sin^2 {\frac{f}{2}} = 1 - \cos{f},得:

\sin{\frac{f}{2}} = \sqrt{\frac{(e-1)\sinh^2 {\frac{F}{2}}}{1+e \cosh{F}}}

類似地,利用2 \cos^2 {\frac{f}{2}} = 1 + \cos{f},可得:

\cos{\frac{f}{2}} = \sqrt{\frac{(e+1)\cosh^2 {\frac{F}{2}}}{1+e \cosh{F}}}


\tan{\frac{f}{2}} = \sqrt{\frac{e-1}{e+1}}\tanh{\frac{F}{2}}——(10)



Solution to IDL Memory Allocation Problem

Long time since the last update of this blog.

I’ve confronted with the following error when trying to processing a large array made up of SECCHI COR2 data, and it had been extremely painstaking before the solution was found. Many of the solutions to the problem told you to simply switch a platform, i.e. from Win 32-bit to Linux, in that the problem inherits from discontiguous distribution of  memory allocation. For a better understanding, let’s suppose you get 1000MB free space in total. However, it may well be that summing up several, i.e., five 200MB discontinuous free space gives you this value. Yet if you want to use 500MB at a time in IDL, the following error warning message pops out:

   % Unable to allocate memory: to make array.
      Not enough space
   % Execution halted at: $MAIN$

I get an excellent alternative way to overcome this problem partly after reading guides from Coyote’s Guide. It works for me! At least I can reduce a median background from the daily stack.

When working with IDL after that aforementioned error message, exit the Workbench first. The next step is to find out a file entitled idlde.ini in one of the files storing IDL. Then use some text editor to open it. It may read like following:


Make sure that {VM_DIR} remains unchanged; it’s the path to Java JVM on your machine and malfunctions may occur otherwise. Change the lines:


Restart your IDL Workbench and you may see difference. Yet please don’t get too exhilrated. If you still use an oversized array, this may not have your problem resolve. Just never be too careful when running IDL in Win 32-bit machine. Frankly speaking it’s better recommended to running IDL in non-Win machines, on which It’s said that more powerfulness of IDL will be functioned, but I can’t tell more as lacking such experience.

Solution to IDL Memory Allocation Problem