ZiHao

On 11/10/2023 at 05:19, athornett said:

Hi Zihao,

Thanks for getting back to me.

I am VERY excited = I have compared my data plots at different galactic longitudes with yours and there is close match for areas I have so far been able to study. Wow!

I would love to go the next step and create the same plots as you have done labelled "Milky Way Structure" - I don't mind if your scripts are not perfect! Any chance you could PM me with info on how to do it? Perhaps I could help polish them up a bit, if that would help in any way?

By the way, although I live in England, my wife is from Penang - so its definitely a small world and wonderful to have folks from Malaysia on this forum!! I am an associate member of the Astronomical Society of Penang, who are a wonderful group doing marvellous outreach.

Andy

What a coincidence! I am not a member of the society but I was a participant of their astrophotography competition in 2019.

The first step would be fitting multiple gaussians to the spectrum, you can do this in either the Origin software (https://www.youtube.com/watch?v=WX_X3XlbKRA) or the python script below where you need to manually enter the guesses for the parameters of the gaussians.

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
import scipy as sp
import scipy.signal

def gaussian(x, A, x0, sig):
    return A*np.exp(-(x-x0)**2/(2*sig**2))

def multi_gaussian(x,*pars):
    total = 0
    offset = pars[-1]
    i = 0
    for j in range(int(len(guess)/3):
        single_g = gaussian(x, pars[i+j], pars[i+j+1], pars[i+j+2])
        total+=single_g
        i+=2
    return total + offset

para = []
data = #yaxis data
data = sp.signal.medfilt(data,5)
velocity = #xaxis data
guess = #guesses for parameters
    
popt, pcov = curve_fit(multi_gaussian, velocity, data, guess)

j = 0
for k in range(int(len(guess)/3)):
    para.append(popt[j+k+1])
#save the parameters

plt.plot(velocity, data, 'b')
plt.plot(velocity, multi_gaussian(velocity, *popt), 'r')
plt.show()

The second step is to obtain the x-y coordinates in the milky way structure using the mean value (which refers to the radial velocity) from the parameters.

Vo = 220
Ro = 8.5
x = []
y = []

l = #longitude in radians
para = #parameters

for a in range(len(para)):
    Vr = para[a]
    Rs = Ro*Vo*np.sin(l)/(Vo*np.sin(l)+Vr)

    if Rs**2-Ro**2*(np.sin(l))**2 < 0:
        continue
        
    rpos = np.sqrt(Rs**2-Ro**2*(np.sin(l))**2) + Ro*np.cos(l)
    rneg = -np.sqrt(Rs**2-Ro**2*(np.sin(l))**2) + Ro*np.cos(l)

    if rpos > 0 and rneg > 0:
        continue

    elif rpos < 0:
        continue

    else:
        x.append(rpos*np.sin(l))
        y.append(-rpos*np.cos(l) + Ro)

The final step is superimposing our data on top of the galactic structure plotted using the parameters provided from [11].

ro = [3.48, 3.48, 4.90, 4.90]
theta_o = [1.57, 4.71, 4.09, 0.95]
k = [4.25, 4.25, 4.89, 4.89]

for i in range(4):
    theta = np.linspace(theta_o[i], theta_o[i] + 2*np.pi, num = 500)
    r = ro[i]*np.exp((theta - theta_o[i])/k[i])
    xm = r*np.cos(theta)
    ym = r*np.sin(theta)
    wains = plt.plot(xm, ym, 'k')
plt.show()

18 hours ago, athornett said:

This is just what I need to help me with my hydrogen line radio observations - particularly your atlas of plots at various galactic longitudes gives me something to compare = I am on holiday at present but as soon as I get home I will compare I have so far with yours.

do you know of any more comprehensive resource that shows expected plots and graphs?
also what software you use to collect data, process and plot it and then superimpose your results on galactic structure?
I am currently using ezRA suite from Ted a line in SARA.

thanks for posting!
Andy

I am writing my own python code for data acquisition, calibration and representation. Plotting the data is done with matplotlib and some formulae that you can refer to here. I will try to take a look at my code again as it doesn't look very user friendly for now, and if I have the time, I will publish it on github.

4 hours ago, Victor Boesen said:

Excellent work ZiHao!!! Good to see you're still working on H-line stuff I haven't had much time to do the same lately, but hopefully I will in the future. I'm currently working on refurbishing the ground station at my university. If I succeed I should have a 3m dish to do radio astronomy with

Victor

Thank you! Not too sure when I can start a new radio astronomy project since I have started university recently. Little space to work with in campus, but I will try. All the best with the refurbishment, and I look forward to the 3m dish. It will open up many opportunities for challenging radio astronomy projects!

The 21cm hydrogen line is a spectral line produced by atomic hydrogen. Due to the abundance of hydrogen atoms and the nature of radio waves to penetrate dust clouds in our galaxy, we can study the milky way by observing the doppler-shifted spectral line.

Construction of the geodesic parabolic dish developed by JA6XKQ and 1.42Ghz cantenna for the hydrogen line project:

For the introduction of the design, please read this article [1] written by JA6XKQ. For those who are interested in the details of the derivation and surface accuracy of the structure, you can refer to [2]. Another important thing is the supporting structures such as adding supporting rods behind the dish or additional inner meshes to prevent the heavy feed from deforming the structure. See [3] and [4].
This is a 1.4m f/D 0.345 dish. Three L shape aluminium bar as supporting structure is attached to a stainless steel bowl. The same material is used for the feed supports and the reflector is made from galvanized wire mesh. This picture is taken months ago and the mesh is now covered with rust, I plan to apply rust converter to the surface and spray paint at the end. An alternative is to use stainless steel mesh at the beginning, but the cost will be high. If you are interested, please refer to the excel sheet at [5] to obtain the specific dimensions.

Cantenna's diameter: 14.5cm, length: 18cm, probe's distance from the back (1/3 of the length): 6cm, probe's length: I started with 0.25*λ, and it is trimmed until a good SWR is obtained with a NanoVNA. These dimensions are based off from the design from Figure 3 at [6] If you want to determine the diameter based on your desired edge taper, you can see page 3 at [7]. A N type female connector is used to attach the copper probe. No soldering needed as there's a hole where the probe can be slid in tightly. To compute the dimensions of the cantenna and a choke ring, refer to the excel sheet from SETI League at [8].

Special thanks to my dad's friend for providing material for the dish and fabricating the feedhorn and its bracket.

Observations and results:

Spectra below are along the plane of our galaxy from galactic longitude 0 to 360 degrees at latitude 0 degrees and they are the combination of multiple gaussian curves and the corresponding radial velocity of the individual peaks are then used to construct the structure and rotation curve of the milky way. Radial velocities are corrected to the local standard of rest (LSR). Individual spectrum is multiplied by a constant obtained from the calibration using the S7 region's peak brightness temperature value from the LAB survey, at [9].

Observations for the hydrogen line map from declination -60 to 60 degrees are done by pointing the dish at a specific declination and different regions of the sky will be captured thanks to Earth's rotation. This is a technique known as drift scan. Hydrogen column density is calculated by integrating the area under the spectrum from radial velocity -150km/s to 150km/s and times 1.82x10^18, given at [10].

Milky way structure. Data seemed to only match with the Perseus and Carina-Sagittarius arms. Most of the points are in between the Perseus and Norma arm, this is probably due to unresolved peaks in the spectrum from smaller radio telescope. Parameters for the spirals arm plot can be found on page 7 at [11].

Lastly the rotation curve fitted with a logarithmic curve. The rotation speed is obtained by choosing the most red-shifted and blue-shifted radial velocities in the spectra obtained from Quadrant I and IV respectively. The rotation speed are lower compared to the published data from Clemens (1985) attached, because the most red and blue-shifted peaks are often weak and difficult to detect with good SNR and hence the chosen peaks will be not be the maximum and minimum radial velocities. But the general trend can be seen clearly in the graph, a flat curve as you go further away from the center of milky way, suggesting that something is giving them extra speed and most of the mass might not be concentrated in the center. Thanks for reading!

References:
[1]http://www.terra.dti.ne.jp/~takeyasu/Geodesic_Parabola_Antenna_2_1.pdf
[2]https://f4buc.pagesperso-orange.fr/parabole_geodesique2.htm
[3]http://www.terra.dti.ne.jp/~takeyasu/PhotoGallery.html#new_photo
[4]https://www.youtube.com/watch?v=od2SEhq8DR4&t=11s
[5]http://www.terra.dti.ne.jp/~takeyasu/
[6]https://www.semanticscholar.org/paper/A-compact-radio-telescope-for-the-21cm-line%3B-za-21-Saje-Vidmar/e5baaf7450e0aea72e7926e1bdd75d3175f23424
[7]http://www.w1ghz.org/antbook/chap6-3.pdf
[8]http://www.setileague.org/hardware/feedchok.htm
[9]https://www.astro.uni-bonn.de/hisurvey/euhou/LABprofile/
[10]https://www.cv.nrao.edu/~sransom/web/Ch7.html
[11]https://iopscience.iop.org/article/10.1086/501516/pdf

11 hours ago, Steve714 said:

Great 21cm work! Also from Sussex (East), I thought some measurements at ~3.3GHz on usual suspects (thermal (Moon)and supernova remnant(Taurus A)) might be of interest. This drift scan was made using a 1.9m diameter dish and LimeSDR receiver.  The moon was 1.4 degrees lower in the sky than the dish elevation which was set for Taurus A. Happy to provide details of system (LNAs, filters, software etc).

 

Nice observation. I think it's quite rare to see someone observing in the S band so I am interested in the details of the hardware and software used! 

12 hours ago, SteveBz said:

And what do you get for R, X & Y?

Using rv=11,
R=11.9, r=3.47, X=-0.603, Y=3.42
X=rsinl  Y=-rcosl
Add 8.5 to Y if you want the galactic centre to be (0,0).

 

3 hours ago, SteveBz said:

Yes it's here.  Thanks for looking.

Steve.

debug(ra=94.1,dec=21.1).json 131.77 kB · 1 download

Alright, I tried plotting the results just now and I have the same results, rv=11. This is interesting, are your rv for spectrum of other locations off by few kilometers as well? If yes, then I think it will be good to check with @Victor Boesen LSR correction and compare the results.

7 hours ago, SteveBz said:

Hi ZiHao,

I think you're right, it is the LSR.  However, I can't see how it's wrong.  Here's a plot:

My values are in the red box.  rv = 11 km/s and R is negative.  Obviously that is not going to work out well.  l = 190 deg is quadrant 3 and should have positive X and negative Y. The opposite is true.

The corrected LAB value is 6 km/s.  The LSR and planetary compound correction value is -41 km/s, quite large in comparison.  The orange curve is the Gaussian-fitted curve.  If anything it's skewed to lower values.

I started by calculating the corrections independently (-29 km/s and -12 km/s) and then togther by feeding the planetary correction onto the LSR calculation.  Both give -41 km/s.

Noty sure what to do next.

Steve.

I think there's some mistake with the formula that you are using, I tried plugging in rv=6.18, as well as rv=11, I get positive R. You can refer to the notes I attached earlier or the equation you corrected in the pdf. Is it possible for you to send the raw data for this spectrum, l=190? I will try to correct it to LSR later and see if it matches with LAB survey.

10 minutes ago, SteveBz said:

Here I've replaced the marker with the quadrant number for debugging purposes.  You can see that some of it fine and some is very strange.

 

Hmm...Are the LSR corrections correct? I try substituting the values from the LAB survey, they seems to give positive R. If the corrections are correct, then I guess you have to reject the values obtained near l=180, from the earlier pdf paper you shared, it seems like there is no data points around that area also.

1 hour ago, SteveBz said:

Another thing is that the approximation we make around differential rotation (ie V0=V) also seems to be wrong.  The actual answer for Q1 & Q2 is V=V0, but the answer for Q3 & Q4 seems to be V = -V0, because it's roating towards us.  That's why @ZiHao got the answer he did.

Furthermore, in areas where V0*Sin(l) ~ Vr, then the denominator has a very high error rate and is unreliable, especially, it seems to me, where 'l is near 180.  Sin(l) becomes very small and the R calculation diverges.  What do you think?

I don't think you should worry about the sign in this case. If you want to substitute V=-V0 into the equation(the one that differs with a negative sign) for Q3 and Q4, we will still get the same equation. You can also try playing with the signs in the radial velocity formula and the patterns of receding, approaching, etc in the quadrants will match with what we expected.

At locations for l=0,90,180,270, I believe there will be a peak centered more or less at 0km/s because the velocity of component is parallel or perpendicular to us so radial velocity is zero. The R at l=0,180 is not easy to obtain with this method. I am not sure about what are the other methods to obtain the particular R at those locations, but perhaps the other peaks around the 0km/s peak should give a good enough indication of the vicinity of the region that you are observing.

27 minutes ago, SteveBz said:

Exactly.  So maybe we can't share until all the calibrations have been done.

👍It would be a very exciting project when we can assemble data collected from different parts of the world to create a rather complete hydrogen intensity map of the sky! 

9 minutes ago, ZiHao said:

You can also convert the corrected rvs back to frequencies, but I don't recall seeing that in observations reports, maybe there are other reasons for doing that🧐

And like what you have mentioned in the previous posts, working on the raw frequency spectrums will be easier for baseline subtraction, noise etc. It should be the final spectrum that we are using to convert to LSR.

35 minutes ago, SteveBz said:

Both radial velocities and frequencies or just RVs?

There are corrections you can only apply if you haven't transformed to LSR.

You can also convert the corrected rvs back to frequencies, but I don't recall seeing that in observations reports, maybe there are other reasons for doing that🧐

6 hours ago, SteveBz said:

Thanks for this.  I've now updated it and it seems to work, but I've no idea why.  Here is the update:

Hmmm. What did you get when you enter radial_velocity=rv_corr1 instead for this?

observation = ICRS(ra=ra*u.deg, dec=dec*u.deg, \
                    pm_ra_cosdec=0*u.mas/u.yr, pm_dec=0*u.mas/u.yr, \
                    radial_velocity=0 *u.km/u.s, distance = 1*u.pc)

And rvs=rvs+rv_corr2. As we need to correct our barycentric radial velocity to the LSR, according to the post on stackexchange. I also need to test this procedure on other parts of the milky way, perhaps after my exams in June.

6 hours ago, SteveBz said:

I'm happy with the idea of .json, but I want to know if @ZiHao likes it too, or if he'd rather change it.

I will try to convert the np arrays to json type if there is sharing of data required, if that's you meant by changing😅

6 hours ago, SteveBz said:

But then we have a number of questions around radial velocity or frequency and should the velocity be corrected for LSR or not.  It seems to me if you put RV you have to correct it for LSR.  If you just have frequency, then not.

We always correct to LSR so that we can compare our data easily with everyone observing at different locations. Since we are working with velocities most of the time, I think that is probably why we usually have velocity as the x-axis.

6 hours ago, SteveBz said:

Hi Gents,

On a second point, my LSR correction does not seem to be working. Still.  At certain declinations, I'm about 10 km/s out, which means that slow-moving clouds can go from being moving towards to moving away (eg for 270 > l > 90 this does not make sense). 

Can I ask:

a) what declinations or values of galactic longitude are you able to access?  And what are your errors against the LAB survey?

b) Could you post your code here again please, so I can check mine against yours?

And here is mine:

def correct4LSR(self):
        # Define RV of 29.8 km/s
        rv = 29.8 * u.km / u.s 
        # Define V0 of 236 km/s
        V0 = 236 * u.km / u.s 
        # Define V0 of 220 km/s
        R0 = 8.15 * u.kpc
        # Define c0 of 299792.458 # km/s
        c0 = 299792.458 * u.km / u.s 
        for coords, record in self.raws.items(): 
            # Create Array of RAs
            ra=record['ra']
            #self.ras.append(ra)  # RA
            
            # Create Array of DECs
            dec=record['dec']
            #self.decs.append(dec) # DEC
            
            # Create Array of Peak signals
            sig_dB=max(record['sig_dBs'])
            #self.sig_dBs.append(sig_dB.max()) # peak dB value
            
            # Create Array of Speeds Lists
            speeds=record['rvs']
            #speedspeed0=speeds
            # Get time
            DT=record['DT']
            
            #Get time of exposure
            date_time_str=DT[:4]+'/'+DT[4:6]+'/'+DT[6:8]+' '+DT[8:10]+':'+DT[10:12]+':'+DT[12:14]
            dt=datetime.strptime(date_time_str,'%Y/%m/%d %H:%M:%S')
            t = Time(dt)

            loc = EarthLocation(lon=0*u.deg,lat=51*u.deg,height=88*u.m)

            #Correct for motion round the Sun
            sc = SkyCoord(record['ra']*u.deg,record['dec']*u.deg)
            vcorr = sc.radial_velocity_correction(kind='barycentric', obstime=t, location=loc)  
            vcorr= (float(str(vcorr)[:10])/1000)
            #print(vcorr)
            corrected_rv = speeds+vcorr

            #Correct for local group of stars
            observation = ICRS(ra=record['ra']*u.deg, dec=record['dec']*u.deg, \
                    pm_ra_cosdec=0*u.mas/u.yr, pm_dec=0*u.mas/u.yr, \
                    radial_velocity=0 *u.km/u.s, distance = 1*u.pc)

            new_rv = observation.transform_to(LSR()).radial_velocity 
            new_rv= (float(str(new_rv)[:10]))
            #print(new_rv)
            # Adjust Doppler velocity
            corrected_vel=float(str(new_rv)[:5])+corrected_rv
            self.raws[coords]['rvs']=corrected_vel
        return

Thanks.

Steve.

This line

observation = ICRS(ra=record['ra']*u.deg, dec=record['dec']*u.deg, \
                    pm_ra_cosdec=0*u.mas/u.yr, pm_dec=0*u.mas/u.yr, \
                    radial_velocity=0 *u.km/u.s, distance = 1*u.pc)

I think you should enter radial_velocity=vcorr. Try not changing the 'speeds' array first, when you get the final new_rv, only then you add the value to speeds. So this line has to be removed

 corrected_rv = speeds+vcorr

26 minutes ago, Victor Boesen said:

Thanks!! However, did you use my software to receive it? The DC spike is reaaaalllllyyyy bad 😅 Since it's so strong compared to the rest of the spectrum the std is not a good representation for the rest of the spectrum.

So it's not finding any of the other noise spikes.

Data is received with my python code, the raw data is first filtered by a median filter of size 30, or even higher IIRC, but I didn't use those in my attached data. I gave up on this clone as even though the sharp spikes like the DC spikes is removed, I am left with some broadband peaks which is very annoying to work with. Now with the genuine rtlsdr, I am using a median filter of size 15, and it got rid of the DC spike, leaving a clean spectrum. 

5 hours ago, Victor Boesen said:

I need your data 

The data is here! You might want to limit the y axis from -0.0005 to 0.0040, and the data is...you will know once you open it🥴, this is taken with the clone rtlsdr, no filters applied.

source20.npy velocity.npy

1 hour ago, Victor Boesen said:

I need your data 

I just wrote a simple function to identify noise spikes in the spectrum, and have tested it with some of my own data from yesterday. It works great, however, I'd like to test it on some other data than my own to confirm it works before I do any more tweaking  

If you could possibly send some debug files (preferably with the new version on GitHub) I'll do some testing on your data. Here are the first results with my own data:

In none of the cases were any bins from the hydrogen line itself replaced/modified, only the noise spike on the left!!

EDIT: Just updated the slant correction too. Now it's corrected for slant, after the noise spikes are removed. This yields the following result:

Victor

Very nice! I might send some old data taken few months ago later. Curious to know, did you locate approximately the frequency bins with the noise spikes and replace the spikes with the mean of nearby values?

7 hours ago, SteveBz said:

I'm still struggling with this.  

Line 1 has a typo 'rest_freq/rest_freq' .  I think it should be freq/rest_freq, no?

1) rv = raw radial velocity from observation, uncorrected.

2) vcorr = the speed of the Earth round the Sun, more or less.

3) rv = rv + rcorr ... Adjusted for speed of Earth round Sun

4) my_observation = plug rv into Astropy.

5) new_rv = observation adjusted for movement of local star group

6) corrected_vel=corrected_vel+new_rv  - what's happening here?  Isn't new_rv in 5) the answer?

Tx 

Steve.

 

I did rest_freq/rest_freq so that the raw uncorrected velocity is zero, and the rest of the code is just to focus on correcting the zero raw uncorrected velocity, and the final new_rv will just be a single value of the zero corrected velocity. And finally we just have to add this new_rv value to the corrected_vel array which originally was an array with raw uncorrected velocities, this shifts the spectrum. You can also correct for individual raw uncorrected velocity but I find that it takes long time to compute, therefore I use the zero point so that I can directly add to the raw array to shift all the values easily. Sry for the confusion!

On 18/04/2022 at 03:35, SteveBz said:

Perfect. That's what I got too.

What I don't have is the derivations for the other quadrants or outside the solar circle. Do you think they are different?

Alright, I have tried playing around with the quadrants and came up with these formulae:

Quadrant I: V_r=V_o*sinl((R_o/R_s )-1)  (basically just the same of what we have)
Quadrant II: same as Q1
Quadrant III: V_r=V_o*sinl((-R_o/R_s)+1)
Quadrant IV: same as Q3

We get these formulae if we always take the velocity component of source minus the velocity component of sun. The formulae for Q3 and Q4 differs by a negative sign, which means that the convention of redshift being positive and blueshift being negative is flipped now, you can still get the correct R_s but you would have to substitute a 'negative redshift into the equation' which is quite confusing. So just use the original formula and apply it to all quadrants, as it follows our convention of the signs.

If we were to take into account the arrows or the direction of the resolved velocity, and always take the velocity component of the object being 'chased' minus the velocity component of the object 'following', be it the sun or the source, then we will get the same formula for all quadrants.
 

7 minutes ago, Victor Boesen said:

Alright, thank you. I think I understand it now 😅 So, the barycentric_rv = rv + v_correction just corrects for the Earth's orbital velocity around the sun, whereas the LSR considers the Sun's motion relative to the observed area. In my case, developing my software, I am only interested in the barycentric radial velocity since the LSR velocity must depend on the distance to the area. For the majority of users (including myself) this is unknown for the most part. Thanks for sharing your resources on here! They are very helpful Do you have a GitHub account by any chance where I can have a look at some of your stuff?

Actually, the transformation to the LSR velocity is somehow independent of the distance to the area, as mentioned by the guy on stackexchange. I tried playing with the distance in the code just now and the final result did not change! 
I always wanted to gather all these information such as code on a website like blogspot, github but didn't do it because I have only two 'main pieces' of code: solar drift scan and hydrogen line, perhaps I should try to add all these in no matter the quantity. I will share it here once they are ready👍

11 minutes ago, Victor Boesen said:

@ZiHao Am I correct in assuming that the radial_velocity_correction function returns the Earth's velocity in the direction of the observed area?

In the sketch above the radial_velocity_correction would then return "-velocity correction"??

Victor

I think the general idea is to first correct earth's velocity to the solar system barycenter, and then we convert this to the local standard of rest, very vague explanation as I am struggling to visualize this also.
After computing the barycentric velocity correction, simply add this to your radial velocity, barycentric_rv = rv + v_correction. The text listed here https://astronomy.stackexchange.com/questions/41196/looking-for-a-routine-to-correct-for-local-standard-of-rest-lsr seems to be the only explanations I can find.

48 minutes ago, SteveBz said:

Actually, it seems to be nearly linear in Dec.  

Interesting pattern, might be the pointing error, or the code itself, not too sure... I have only one set of data for now so can't really test it reliably. We need other sets of data🧐

1 hour ago, SteveBz said:

OK, I just used your code.  It improved my code by about 7 km/s in the right direction, but I think I'm still about 10 km/sec too fast.  See here (the blue is corrected and the orange uncorrected).

Tx

Steve.

Could it be pointing error in the mount? Maybe you could try pointing in the zenith, I think it's much easier.

1 hour ago, SteveBz said:

I wonder @Victor Boesen, should I have calibrated my rtl-sdr?  The guy in the video suggested by @ZiHao talked about calibration.  What do you think? ZiHo, did you calibrate yours?

Steve

Nope, I didn't calibrate mine, like what Victor had mentioned, 1-2ppm shouldn't cause that much of a difference.

46 minutes ago, Victor Boesen said:

I have two RTL SDR's and both of them are accurate within 2PPM (probably even lower). However, you can definitely try to calibrate in on a known signal anyways.

Btw Victor, for the pyrtlsdr wrapper, in sdr.freq_correction, it seems that a value has to be entered, how did you deal with this?