SteveBz

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?

18 hours ago, Victor Boesen said:

Yes, that's the mistake. It basically means the "raw" uncorrected radial velocity is used as the radial velocity as it would be from the barycenter.

Which, like shown in the beautiful drawing above, is not the case!

Yes, that's right.  The only issue I'm having is that the calculations seem to return numbers with different sign conventions, or so it seems.  To get the right answer you seem to have to (and I may have got this wrong), but you seem to have to add the one and subtract the other 🥴

9 hours ago, Victor Boesen said:

I'm handling the observed radial velocity as the radial velocity with respect to the barycentric frame of reference (ICRS).

I'm not sure what that means? Is that the mistake? You mean you haven't added the earth round the sun?

57 minutes ago, ZiHao said:

👍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! 

It would be great.

4 minutes ago, ZiHao said:

It should be the final spectrum that we are using to convert to LSR.

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

10 hours ago, ZiHao said:

We always correct to LSR

Both radial velocities and frequencies or just RVs?

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

10 hours ago, ZiHao said:

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

I was suggesting that we all try to use the same file type and the same format to share data.

But actually,  I'm now thinking that I'll just write a conversion routine to switch between formats. It's unlikely that the LAB site at Bonn will want to change 😂.

Kind regards,

Steve.

I just downloaded a LAB file. The format is very old fashioned and I don't recommend it, however it's interesting to note the info they included in the file:

v_lsr [km/s], T_B [K], freq. [Mhz], wavel. [cm]
    -399.83      0.13   1422.300   21.07800
    -398.80      0.01   1422.295   21.07808

Steve

19 minutes ago, Victor Boesen said:

All were received with the same parameters. Besides the default parameters the following parameters were used:

50000 FFTs
10 Median

 

Some of your pngs have better noise than others? Why is that?

14 minutes ago, Victor Boesen said:

I'm assuming that this because I haven't calculated the LSR for each observation and only corrected for the Earth's velocity. Will do some more coding to try and calculate the LSR!!

The LSR code that I posted above seems to work. The corrections for the area I was looking at had LSR = 5 and planetary movement - 5 - or the other way round, but they were both about the same magnetude, but of different sign. To get it to match I had to reverse one of the signs, then it worked beautifully.

Steve.

12 minutes ago, Victor Boesen said:

 

Wow, really cool. How did you track the galactic plane?

1 hour ago, Victor Boesen said:

Great idea! I am still pondering about improving the data included in the debug file.

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 also think it's not a debug file, because it's everything.  It's really the output data of the whole thing, so a new name would be a start

I agree with your idea if putting all the run parameters in it.

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.

Maybe it could be a runtime option (raw frequency or corrected RV).  Frequency is easier to filter for electronic noise because the spikes are always in the same place. In fact I make some artificial 'flats' from the scans in the areas of empty space and then apply them to the other areas.  It would probably work with your big spike too.

In the longer term, I'll probably want to use a big-data approach with Pandas Dataframes for the work I do.  It would be much faster, but a learning curve for any new users.

I'd like to hear other views.

Kind regards

Steve.

1 hour ago, ZiHao said:

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

 

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

  
    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 Speeds Lists
            rvs=record['rvs']
            rvs0=rvs
                        
            #Get date, time and location of exposure
            DT=record['DT']
            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)

            # Create Array of RAs
            ra=record['ra']
            # Create Array of DECs
            dec=record['dec']
                        
            #Correct for motion round the Sun
            sc = SkyCoord(ra*u.deg,dec*u.deg)
            rv_corr1 = sc.radial_velocity_correction(kind='barycentric', obstime=t, location=loc)  
            rv_corr1= (float(str(rv_corr1)[:10])/1000)
            #print(vcorr)
            #Correct for local group of stars
            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)

            rv_corr2 = observation.transform_to(LSR()).radial_velocity 
            rv_corr2= (float(str(rv_corr2)[:10]))
            # Adjust Doppler velocity
            rvs = rvs+rv_corr1-rv_corr2
            self.raws[coords]['rvs']=rvs
            
        return

The key line is

 rvs = rvs+rv_corr1-rv_corr2

Strange, no?  Also, once I plotted the final spiral plot, I had to do:

X=-X

because the Galaxy was rotating in the wrong direction and also

rvs=-rvs.

to change the blue/red shift colours. So I really need to chase these through and understand them properly and make sure they are legitimate.

Kind regards

Steve.

So here it is.  Still got a few funnies.  Most of the numbers were back to front and I need to trace it through and there's a strange little spur South of Sol, but otherwise nearly there.  It covers 2 quadrants (IV & III, and actually I can see quadrant II as well). So there's still a bit more work to be done, but I have been surprised at how well it's come out:

You can easily see the Orion spur leading into the Perseus arm and the outer arm.  There is also that blue arm, that appears on the map without a name.

Kind regards

Steve.

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.

Hi,

@Victor Boesen and @ZiHao, After we collect our data we store it and analyse it.  At the moment I think we're all using different and incompatible data storage mechanisms and formats, so that when, for instance, Victor wanted to borrow ZiHao's data it was in a different format.  My data is in Victor's old format, but modified, so it was also awkward.  Can I suggest we all agree to a common standard and format?  It wouldn't be hard for any of us to change I think, but it would make our data more shareable. At the moment ZiHao appears to be using Python's NumPy format Victor is using a newer json format and I am using an older, but modified, json format.   

There are many formats out there.  I don't know NumPy format, but it might be useful.  I'm comfortable with json, but there are even other Pandas-style formats like pickle, which might be useful.

Can we have a discussion about what container (ie Numpy, json, csv, pickle etc to use) and then, separately, what format to use within that container?

While I'm happy with json, I'd also be happy to consider an alternative.

What do you think?

Kind regards,

Steve.

2 minutes ago, Victor Boesen said:

Excellent, thanks a lot!! This will keep me busy for the next year or so 😅 Really helpful!

There's some .pngs there that you don't really need.  You can delete them.  I've altered the 'parameters' key to hold:

 

{
    "Parameters": {
        "sample_rate": 2400000,
        "ppm": 0,
        "resolution": 11,
        "num_FFT": 10000,
        "num_med": 75,
        "ra": 57.6,
        "dec": 61.0,
        "l": "142d58m28.6283s",
        "b": "5d17m52.6504s",
        "const": "Camelopardalis",
        "DT": "20220409150351"
    },

In particular, the DT (Datetime) allows you to use @ZiHao's trick to convert to the LSR of the sample date and time.  My lat and long are about 51 N, 0 E, 88 m.

Here it is zipped up.  Not sure about wetransfer, I don't have a subcription.

Sorry

Steve.

ra16-61v1.tar.gz

8 minutes ago, Victor Boesen said:

Should be perfect!! Perhaps it would be the easiest to upload it to wetransfer?

Can you PM me your email adress then?

Steve.

On 20/04/2022 at 08:32, Victor Boesen said:

I need your data 

Victor, just seen this.  You can have as much as you like of mine, but it's the old format. It that OK?

Kind regards

Steve.

So here's a small step.  Broadly the arcs are 'spiral-like'.  Broadly, too, they centre at the Galactic Centre.  I think I need to step through the calculations again, because the distances don't look right and the range of velocities doesn't seem large enough.  The signs of the velocities, however, look about right.  I've used automatic Gaussian fitting with a fixed Sigma, which is probably not realistic.  I hoped it might work because I hoped the thermodynamic temperature range of the bands would be similar.

There's some funny stuff in the X is negative range.

Comments welcome.

Steve.

On 17/04/2022 at 16:39, ZiHao said:

rv=c*((rest_freq/rest_freq)-1)*u.km/u.s

t = Time(, format='mjd')

loc = EarthLocation(lon=,lat=,height=)
ra=*u.deg
dec=*u.deg
sc = SkyCoord(ra,dec)

vcorr = sc.radial_velocity_correction(kind='barycentric', obstime=t, location=loc)  
rv = rv + vcorr

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

new_rv = my_observation.transform_to(LSR()).radial_velocity

corrected_vel=corrected_vel+new_rv

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.

9 hours ago, Victor Boesen said:

Thank you Steve!! Still got some work to do with the noise spikes removal. No need to apologize for not testing it yet!

I do my observing besides multiple apartments in the city so far from ideal. Parameters used are the following:

"DSP Parameters": {
        "number_of_fft": 50000,
        "resolution": 11,
        "median": 10
}

Victor

Have you tried a higher median,  like 50? It won't get ride of the spike, but I think it may smooth the noise? If your curve is smoother it's easier to do Gaussian fitting.

If you can reliably do max(freqs) you can get the mean of the highest peak and then subtract a fake gaussian (norm.pdr) from the curve. I'll post some code later.

20 minutes ago, Victor Boesen said:

I managed to do some testing of my new rewritten software. After seeing @ZiHao using Astropy I decided to try and switch to astropy from pyephem, which made things easier when correcting for radial velocity. Here are three observations from today and a comparison to the same area from the H-line survey.

And from H-line survey. Note that it's just mirrored since they plot it as a function of increasing radial velocity, whereas I plot mine as a function of increasing frequency.

Next target:

And final target which greatly shows the correlation between my spectrum and theirs. Look at the peak which is basically centered at 0 in both cases.

Annoyingly, the noise spike is just inside the filter I apply when searching for the hydrogen line peak, so next thing on my todo-list will be removing these large spikes! This will, however, be a lot easier to work on now that I have some debug files I can test on

I was also able to test out the new UI, which I really liked, and everything worked as it should. Therefor, I have merged the rewrite branch with the main branch on GitHub
https://github.com/byggemandboesen/H-line-software

Victor

Really nice, Victor.  Very good correlation.  Once I've stabilised my code again, I'll refresh my copies of H-line and give you some feedback. I'm sorry I've not done it so far.

Where are you based in terms of other houses and stuff?  Are you in the country or town?  How many FFTs are you doing and what about the Median smoothing parameter?

The stars are out tonight.  It's quite hard to align your dish!  At my estimate the dish is about 1.5 degrees (my thumb's width)  to the East of Polaris and Polaris itself is about .7 degrees West of NCP. Therefore the dish is within 1 degree of NCP.  It's good enough. So I don't really understand why I have this structural error.

On top of all of this, the code is becoming a nightmare and I really need to refactor it 😅

Steve.