|
From: | Joe Koski |
Subject: | Re: Problems loading a program in octave |
Date: | Thu, 15 Jun 2006 16:14:49 -0600 |
User-agent: | Microsoft-Entourage/11.2.3.060209 |
Hi,Hayden,
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
_______________________________________________
Help-octave mailing list
address@hidden
https://www.cae.wisc.edu/mailman/listinfo/help-octave
[Prev in Thread] | Current Thread | [Next in Thread] |