Ho trovato la causa dell'errore.
Contrariamente a quanto ho scritto nelle risposte precedenti (e mi scuso per il mio errore) il problema risiede solo indirettamente nel numero di righe della matrice "
p" e, principalmente nell'algoritmo:
[*] al termine di ogni iterazione del loop interno, il vettore "
Hv viene assegnato al cellarray "
Hvf"
Hvf{i}=Hv;
[*] successivamente, tramite la funzione "
cellfun" viene calcolato il valore
massimo del cellarray
[max_Hvf,idx_max]=cellfun(@max,Hvf);
[*] Il valore dell'
indice-posizione del valore massimo viene usato per accedere al corrispondente elemento della prima colonna della matrice "
p"
Ti_max=p(idx_max,1);
[*] Al termine della prima iterazione del loop esterno, non ci sono problemi
Il problema sorge a partire dalla seconda iterazione.
[*] al termine della seconda iterazione, il
cellarray Hvf contiene due vettori (i due vettori
Hv calcolati nel loop interno
[*] l'istruzione
[max_Hvf,idx_max]=cellfun(@max,Hvf);
ritorna in
idx_max due valori: gli indici-posizione dei valori massimi dei due vettori contenuti nel cellarray "Hvf"
[*] quando nell'istruzione suvvessiva, "
idx_max" viene usato per accedere elementi della prima colonna della matrice "
p" si genera l'errore.
Supponendo che "idx_max" contenga i valori (15444, 4748), l'istruzione
Ti_max=p(idx_max,1);
si traduce in
Ti_max=p([15444, 4748],1);
il che significa "prendere i valori presenti nelle righe "15444" e "4748" della prima colonna della matrice "
p" ed assegnarli a Ti_max.
Il primo valore (15444) identifica la posizione del massimo del primo elemento del
cellarray.
Siccome ad ogni iterazione del loop principale la matrice "
p" viene aggiornata, se nel corso della seconda iterazione la matrice "
p" non ha un numero di righe uguale (o maggiore) a quello della prima iterazione, l'istruzione "cerca" di accedere alla riga "15444" della matrice "
p" che "non esiste", perchè la matrice "
p", nella seconda iterazione, ha meno di "15444" righe.
Questo è il motivo per il quale viene viene generato il messaggi di errore.
La soluzione dipende dall'algoritmo che vuoi implementare.
Quando, al termine del loop interno cerchi il valore massimo all'interno del cellarray "Hvf" devi cercarlo tra tutti gli elementi del cellarray calcolati fino a quel moneto o solo nell'ultimo?
[*] nel primo caso avrai tanti "valori massimi" quanti sono i cicli effettuati fino al quel momento dal loop esterno.
Se ti serve un solo valore (per accedere agli elementi della prima colonna della matrice "p")
devi gestire questi "n" valori opportunamente considerando il numero delle righe della matrice "p" nel corso della iterazione attuale.
[*] nel secondo caso, a meno che il cellarray con tutti i vettori "Hv" non ti serva in un altra parte del codice che non hai pubblicato puoi semplicemente "non usare il cellarray, ma identificare il massimo direttamente nel vettore "
Hv" per cui le due istruzioni:
Hvf{i}=Hv;
[max_Hvf,idx_max]=cellfun(@max,Hvf);
possono diventare:
% Hvf{i}=Hv;
[max_Hv,idx_max]=max(Hvf);
Spero di essere riuscito a spiegarmi, la spiegazione è un po' complicata da scrivere.
In sostanza devi rivedere la seconda parte dell'algoritmo e, in particolare, come usare i valori memorizzarti nel cellarray "Hvf".
In questo, purtroppo, non posso aiutarti.
Posso, invece, indicarti alcuni miglioramenti che potresti apportare al codice per renderlo un po' più veloce:
Funzione "Parametri"
Nella funzione "Parametri.m" calcoli i valori dei parametri "
s1", "
s2", "
s3" ecc.
Dal momento che sono delle costanti, non è necessario calcolarli "tutte" le volte che invochi la funzione "Parametri" (nel caso dell file di input "frequenza05mhz.mat" li calcoli 4 milioni e mezzo di volte).
Puoi spostare quella sezione fuori dalla funzione, all'inizio dello script principale e passarli come parametri alla funzione.
Per "limitare" il numero di parametri da passare, potresti anche definirli cole campi di una struttura e passare alla funzione la struttura.
Nella funzione c'è l'"
if"
if Tc==Th %DeltaT=0
Tc e
Th sono
numeri reali, considerando le approssimazioni floating point, la possibilità che due numeri reali siano "
esattamente uguali" (condizione che testi nell'"
if" è estremamente bassa (per non dire "quasi impossibile") per cui le possibilità che il codice esegua le istruzioni contenute nell'"
if"
SM= s1 + s2*Tc + s3*Tc^2 + s4*Tc^3;
RM= r1 + r2*Tc + r3*Tc^2 + r4*Tc^3;
KM= k1 + k2*Tc + k3*Tc^2 + k4*Tc^3;
sono di conseguenza estremamente basse.
Per curiosità, prova a a fare questo test con Matlab e vedi ... cosa succede:
a=0.1 + 0.1 + 0.1
b=0.3
if(a == b)
disp('"a" è uguale a "b"')
else
disp('"a" è diverso da "b"')
end
a-b
Il consiglio, nel caso tu debba confrontare numeri reali, è utilizzare un valore di soglia (opportunamente scelto) come elemento discriminante.
L'"
if" in questione potrebbe essere riscritto così:
soglia=1.0e-05
if(abs(Tc-Th) <= soglia)
SM= s1 + s2*Tc + s3*Tc^2 + s4*Tc^3;
RM= r1 + r2*Tc + r3*Tc^2 + r4*Tc^3;
KM= k1 + k2*Tc + k3*Tc^2 + k4*Tc^3;
else
SMth= s1*Th + (s2*Th^2)/2 + (s3*Th^3)/3 + (s4*Th^4)/4;
SMtc= s1*Tc + (s2*Tc^2)/2 + (s3*Tc^3)/3 + (s4*Tc^4)/4;
[ ... ]
end
Per identificare il valore della soglia (vedi sopra) puoi inserire, prima dell'inizio del loop interno:
tmp_Tc=p(:,2);
tmp_Th=p(:,3);
min_diff=min(abs(tmp_Tc-tmp_Th))
soglia=0.0001; % questa istruzione dovrebbe essere spostata fuori dal loop
if(min_diff <= soglia)
disp(['Input file: ' files,freq{i},misure,'.mat'])
disp(['min_diff = ' num2str(min_diff)])
disp([' Possibile Tc == Th non catturato in "Parametri.m"'])
end
continue
e far girare lo script principale; l'istruzione continue, evita di eseguire il loop interno.
Una volta identificata la soglia, puoi rimuovere questa porzione di codice (e, soprattutto l'istruzione "
continue")
script principale
Puoi inserire, prima dei due loop, il calcolo dei parametri "
s1", "
s2", "
s3" ecc.
%
% Definizione parametri delle equazioni (from "Parametri.m")
%
param.s1=1.33450*10^(-2); %questi parametri sono per una cella con 71 coppie e 6-ampere module
param.s2=-5.37574*10^(-5);
param.s3=7.42731*10^(-7);
param.s4=-1.27141*10^(-9);
param.r1=2.08317;
param.r2=-1.98763*10^(-2);
param.r3=8.53832*10^(-5);
param.r4=-9.03143*10^(-8);
param.k1=4.76218*10^(-1);
param.k2=-3.89821*10^(-6);
param.k3=-8.64864*10^(-6);
param.k4=2.20869*10^(-8);
e passare la struttura "param" alla funzione "Parametri"
[SM,RM,KM]=Parametri(param,Imax,N,Tc,Th);
nella funzione puoi assegnare alle variabili "
s1", "
s2", "
s3" ecc. i valori dei corrispondenti campi della struttura o usare direttamente la struttura.
La sezione di lettura dei files di input è ridondante e "spreca" memoria":
load(cat(2,files,freq{i},misure));
data=eval(cat(2,files,freq{i},misure));
p=data;
"
load" legge già il file, per cui, la seconda istruzione è inutile (e, inoltre, utilizza la funzione "
eval", cosa assolutamente sconsigliata).
L'istruzione "p=data" è superflua, perchè non scrivere direttamente
"
p=eval(cat(2,files,freq{i},misure))"
Volendo gestire in modo dinamico i nomi dei files di input e risparmiare memoria, puoi leggere i files in questo modo:
%
% Il file di input viene caricato nella struttura "tmp"
%
tmp=load([files,freq{i},misure,'.mat']);
%
% Istruzione non necessari, il file di input viene già letto dalla
% istruzione precedente
%
% % % data=eval(cat(2,files,freq{i},misure));
%
% Identificazione del none del campo della struttura
%
f_name=fieldnames(tmp);
%
% Assegnazione del contenuto del campo alla matrice "p"
%
p=getfield(tmp,f_name{1});
%
% Cancellazione della variabile temporanea "tmp"
%
clear tmp
Prima del loop interno, conviene inizializzare il vettore "
Hv"; non è strettamente necessario, MatLab gestisce l'incremento della dimensione delle struttura dati automaticamente, ma questo rallenta in modo significativo l'esecuzione dello script.
%
% Inizializzazione del vettore "Hv"
%
Hv=zeros(length(p),1);
Nel loop interno, i vettori "
I", "
Th" e "
Tc" non sono utilizzati come tali, quindi, a meno che non servano in altre parti del codice che non hai pubblicato, puoi sostituire "
I(v)", "
Th(v)" e "
Tc(v)", semplicemente con "
I", "
Th" e "
Tc".
Nel caso invece ti servissero come vettori, dovresti inizializzarli prima del loop (per le ragioni sopra espresse).
Analogo discorso per i vettori
SMrec,
RMrec e
KMrec che sembrano inutilizzati (e non inizializzati).
Per quanto riguarda le ultime quattro istruzioni:
fasesens1=((Ti_max-To_maxs1)/(1/frequenza(i)));
fasesens2=((Ti_max-To_maxs2)/(1/frequenza(i)));
fasesens3=((Ti_max-To_maxs3)/(1/frequenza(i)));
fasesens4=((Ti_max-To_maxs4)/(1/frequenza(i)));
all'inizio dello script le variabili "
fasesens1", "
fasesens2", "
fasesens3" e "
fasesens4" sono definite come vettori, ma poi sono usate come variabili "semplici"
fasesens1=zeros(1,16);
fasesens2=zeros(1,16);
fasesens3=zeros(1,16);
fasesens4=zeros(1,16);
Dovresti verificare.
Ribadito che ho trovato l'errore, ma non la soluzione che dipende dall'algoritmo che devi implementare, inserisco di seguito i codici della funzione "Parametri" e dello script principale con le modifiche che ti ho suggerito (QUINDI L'ERRORE "C'E ancora)
OVVIAMENTE DECLINANDO QUALUNQUE RESPONSABILITA sull'uso e sui risultati
Funzione "Parametri"
% function [ SM,RM,KM ] = Parametri( Imax,N,Tc,Th)
function [ SM,RM,KM ] = Parametri(param,Imax,N,Tc,Th)
%
% Aggiunto parametro di input "param" => struttura contenente i valori
% delle costanti "s1", "s2", "s3", ...
%
%PARAMETRI: calcolo del coefficiente di Seedback SM, della resistenza
%elettrica RM e della conduttanza termica KM. I valori di N e Imax servono
%per adattare i parametri ottenuti nel sistema considerato. Tesi pag 12...
%parametri delle equazioni:
%
% Calcolo delle costanti spostato nel "main" in modo da essere fatto una
% volta sola
%
s1=param.s1; %=1.33450*10^(-2); %questi parametri sono per una cella con 71 coppie e 6-ampere module
s2=param.s2; %=-5.37574*10^(-5);
s3=param.s3; %=7.42731*10^(-7);
s4=param.s4; %=-1.27141*10^(-9);
r1=param.r1; %=2.08317;
r2=param.r2; %=-1.98763*10^(-2);
r3=param.r3; %=8.53832*10^(-5);
r4=param.r4; %=-9.03143*10^(-8);
k1=param.k1; %=4.76218*10^(-1);
k2=param.k2; %=-3.89821*10^(-6);
k3=param.k3; %=-8.64864*10^(-6);
k4=param.k4; %=2.20869*10^(-8);
if Tc==Th %DeltaT=0
SM= s1 + s2*Tc + s3*Tc^2 + s4*Tc^3;
RM= r1 + r2*Tc + r3*Tc^2 + r4*Tc^3;
KM= k1 + k2*Tc + k3*Tc^2 + k4*Tc^3;
else
SMth= s1*Th + (s2*Th^2)/2 + (s3*Th^3)/3 + (s4*Th^4)/4;
SMtc= s1*Tc + (s2*Tc^2)/2 + (s3*Tc^3)/3 + (s4*Tc^4)/4;
SM = (SMth-SMtc)/(Th-Tc);
RMth= r1*Th + (r2*Th^2)/2 + (r3*Th^3)/3 + (r4*Th^4)/4;
RMtc= r1*Tc + (r2*Tc^2)/2 + (r3*Tc^3)/3 + (r4*Tc^4)/4;
RM = (RMth-RMtc)/(Th-Tc);
KMth= k1*Th + (k2*Th^2)/2 + (k3*Th^3)/3 + (k4*Th^4)/4;
KMtc= k1*Tc + (k2*Tc^2)/2 + (k3*Tc^3)/3 + (k4*Tc^4)/4;
KM = (KMth-KMtc)/(Th-Tc);
end
%for the configuration considered with N and Imax, the following
%normalization must be done (pag 16 pdf tesi):
SM=SM*N/71;
RM=RM*6*N/(Imax*71);
KM=KM*Imax*N/(6*71);
Script principale
misure='mhz';
files='frequenza';
freq={'05' '07' '09' '1' '2' '3' '5' '6' '7' '8' '9' '10' '20' '30' '40' '50' };
Hvf=[];
frequenza=[0.5, 0.7, 0.9, 1, 2, 3, 5, 6, 7, 8, 9, 10, 20, 30, 40, 50]*10^-3;
modulosensore1=zeros(1,16);
modulosensore2=zeros(1,16);
modulosensore3=zeros(1,16);
modulosensore4=zeros(1,16);
fasesens1=zeros(1,16);
fasesens2=zeros(1,16);
fasesens3=zeros(1,16);
fasesens4=zeros(1,16);
%
% Definizione parametri delle equazioni (from "Parametri.m")
%
param.s1=1.33450*10^(-2); %questi parametri sono per una cella con 71 coppie e 6-ampere module
param.s2=-5.37574*10^(-5);
param.s3=7.42731*10^(-7);
param.s4=-1.27141*10^(-9);
param.r1=2.08317;
param.r2=-1.98763*10^(-2);
param.r3=8.53832*10^(-5);
param.r4=-9.03143*10^(-8);
param.k1=4.76218*10^(-1);
param.k2=-3.89821*10^(-6);
param.k3=-8.64864*10^(-6);
param.k4=2.20869*10^(-8);
%
for i=1:(size(freq,2)); %ciclo atto a scorrere tutte le acquisizioni fatte nella condizione scelta
% % % load(cat(2,files,freq{i},misure));
%
% Il file di input viene caricato nella struttura "tmp"
%
tmp=load([files,freq{i},misure,'.mat']);
%
% Istruzione non necessari, il file di input viene già letto dalla
% istruzione precedente
%
% % % data=eval(cat(2,files,freq{i},misure));
%
% Identificazione del none del campo della struttura
%
f_name=fieldnames(tmp);
%
% Assegnazione del contenuto del campo alla matrice "p"
%
p=getfield(tmp,f_name{1});
%
% Cancellazione della variabile temporanea "tmp"
%
clear tmp
%
% % % Hv=[];
%
% Inizializzazione del vettore "Hv"
%
Hv=zeros(length(p),1);
%
Imax=5;
N=127;
%
% Inizializzazine calcolo durata del loop
%
tic
%
% Verifica della possibilità che "Tc" sia (quasi) uguale a "Th"
%
tmp_Tc=p(:,2);
tmp_Th=p(:,3);
min_diff=min(abs(tmp_Tc-tmp_Th))
soglia=0.0001; % questa istruzione dovrebbe essere spostata fuori dal loop
if(min_diff <= soglia)
disp(['Input file: ' files,freq{i},misure,'.mat'])
disp(['min_diff = ' num2str(min_diff)])
disp([' Possibile Tc == Th non catturato in "Parametri.m"'])
end
%
for v=1:length(p(:,1))
% % % I(v)=p(v,8);
I=p(v,8);
% % % Th(v)=p(v,3); %Valore della temperatura del lato caldo (sensor8)
Th=p(v,3); %Valore della temperatura del lato caldo (sensor8)
% % % Tc(v)=p(v,2); %Valore della temperatura del lato freddo (sensor2)
Tc=p(v,2); %Valore della temperatura del lato freddo (sensor2)
%calculation parameters
% % % [SM,RM,KM]=Parametri(Imax,N,Tc(v),Th(v));
[SM,RM,KM]=Parametri(param,Imax,N,Tc,Th);
% % % SMrec(v)=SM;
% % % RMrec(v)=RM;
% % % KMrec(v)=KM;
%calculation heat flux
% % % Hv(v,1)=SM*I(v)*Th(v)+0.5*RM*I(v)^2-KM*(Th(v)-Tc(v)); %calcolo di Qh pag 15
Hv(v)=SM*I*Th+0.5*RM*I^2-KM*(Th-Tc); %calcolo di Qh pag 15
end
%
% Stampa della durata del loop nella Command Window
%
toc
%
Hvf{i}=Hv;
[max_Hvf,idx_max]=cellfun(@max,Hvf);
Ti_max=p(idx_max,1);
min_Hvf=cellfun(@min,Hvf);
[max_p1,idxT1]= max(p(:,4));
To_maxs1=p(idxT1,1);
[max_p2,idxT2]= max(p(:,5));
To_maxs2=p(idxT2,1);
[max_p3,idxT3]= max(p(:,6));
To_maxs3=p(idxT3,1);
[max_p4,idxT4]= max(p(:,7));
To_maxs4=p(idxT4,1);
modulosensore1(i)=20*log10(max(p(:,4))-max_Hvf)/(min(p(:,4))-min_Hvf);
modulosensore2(i)=20*log10(max(p(:,5))-max_Hvf)/(min(p(:,5))-min_Hvf);
modulosensore3(i)=20*log10(max(p(:,6))-max_Hvf)/(min(p(:,6))-min_Hvf);
modulosensore4(i)=20*log10(max(p(:,7))-max_Hvf)/(min(p(:,7))-min_Hvf);
fasesens1=((Ti_max-To_maxs1)/(1/frequenza(i)));
fasesens2=((Ti_max-To_maxs2)/(1/frequenza(i)));
fasesens3=((Ti_max-To_maxs3)/(1/frequenza(i)));
fasesens4=((Ti_max-To_maxs4)/(1/frequenza(i)));
end
Hope this helps.