function srtdata=readsrt(filename); % function srtdata=readsrt(filename); % version 1.5 18/01/11 % version 1.4 28/02/07 % version 1.3 06/02/07 % version 1.2 30/01/07 % Usage % srtdata=readsrt() % % This function returns a structure containing the data in . It contains the fields: % % time: [n x 1 double] => UTC time in matlab serial date number format (help datenum) % timestring: [n x 20 char] => UTC time in easy to read string form '19-Oct-2005 23:11:40' % timeelapsed: [n x 1 double] => time elapsed in seconds since the first observation % lst: [n x 1 double] => local sidereal time (day fraction) in matlab serial date number format % ra: [n x 1 double] => right ascension of observation in hours % dec: [n x 1 double] => declination in degrees % az: [n x 1 double] => azumith % el: [n x 1 double] => elevation % az_off: [n x 1 double] => azumith offset % el_off: [n x 1 double] => elevation offset % freq: [n x 1 double] => frequency % freq_bw: [n x 1 double] => frequency bandwidth % mode: [n x 1 double] => mode % nchannels: [n x 1 double] => number of channels % vlsr: [n x 1 double] or []=> vlsr data, if present. If not empty array returned % data: [n x nchannels double] => all channel data % datablock: [n x 1 double] => block number - see below % blockheader: {1 x nblocks cell} => block header - see below % (n=number of observations) % % Notes % ------ % To access a field use eg: % t=srtdata.time; %extract a vector containg the time of each observation % imagesc(srtdata.data); % % The field contains which observing block the data belongs to. If you are % running an srt script each new command generates a new block. The telescope command that % caused the new block is recorded in the field which is a CELL array. You % need to use curly brackets to get its contents: % eg: % disp(srtdata.blockheader{1}) % * STATION LAT= 56.00 DEG LONGW= 4.00 % * srt.cmd: line 15 : record moonman2.rad % % To select data just from a specific block use the FIND command to get an index of all the % observations with this block number. % For example to get just block 6 data: % % I=find(srtdata.datablock==6); % I is a vector containing the row indicies % data6=srtdata.data(I,:); % choose rows I and all colums form the data set % % Advanced example: % Plot the time average of all the channels in block 4: % % plot(mean(srtdata.data(find(srtdata.datablock==4),:),1)) % (hint - the final '1' means take a mean in the fist (row) dimension) % % See also DATENUM, DATEVEC, DATETICK % % Any problems - email hugh@astro.gla.ac.uk % 30/11/05 - added reading of vlsr data % 18/01/11 - corrected (minor) error in JD calculation % -------------------------------------------------------------------------------------- if not(exist(filename,'file')), error(['File not found: ' filename]); end; fid=fopen(filename); blocknumber=0; ndatalines=0; blocknumberlist=[]; nfileline=0; previouscommentline=NaN; while not(feof(fid)), % read the file line by line - do this all together for speed fileline = fgetl(fid); nfileline = nfileline+1; %does the line start with a *? if so its a comment if fileline(1)=='*' %More than one consecutive comment - add to the previous one if (previouscommentline+1) == nfileline, comments{blocknumber}=strvcat(comments{blocknumber},fileline); else %otherwise its a new data block blocknumber=blocknumber+1; comments{blocknumber}=fileline; previouscommentline=nfileline; end; else % if not save to a cell array, good cos cell arrays dont need contiguous memory space ndatalines=ndatalines+1; blocknumberlist(ndatalines)=blocknumber; %for looking up the block comment later datalines{ndatalines}=fileline; end; end; fclose(fid); %find out how many columns the original data has - assume it is constant through the file % if the velocity is in the file the penultimate column has 'vlsr in it dtest=(textscan(datalines{1},'%n','delimiter',':| |(vlsr)','multipleDelimsAsOne',1)); numcols=numel(dtest{1}); if numel(findstr(datalines{1},'vlsr'))==0, lastdatacol=numcols; % number of columns in first line of data vlsr=false; else lastdatacol=numcols-1; % the last col is the vlsr data in this case vlsr=true; end; rawdata=zeros(ndatalines,numcols); %reserve space for data arrray - for speed % now loop through the string data converting it to numbers for i=1:numel(datalines), dataline=textscan(datalines{i},'%n','delimiter',':| |(vlsr)','multipleDelimsAsOne',1); if numel(dataline{1})~=numcols, error('readsrt.m:columnerror','The number of columns in the data changed. You''re on your own here I''m afraid'); else rawdata(i,:)=dataline{1}; end; end; [nr,nc]=size(rawdata); %make a structure containing the data % first make the time: year=rawdata(:,1); dayofyear=rawdata(:,2); h=rawdata(:,3); m=rawdata(:,4); s=rawdata(:,5); srtdata.time=(datenum(year,1,0)+dayofyear+datenum(0,0,0,h,m,s)); srtdata.timestring=datestr(srtdata.time); srtdata.timeelapsed=(srtdata.time-srtdata.time(1))*60*60*24; srtdata.az=rawdata(:,6); srtdata.el=rawdata(:,7); srtdata.az_off=rawdata(:,8); srtdata.el_off=rawdata(:,9); [srtdata.lst,srtdata.ra,srtdata.dec]=readsrt_utc2sidereal(srtdata.time,srtdata.el+srtdata.el_off,srtdata.az+srtdata.az_off); srtdata.freq=rawdata(:,10); srtdata.freq_bw=rawdata(:,11); srtdata.mode=rawdata(:,12); srtdata.nchannels=rawdata(:,13); if vlsr, srtdata.vlsr=rawdata(:,numcols); else srtdata.vlsr=[]; end; srtdata.data=rawdata(:,14:lastdatacol); srtdata.datablock=blocknumberlist'; srtdata.blockheader=comments; % ----------------------------------------------------------------------------------------------------------- function [LST,RA,dec]=readsrt_utc2sidereal(UTC,alt,az); % UTC to sideal time conversion, correct for Glasgow university observatory (55?54'8.646" N) % UTC is an array of matlab datenumbers - create using DATENUM % alt (altitude in degrees) and az (azimuth in degrees (east from north)) are the pointing of the % telescope - get from the header of the .rad file % LST good to about a second until 2099! % eg. % [lst,ra,dec]=utc2sidereal(datenum(2001,12,19,18,0,0),10.0607,173.7908); % [lst,ra,dec]=utc2sidereal(datenum('19.05.2000 12:00:00','dd.mm.yyyy HH:MM:SS'),10.0607,173.7908) % % also az/alt/ST to ra/dec (apparent, i.e., unprecessed) % See also DATENUM, DATEVEC obs_long=-(4+18/60 + 24.928/3600); % longitude of observatory is -4?18'24.928" E obs_lat =(55+54/60+8.646/3600); % latitude of observatory is 55?54'8.646" N [YY,MM,DD,h,m,s]=datevec(UTC); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % now the calculation... % calculate the Julian Day number, JD (note this starts at noon) % ## this version has an error due to fix issues... % j2 = fix((j-14)/12); % JD = DD - 32075 + fix(1461*(YY+4800+j2)/4) ... % + fix(367*(MM-2-j2*12)/12) ... % - fix(3*fix((YY+4900+j2)/100)/4); % corrected 18/01/2011 JD_2000= 2451545; % julian day num at 1200 (midday), 1/1/2000 MD_2000= 730486; % matlab day num at 1200, 1/1/2000 MD=floor(datenum(UTC-12/24)); %matlab day num with offset changing at midday JD=MD-MD_2000+JD_2000; % calculate Julian centuries, JC, between 12h 1/1/2000 and 0h on % the calendar day in question JC = (JD-2451545.5)/36525; % calculate the Greenwich mean sidereal time at 0h, in seconds GMST = 6*3600 + 41*60 + 50.54841 ... + 8640184.812866*JC ... + 0.093104*JC.^2 ... - 0.0000062*JC.^3; % add on the seconds in the rest of the day GMST = GMST + (h*3600 + m*60 + s)/0.9972696; % convert to fraction of a day GMST = GMST/(3600*24) - floor(GMST/(3600*24)); % calculate the local sidereal time in day fraction LST = GMST + obs_long/360; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %now calculate the hour angle, HH, and dec, in degrees, from the alt and az: az = az/180*pi; alt = alt/180*pi; obs_lat = obs_lat/180*pi; % convert to radians dec = 180/pi*asin(sin(alt).*sin(obs_lat)+cos(alt).*cos(obs_lat).*cos(az)); HH = 180/pi*(atan2(-cos(alt).*sin(az),(sin(alt)*cos(obs_lat)-sin(obs_lat)*cos(alt).*cos(az)))); %finally calculate right ascension, RA, in hours, within principal range RA = (LST - HH/360)*24; RA = mod(RA,24);