function [ I, CC ] = CC_generaIrraggiameto( CC )



%% 1) genera la nuvolosità come un rumore colorato saturato
% ho assunto: 1 per sole pieno, MIN per sole coperto
% di notte o quando la facciata è in ombra, nessuna insolazione ma questo
% viene calcolato successivamente

tempi = CC.inizioNum:CC.dt:CC.fineNum;
lumin = lsim(   zpk([1 .2+0.05*i .2-0.05*i],[.9995 1 .1+0.1*i .1-0.1*i],1,1),...
    randn(1,CC.Nsim));
lumin = lumin - min(lumin); % rendilo positivo
[yy,xx] = hist(lumin,512); % inizia il trattamento di saturazione
YY = cumsum(yy)./length(lumin);
SUP = xx(find(YY > (1-CC.frazioneSolePieno),1));
INF = xx(find(YY > CC.frazioneCoperto,1));
lumin(lumin<INF) = INF;
lumin(lumin>SUP) = SUP;
lumin = (lumin-INF)./(SUP-INF);
lumin = (1-CC.min)*lumin+CC.min;
CC.lumin = lumin;



%% 2 pre) genera i tempi su un giorno

fattore = 1;
tempi2 = CC.giornoMedio:CC.dt*fattore:(CC.giornoMedio+datenum(0,0,0,24,0,0)-CC.dt*fattore);
Nsim2 = length(tempi2);
tempi3 = CC.giornoMedio:CC.dt:(CC.giornoMedio+datenum(0,0,0,24,0,0)-CC.dt);

%% 2) genera il vettore congiungente il sole con un punto sulla superficie
% della terra
k = 1;
for t = tempi2
    v_sole(k,:) = versore_al_sole(t,CC.location,CC.UTC);
    k = k + 1;
end



%% 3) carica la mesh, calcola i baricentri

load('MeshSuperficie.mat');
% NodiSuperficiali = coordinate xyz vertici triangolo
% MeshSuperficiale = indici ai vettori dei vertici, 3 indici fanno un
% triangolo

% -- ROTAZIONE --
R = [];
load R
NodiSuperficiali = NodiSuperficiali*R';
NormaleFacce = NormaleFacce*R';
% -- FINE ROTAZIONE --

CC.Ntri = length(MeshSuperficiale);

for k = 1:CC.Ntri
    bar(k,:) = mean(NodiSuperficiali(MeshSuperficiale(k,:),:));
end



%% 3 bis) Aggiungi alberi come ambiente

ix = 1;
for y = -3:2
    for z = 0:9
        amb(ix,:) = [8 y z];
        ix = ix + 1;
    end
end



%% 4) calcola la luminosità con ombre incluse

% velocizza i conti
I2 = zeros(CC.Ntri,Nsim2,'single');
v_sole = single(v_sole');
bar = single(bar);
amb = single(amb);
NormaleFacce = single(NormaleFacce);

% quelli sottoterra sono sempre in ombra totale
CC.sempreInOmbra = find((  bar(:,3) <= 0  ));
insiemeIndici = setdiff(1:CC.Ntri,CC.sempreInOmbra);
I2(CC.sempreInOmbra,:) = zeros(length(CC.sempreInOmbra),Nsim2);

wb = waitbar(0,'Calcolo ombre');

for k = insiemeIndici
    tic
    waitbar(k/CC.Ntri,wb);
    % vettRaggi unisce il baricentro corrente con tutti gli altri baricentri
    % tranne se stesso
    vettRaggi = [   ( bar([(1:k-1) (k+1):end],:)-repmat(bar(k,:),length(bar)-1,1) ) ;
        ( amb - repmat(bar(k,:),length(amb),1) ) ];
    % vettRaggi è normalizzato
    vettRaggi = vettRaggi./repmat(sqrt(sum(vettRaggi.^2,2)),1,3);
    % calcolo il prodotto scalare che corrisp al coseno essendo
    % normalizzati
    coseno = vettRaggi*v_sole;
    % l'irraggiamento è proporzionale al coseno tra la normale e il versore
    % al sole ammesso che
    % ~any(coseno > BC.cosenoOmbra): non sia oscurato da altri baricentri
    % (v_sole(:,3) > 0): non sia notte
    I2(k,:) = (NormaleFacce(k,:)*v_sole).*...
        ~any(coseno > CC.cosenoOmbra).*(v_sole(3,:) > 0);
    toc
end
% se l'angolo tra normale e sole è negativo, allora è in ombra
I2(I2 < 0) = 0;
close(wb)

%%
if CC.show == 1
    tempo = 1;
    VisualizzaSurf2(NodiSuperficiali,MeshSuperficiale,bar,NormaleFacce,I2(:,tempo),v_sole(tempo,:),amb,campos,axis);
    title(datestr(tempi2(tempo)))
    drawnow
    axis equal
    for tempo = 1:length(tempi2)
        cp = campos;
        ax = axis;
        VisualizzaSurf2(NodiSuperficiali,MeshSuperficiale,bar,NormaleFacce,I2(:,tempo),[0 0 0],amb,cp,ax);
        title(datestr(tempi2(tempo)))
        colormap hot
        colorbar
        caxis([0 1])
        campos(cp)
        axis(ax)
        drawnow
    end
end

%% 5) Interpola (se fattore > 1) e affianca per il periodo necessario e annuvola

I = zeros(CC.Ntri,CC.Nsim,'single');
for k = 1:CC.Ntri
    I(k,:) = repmat(interp1(tempi2,I2(k,:),tempi3),1,CC.giorni).*CC.lumin';
end




%%

    function v = versore_al_sole( time, location, UTC )
        %VERSORE_AL_SOLE versore verso la posizione del sole
        %   USCITA
        %   un versore [x y z] con convenzione
        %       NORD   = asse Y
        %       EST    = asse X
        %       ZENITH = asse Z
        %   INGRESSI
        %   location come in sun_position
        %   time � in format matlab (vedi datenum)
        %   UTC � la correzione da apportare (in format matlab)
        %       in Italy, in inverno siamo a UTC+1 e in estate a UTC+2
        %       quindi fornire UTC=1 in inverno e UTC=2 in estate
        %       il programma provvede a effettuare la correzione

        sun = sun_position(time-UTC, location);
        sun.azimuth = sun.azimuth*pi/180;
        sun.elevation = (90-sun.zenith)*pi/180;
        v = [
            cos(sun.elevation)*sin(sun.azimuth);
            cos(sun.elevation)*cos(sun.azimuth);
            sin(sun.elevation)
            ];
    end

    function VisualizzaSurf2(p,t,bar,nor,I,v,amb,cp,ax)
        hold off
        h=trisurf(t,p(:,1),p(:,2),p(:,3),I);
        hold on
        quiver3(bar(:,1),bar(:,2),bar(:,3),nor(:,1),nor(:,2),nor(:,3))
        plot3(amb(:,1),amb(:,2),amb(:,3),'+')
        v = v*10;
        quiver3(0,0,0,v(1),v(2),v(3))

    end


end
