help-octave
[Top][All Lists]
Advanced

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

Fwd: Problems loading a program in octave


From: Hayden Rampadarath
Subject: Fwd: Problems loading a program in octave
Date: Thu, 15 Jun 2006 23:16:03 +0100 (BST)



Note: forwarded message attached.

Send instant messages to your online friends http://uk.messenger.yahoo.com

--- Begin Message --- Subject: Problems loading a program in octave Date: Thu, 15 Jun 2006 22:51:40 +0100 (BST)
Hi,
 
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???
 
Please help!!!!!!
 
Hayden
 
p.s. I have attached the program

Send instant messages to your online friends http://uk.messenger.yahoo.com

# 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

replot

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

[u,s,v]=svd(matrix_a,0);

# 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.)
#
s;

# 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")
 pause()
 
 # 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")
 pause()
 
oneplot();  # back to one plot per page

ylabel("")


# 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


--- End Message ---

reply via email to

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