[Top][All Lists]

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

Re: Problems loading a program in octave

From: Joe Koski
Subject: Re: Problems loading a program in octave
Date: Thu, 15 Jun 2006 16:14:49 -0600
User-agent: Microsoft-Entourage/

on 6/15/06 3:51 PM, Hayden Rampadarath at address@hidden wrote:

I am having problems loading a program.
Everytime I enter 'load C:\plates\read_hayden.m'
I get
'error: load: C:\plates\read_hayden.m: inconsistent number of columns near line 2'
'error: load: unable to extract matrix size from file 'C:\plates\read_hayden.m'
any ideas why???

I had the same problem this morning on my Mac. It usually indicates that there is something wrong with the data file format. With an editor, check to make sure that you have the correct (DOS or UNIX) line endings. This morning, I had to delete some hidden characters at the end of the file, and then resave my data with a new blank line at the end.


Please help!!!!!!
p.s. I have attached the program
 Send instant messages to your online friends

# short program to process Hayden's files.

# Hi Hayden - I have modified the format of your data files to make it easier
# to read in Octave - a file of numbers is easy to read with the 'load'
# command

# The data are processed in the same way as I did for CF Oct

# I had lost the original program, but I did this one this morning
# I hope it works.  Let me know what you get with the real data.
# I've run it with 'made up data'

# Good luck.  John.

# firstly read in the B and B-V data for the field stars
# using Octave's load command - a quick way to get a file of numbers
# into Octave

load /cygdrive/c/plates/B_and_B_V_list.txt

# I am using Octave via cygwin, so the file path looks like a linux path.
# In windows it would just be:
#     c:\john_inn\astronomy\hayden\B_and_B_V_list.txt

# this makes a matrix called B_and_B_V_list

# pull out the B data for the field stars by taking the 1st column
# and the B-V data from the 2nd column

B_field = B_and_B_V_list(:,1);
B_V_field = B_and_B_V_list(:,2);

plot (B_field, B_V_field, "33");

# B-V of the target object, for later use

OJ287_B_V = 0.45;

# load the modified data file

load /cygdrive/c/plates/modified_plate_format.txt

size_mod = size(modified_plate_format);  # finds how much data we read

n_plates = size_mod(1); # this is the number of plates in the file

# pull out the data by column

plate_number = modified_plate_format(:,1);
datep = modified_plate_format(:,2);

# f1 - f12 are field star counts
# OJ is the target counts
# I_sky is the sky counts

f1 = modified_plate_format(:,3);
f2 = modified_plate_format(:,4);
f3 = modified_plate_format(:,5);
f4 = modified_plate_format(:,6);
f5 = modified_plate_format(:,7);
f6 = modified_plate_format(:,8);
f7 = modified_plate_format(:,9);
f8 = modified_plate_format(:,10);
f9 = modified_plate_format(:,11);
f10 = modified_plate_format(:,12);
f11 = modified_plate_format(:,13);
f12 = modified_plate_format(:,14);
OJ = modified_plate_format(:,15);
I_sky = modified_plate_format(:,16);

# clear the graphics window in case there is anything there already

clg ;

plot (datep-1968., f1, "01")  ; # 01 means black crosses in my octave

# 'hold on' means the next plot is plotted on top of previous plots

hold on

plot (datep-1968., OJ, "04") ; # 04 means black squares in my octave

# set up the labels for the axis

xlabel ("Date - 1968.0")
ylabel ("Raw counts from Maxim DL")

# replot redraws the plot after something has changed - e.g. the labels in
# this case


fprintf(stderr,"Pausing so you can look at this plot - hit return to continue\n ");
pause () # Pause

# turn off overplotting

hold off

# Hayden: Now I assume that I have to subtract I_sky from all 'counts'
# and negate the result?  I am not sure, only half remembering what you
# told me in an earlier email

zero_pt = 27.5  ; # this is to get the magnitudes close to the B mags

# Here are the raw magnitudes....

mag_f1 = -2.5*log(I_sky - f1) + zero_pt;
mag_f2 = -2.5*log(I_sky - f2) + zero_pt;
mag_f3 = -2.5*log(I_sky - f3) + zero_pt;
mag_f4 = -2.5*log(I_sky - f4) + zero_pt;
mag_f5 = -2.5*log(I_sky - f5) + zero_pt;
mag_f6 = -2.5*log(I_sky - f6) + zero_pt;
mag_f7 = -2.5*log(I_sky - f7) + zero_pt;
mag_f8 = -2.5*log(I_sky - f8) + zero_pt;
mag_f9 = -2.5*log(I_sky - f9) + zero_pt;
mag_f10 = -2.5*log(I_sky - f10) + zero_pt;
mag_f11 = -2.5*log(I_sky - f11) + zero_pt;
mag_f12 = -2.5*log(I_sky - f12) + zero_pt;
mag_OJ = -2.5*log(I_sky - OJ) + zero_pt;

ylabel ("Raw Differential Magnitude OJ - field star 1");
plot (date, mag_OJ - mag_f1, "04")

fprintf(stderr,"Pausing so you can look at this new plot -hit return to continue\n ");
pause (); # Pause

# Now we want to do the SVD business...

# For each plate (i.e. for each line in the modified_plate_format file)
# we have a set of 12 magnitudes, M.
# We know the true B and B-V of these 12 stars.
# we want to solve for the betas
#  B = beta1*M + beta2*M^2 + beta3(B-V) + beta4

beta=zeros(n_plates,4) ; # makes an array to store the beta values
#                           4 betas per plate

# make an array to store the transformed field star magnitudes

corr_mags_plate=zeros(n_plates, 12);

# make an array to hold the corrected OJ magnitudes

# We can now do each plate, one at a time, in a 'for' loop

for jj=1:n_plates

dd = linspace(1,12,12)*0 + 1.0  ; # i.e. first guess for beta4 is 1

# pull out the mags of the field stars for a given plate

pla_obs = [ mag_f1(jj),mag_f2(jj), mag_f3(jj), mag_f4(jj), \
            mag_f5(jj),mag_f6(jj), mag_f7(jj), mag_f8(jj), \
            mag_f9(jj),mag_f10(jj), mag_f11(jj), mag_f12(jj)];

# form the matrix we want to solve
# i.e. soln dependent on M, M^2, B-V and a zero point term ('dd')

matrix_a = [pla_obs; pla_obs.*pla_obs;B_V_field'; dd];

# do the singular value decomposition of matrix_a to look at the singular
# values


# the singular values are s
# From these data, the s take on 4 values for each plate,
# these are roughly 800, 3, 1, and 0.05
# The 0.05 values are the noise values, and we want to discard solutions
# that use these.  (See the MNRAS paper on CF Oct.)

# find the psuedo inverse of matrix_a, uses tolerance of 0.1 mag,
# which is the approx obs error and which will cut off the smallest singular
# values - i.e those near 0.05

[soln]= pinv(matrix_a, 0.1)  ; #

# calculate the 4 beta values for this plate

 beta(jj,:) =   vec((B_field'*soln))' ;
 # we can get 'corrected' magnitudes for the field stars by taking the pla_obs
 # array and transforming it using the betas we just found...
 corr_mags_plate(jj,:) = (matrix_a'*beta(jj,:)')' ;
 # now do a plot comparing the transformed mags and the B mags
 axis ([15, 11, 15, 11]);  # set up the limits to the plot
 xlabel("B mag")
 ylabel("Transformed Plate mag")
 plot (B_field, corr_mags_plate(jj,:) , "04");
 fprintf(stderr,"Press return to continue \n")
 # Now we need to apply the betas for this plate to the OJ obs
 Corr_OJ(jj) = beta(jj,:)*[mag_OJ(jj); mag_OJ(jj)*mag_OJ(jj); OJ287_B_V; 1.0];
 end;  # end of the loop
 # enter multiplot mode to plot the betas for the different plates
 clg;      # clears the graphics screen
 hold off; # make sure we are not in overplot mode
 axis;     # remove the limits we put onto the plot axes
 automatic_replot = 0; # turns off automatic redraw in my octave
 #                       if there's a problem try setting to 1
 xlabel ("Date - 1968.0");
 subplot (4,1,1)
 ylabel("Beta 1 = M term")
 plot (datep-1968., beta(:,1), "04")
 subplot (4,1,2)
 ylabel("Beta 2 = M^2")
 plot (datep-1968, beta(:,2), "04")
 subplot (4,1,3)
 ylabel("Beta 3 = B-V")
 plot (datep-1968, beta(:,3), "04")
 subplot (4,1,4)
 ylabel("Beta 4 = zero pt")
 plot (datep-1968, beta(:,4), "04")
 fprintf(stderr,"Press return to continue \n")
oneplot();  # back to one plot per page


# now to plot the data - finally!
ylabel("Corrected OJ287 B mag")

title("OJ287 (squares) and field star 10 (+)")

plot (datep-1968.0, Corr_OJ, "04");

hold on

# plot field star 10 corrected mags for comparison
plot (datep-1968.0, corr_mags_plate(:, 10), "01")

hold off

# end of program

Help-octave mailing list

reply via email to

[Prev in Thread] Current Thread [Next in Thread]