\documentclass[oldfontcommands,a4paper]{memoir}
\usepackage[T1]{fontenc}
\usepackage[utf8]{inputenc}
\usepackage[a4paper]{geometry}
\geometry{verbose,tmargin=2cm,bmargin=2cm,lmargin=2.5cm,rmargin=2.5cm}

%\usepackage{color}
\usepackage[usenames,dvipsnames]{xcolor}
\usepackage{mathpazo}%Letra palatino con fuentes para matemáticas

\pagestyle{Ruled}
\usepackage{memhfixc}
\usepackage{mempatch}

\raggedbottom
\sloppybottom
\clubpenalty=10000
\widowpenalty=10000

%\raggedbottomsection
\feetbelowfloat

\captiontitlefont{\itshape}
\captionnamefont{\scshape}
%\captionstyle{\centering}
\hangcaption


\usepackage{url}
\usepackage{hyperref}
\hypersetup{
    bookmarks=true,         % show bookmarks bar?
    unicode=true,          % non-Latin characters in Acrobat’s bookmarks
    bookmarksnumbered=false,
    bookmarksopen=false,
    breaklinks=true,
    backref=true,
    pdftoolbar=true,        % show Acrobat’s toolbar?
    pdfmenubar=true,        % show Acrobat’s menu?
    pdffitwindow=false,     % window fit to page when opened
    pdfstartview={FitH},    % fits the width of the page to the window
    pdftitle={Introduction to solaR},    % title
    pdfauthor={Oscar Perpiñán Lamigueiro},     % author
    pdfsubject={Photovoltaic Solar Energy},   % subject of the document
    pdfcreator={AucTeX/Emacs},   % creator of the document
    pdfproducer={LaTeX}, % producer of the document
    pdfkeywords={solar radiation, photovoltaics, R}, % list of keywords
    pdfnewwindow=true,      % links in new window
    pdfborder={0 0 0},
    colorlinks=true,       % false: boxed links; true: colored links
    linkcolor=Brown,          % color of internal links
    citecolor=BrickRed,        % color of links to bibliography
    filecolor=black,      % color of file links
    urlcolor=Blue           % color of external links 
}

\pretitle{\vfill \begin{flushright} \bfseries \scshape \HUGE \color{BrickRed}}
\posttitle{\par\end{flushright}}

\preauthor{\begin{flushright} \large \scshape}
\postauthor{\par\end{flushright}}

%\date{}
\predate{\vfill \begin{flushright}\large\scshape}
\postdate{\par\end{flushright}\vfill}

\setsecnumdepth{subsection}


% \definecolor{ared}{rgb}{.647,.129,.149}
% \renewcommand{\colorchapnum}{\color{ared}}
% \renewcommand{\colorchaptitle}{\color{ared}}
% \chapterstyle{pedersen}
\chapterstyle{ger}

\setlength{\afterchapskip}{35pt}
\maxtocdepth{section}

%\setcounter{topnumber}{3}
%\setcounter{bottomnumber}{2}
%\setcounter{totalnumber}{4}
\renewcommand{\topfraction}{0.85}
\renewcommand{\bottomfraction}{0.5}
\renewcommand{\textfraction}{0.15}
\renewcommand{\floatpagefraction}{0.7}


%Centra las figuras en los flotantes y los enmarca
\makeatletter
\renewenvironment{figure}[1][]{%
     	\@float{figure}%
		%\begin{framed}    
		\precaption{\rule{\linewidth}{0.4pt}\par}%En las figuras el caption va debajo
		%\hrule\vspace{\onelineskip}
		\centering
		  }{%
		%\end{framed}
		%\postcaption{\rule{\linewidth}{0.4pt}}
		%\vspace{\onelineskip}\hrule
    	\end@float	
}
\makeatother

\makeatletter
\renewenvironment{table}[1][]{%
      	\@float{table}%
		%\begin{framed}    
		\postcaption{\rule{\linewidth}{0.4pt}\par}%En las tablas el caption va encima
		\centering
		  }{%
		%\end{framed}
    	\end@float	
}
\makeatother


\renewcommand{\textfloatsep}{10pt}%Espacio entre el flotante y el texto

 
\makeatletter
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Textclass specific LaTeX commands.
\usepackage[noae]{Sweave}
\newcommand{\Rcode}[1]{{\texttt{#1}}}
\newcommand{\Robject}[1]{{\texttt{#1}}}
\newcommand{\Rcommand}[1]{{\texttt{#1}}}
\newcommand{\Rfunction}[1]{{\texttt{#1}}}
\newcommand{\Rfunarg}[1]{{\textit{#1}}}
\newcommand{\Rpackage}[1]{{\textit{#1}}}
\newcommand{\Rmethod}[1]{{\textit{#1}}}
\newcommand{\Rclass}[1]{{\textit{#1}}}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% User specified LaTeX commands.
%\VignetteIndexEntry{Introduction to solaR}
\usepackage{flafter}
\usepackage{boxedminipage}
\renewenvironment{Schunk}{\begin{center}
    \scriptsize
    \begin{boxedminipage}{0.95\textwidth}}{
    \end{boxedminipage}\end{center}}

\usepackage{siunitx}
\sisetup{per=fraction,fraction=nice, decimalsymbol=comma}
%\usepackage{lscape}



\newunit{\wattpeak}{Wp}
\newunit{\watthour}{Wh}
\newunit{\amperehour}{Ah}
\newunit{\celula}{celula}

\makeatother

\begin{document}

\setkeys{Gin}{width=0.5\textwidth}
\begin{titlingpage}

\title{Introduction to \emph{solaR}}

\author{Oscar Perpiñán Lamigueiro}


\date{06 June 2011}


\maketitle


\end{titlingpage}

\frontmatter

%%\chapter{Introduction}
\chapterprecis{\vfill{}}

\rule[.5ex]{\linewidth}{1pt} 

The  \texttt{solaR} package includes a set of functions which calculate
the solar radiation incident on a photovoltaic generator and simulate the 
performance of several applications of the photovoltaic energy \cite{Perpinan2011}.
This package performs the whole calculation from both \emph{daily} and 
\emph{intra-daily} global horizontal irradiation to the final productivity of 
grid connected PV systems and water pumping PV systems. 
Besides, the package includes several visualization methods based on
the \texttt{lattice} and \texttt{latticeExtra} packages, and tools for the statistical analysis of
the performance of a large PV plant composed of several systems.
The package is constructed with \texttt{S4} classes and methods. 
The time series are constructed with the \texttt{zoo} package \cite{Zeileis.Grothendieck2005}.

Examples of its use are frequently published at
\url{http://procomun.wordpress.com}.

\rule[.5ex]{\linewidth}{1pt} 

\chapterprecis{\vfill{}}

\mainmatter


<<echo=false, results=hide>>= 
library(solaR)
oopt <- options() 
options("digits"=4) 
options("width"=80)
@

\chapter{Solar Geometry}

The apparent movement of the Sun is defined with some equations included in 
the functions  \texttt{fSolD} and \texttt{fSolI}. \texttt{fSolD} computes the daily apparent movement of the Sun from the Earth. This movement is mainly described (for the simulation of photovoltaic systems) by the declination angle, the sunset angle and the daily extra-atmospheric irradiation.
On the other hand, \texttt{fSolI} computes the angles which describe the intra-daily apparent movement of the Sun from the Earth. 

The next example shows these calculations for a certain day:
<<echo=TRUE>>=
#Calculo para un único día 
BTd=fBTd(mode='serie')

lat=37.2
SolD<-fSolD(lat,BTd[100])
SolI<-fSolI(SolD, sample='hour', keep.night=FALSE)
head(SolI)
@

and for a set of days:

<<echo=TRUE>>=
#Calculo para un varios días
SolD<-fSolD(lat,BTd[c(10, 50, 100)])
print(SolD)
@

With the function \texttt{fBTd} it is possible to get time bases with 
different structures. Thus, the calculations for the so called ``average days'' 
need the next piece of code, with the result displayed in the figure \ref{fig:AzimutAltura}. 

<<echo=TRUE>>=
lat=37.2;
SolD<-fSolD(lat,BTd=fBTd(mode='prom'))
SolI<-fSolI(SolD, sample='10 min', keep.night=FALSE)
@

%

\begin{figure}
\begin{centering}
<<echo=TRUE, fig=TRUE>>=
mon=month.abb
p<-xyplot(r2d(AlS)~r2d(AzS),
    groups=month,
    data=SolI, type='l', col='black',
    xlab=expression(psi[s]),ylab=expression(gamma[s]))

plab <- p + glayer({
  idx <- round(length(x)/2+1)
  panel.text(x[idx], y[idx], mon[group.value], pos=3, offset=0.2, cex=0.8)
})

print(plab)

@

\par\end{centering}
\caption{\label{fig:AzimutAltura}Azimuth and height solar angles during the ``average days''.}
\end{figure}


These calculations can also be carried out for the whole year (figure \ref{fig:Declinacion}).

<<echo=TRUE>>=
BTd=fBTd(mode='serie')
solD <- fSolD(lat, BTd)
summary(solD)
@


%
\begin{figure}
\begin{centering}
<<echo=TRUE, fig=TRUE>>=
p <- xyplot(solD$decl)
print(p)
@
\par\end{centering}

\caption{Declination throughout the year\label{fig:Declinacion}}

\end{figure}

These two functions are included in a function,  \texttt{calcSol}. 
This function constructs an object of class  \texttt{Sol} containing in its slots 
the \texttt{zoo} objects created by \texttt{fSolD} and \texttt{fSolI}. 
This class owns methods for getting and displaying information (for example, \texttt{as.zooD, as.zooI, xyplot}).

There are three methods for the sun geometry calculations. These
methods are named as 'cooper' (the default in previous versions),
'spencer', 'michalsky' (the default in this version) and 'strous'. 
<<echo=TRUE>>=
lat=37.2
BTd=fBTd(mode='serie')

sol=calcSol(lat, BTd[100])
print(as.zooD(sol))

solStrous=calcSol(lat, BTd[100], method='strous')
print(as.zooD(solStrous))

solSpencer=calcSol(lat, BTd[100], method='spencer')
print(as.zooD(solSpencer))

solCooper=calcSol(lat, BTd[100], method='cooper')
print(as.zooD(solCooper))

@ 



\chapter{Solar Radiation}

Values of global horizontal irradiation are commonly available either
as monthly averages of daily values or as a time series of daily
during one or several years.  The analysis of the performance of a PV
system starts from the transformation of the global horizontal
irradiation to global, diffuse and direct horizontal irradiance and
irradiation, and then irradiance and irradiation on the generator
surface.


\section{Irradiation and irradiance on the horizontal plane}

The function \texttt{fCompD} extracts the diffuse and direct
components from the daily global irradiation on a horizontal surface
by means of regressions between the clearness index and the diffuse
fraction parameters.  This function need the results from
\texttt{fSolD}, a set of values of global horizontal irradiation
($\si{\watthour\per\meter\squared}$), and the correlation between the
clearness index and the diffuse fraction.  The current version of
\texttt{solaR} includes several correlations (type
\texttt{help(corrFdKt)} for details).  Besides, the user may define a
particular correlation through the argument \texttt{f}.  Once again
for a certain day:

<<echo=TRUE>>=
BTd=fBTd(mode='serie')
SolD<-fSolD(lat, BTd[100])
SolI<-fSolI(SolD, sample='hour')

G0d=zoo(5000, index(SolD))
fCompD(SolD, G0d, corr = "Page") 
fCompD(SolD, G0d, corr = "CPR") 
@ 
and for the ``average days'':

<<echo=TRUE>>=
lat=37.2;
G0dm=c(2.766,3.491,4.494,5.912,6.989,7.742,7.919,7.027,5.369,3.562,2.814,2.179)*1000; 
Rad=readG0dm(G0dm, lat)
solD<-fSolD(lat,fBTd(mode='prom'))
fCompD(solD, Rad, corr = 'Page') 
@ 

Let's use \texttt{corr='user'} to define a function with the
correlation of Collares Pereira and Rabl \cite{Collares-Pereira.Rabl1979}. 
Obviously, we shall obtain the same result as with \texttt{corr='CPR'}.

<<echo=TRUE>>=
fKTd=function(x){(0.99*(x<=0.17))+(x>0.17)*(1.188-2.272*x+9.473*x^2-21.856*x^3+14.648*x^4)}
fCompD(SolD,G0d, corr="user",f=fKTd)
@ 
                      
The daily profile of irradiance is obtained with the function \texttt{fCompI}. 
This function needs the information provided by \texttt{fCompD} and \texttt{fSolI} or \texttt{calcSol}. For example, the profiles for the ``average days'' are obtained with the next code (fig. \ref{fig:ComponentesIrradiancia}).

<<echo=TRUE>>=
lat=37.2;
sol<-calcSol(lat, fBTd(mode='prom'), sample='hour', keep.night=FALSE)

G0dm=c(2.766,3.491,4.494,5.912,6.989,7.742,7.919,7.027,5.369,3.562,2.814,2.179)*1000;
Ta=c(10, 14.1, 15.6, 17.2, 19.3, 21.2, 28.4, 29.9, 24.3, 18.2, 17.2, 15.2)
BD<-readG0dm(G0dm=G0dm, Ta=Ta, lat=37.2)
compD<-fCompD(sol, BD, corr = 'Page') 
compI<-fCompI(sol, compD)
@

%
\begin{figure}
\begin{centering}
<<echo=TRUE, fig=TRUE>>=
p <- xyplot(G0+B0+D0~w|month, data=compI, type='l',auto.key=list(space='right'))
print(p)
@
\par\end{centering}
\caption{Global, diffuse, and direct irradiance during the ``average days''.\label{fig:ComponentesIrradiancia}}
\end{figure}

\subsection{Meteorological data}

There are several functions to construct a \texttt{Meteo} object with radiation and temperature data.
For daily data, if it is stored in a local file or a  \texttt{data.frame}, the functions  \texttt{readBD} and \texttt{df2Meteo} are recommended, while \texttt{readG0dm} is indicated when only 12 monthly means are available. For intradaily data the correspondent functions are \texttt{readBDi} and \texttt{dfI2Meteo}. Besides, \texttt{zoo2Meteo} can construct a \texttt{Meteo} object from a \texttt{zoo} object both for daily and intradaily data.

For example, the \texttt{helios} dataset included in the package, obtained from \url{http://helios.ies-def.upm.es}, can be converted to a \texttt{Meteo} object with the next code:
<<echo=TRUE>>=
data(helios)
names(helios)=c('date', 'G0', 'TempMax', 'TempMin')

bd=df2Meteo(helios, dates.col='date', lat=41, source='helios-IES', format='%Y/%m/%d')

summary(getData(bd))
@ 

On the other hand, the function  \texttt{readSIAR} is able to download the meteorological data 
available at \url{www.marm.es/siar}.
This webpage provides daily measurements from a set of agroclimatic stations 
located in Spain. This function needs the code of the station and its province, 
and the start and end date. The codes of stations and provinces are stored at the
dataset \texttt{RedEstaciones}. For example, there are several stations in Madrid:

<<echo=TRUE>>=
data(RedEstaciones)
Madrid<-subset(RedEstaciones,NomProv=='Madrid')
print(Madrid)
@

\texttt{readSIAR} constructs an object of class \texttt{Meteo}. The data is obtained with 
the method \texttt{getData}.  If only the irradiation series is needed, the method  \texttt{getG0} 
is recommended.

For example, let's obtain the 2009 data from the station at Aranjuez (fig. \ref{fig:Aranjuez}). 
It is important to note that the radiation measurements available at the webpage are
in  $\si{\mega\joule\per\meter\squared}$, but \texttt{readSIAR} converts the values to $\si{\watthour\per\meter\squared}$:

<<echo=TRUE>>=
Aranjuez<-readSIAR(28,3,'01/01/2009','31/12/2009')
@

\begin{figure}
  \centering
<<echo=TRUE, fig=TRUE>>=
p=xyplot(G0~TempMedia|month, data=Aranjuez, type=c('p', 'r'))
print(p)
@ 
  \caption{Daily irradiation and mean temperature in the station of Aranjuez.}
  \label{fig:Aranjuez}
\end{figure}

This database includes information of maximum and minimum values of temperature. 
The function \texttt{fTemp} calculates a profile of the ambient temperature with this information
following the method proposed in \cite{Huld.Suri.ea2006}. 
The evolution of this synthetic temperature during March is displayed in the figure \ref{fig:Ta}.

<<echo=TRUE>>=
lat=41; 
sol=calcSol(lat, BTd=indexD(Aranjuez), sample='hour')
Temp<-fTemp(sol,Aranjuez)
@ 

%
\begin{figure}
\begin{centering}
<<echo=TRUE, fig=TRUE>>=
wTemp=window(Temp, start=as.POSIXct('2009-03-01'), end=as.POSIXct('2009-03-31'))
p=xyplot(wTemp, col='black', ylab='T')+layer_(panel.xblocks(x, DoY, col=c('lightgray', 'white')))
print(p)
@
\par\end{centering}
\caption{Evolution of the ambiente temperature during March 2009 in Aranjuez.\label{fig:Ta}}
\end{figure}



\subsection{The function \texttt{calcG0}}

The previous steps are included in the function  \texttt{calcG0}. For example, with the next code, 
the components of horizontal irradiation and irradiance are obtained from the measurements of the
meteorological station of Aranjuez (figure \ref{fig:calcG0}).

<<>>=
g0<-calcG0(lat=37.2,
	modeRad='siar',
	dataRad=list(prov=28,est=3,
	start='01/01/2009', end='31/12/2009'))
print(g0)
@

%
\begin{figure}
\begin{centering}
<<echo=TRUE, fig=TRUE>>=
p=xyplot(g0)
print(p)
@
\par\end{centering}
\caption{Components of horizontal irradiation calculated with \texttt{calcG0}.\label{fig:calcG0}}
\end{figure}

\texttt{solaR} accepts intradaily irradiation data sources. For example, the Measurement and Instrumentation Data Center of the NREL (NREL-MIDC) provides meteorological data from a variety of stations. We will try the \emph{La Ola - Lanai} station at Hawaii (\url{http://www.nrel.gov/midc/la_ola_lanai/}).
<<echo=TRUE>>=
file='http://www.nrel.gov/midc/apps/plot.pl?site=LANAI&start=20090722&edy=19&emo=11&eyr=2010&zenloc=19&year=2010&month=11&day=1&endyear=2010&endmonth=11&endday=19&time=1&inst=3&inst=4&inst=5&inst=10&type=data&first=3&math=0&second=-1&value=0.0&global=-1&direct=-1&diffuse=-1&user=0&axis=1'

dat <- read.table(file, header=TRUE, sep=',')

lat=20.77
lon=-156.9339
@ 

First, we have to change the names of the columns and calculate the horizontal direct irradiation, since only the normal direct irradiation is included in the file.
<<>>=
names(dat) <- c('date', 'hour', 'G0', 'B', 'D0', 'Ta')
dat$B0 <- dat$G0-dat$D0
@ 

The datalogger program runs using Greenwich Mean Time (GMT), and data is converted to Hawaiian Standard Time (HST) after data collection. With \texttt{local2Solar} we can calculate the Mean Solar Time of the index.
<<>>=
idxLocal <- with(dat, as.POSIXct(paste(date, hour), format='%m/%d/%Y %H:%M', tz='HST'))
idx <- local2Solar(idxLocal, lon=lon)
@ 

Therefore, the  object \texttt{Meteo} is obtained with (figure \ref{fig:NREL-MIDC}):
<<>>=
z <- zoo(dat[,c('G0', 'D0', 'B0', 'Ta')], idx)
NRELMeteo <- zoo2Meteo(z, lat=lat)
@ 

\begin{figure}
  \centering
<<echo=TRUE, fig=TRUE>>=
p <- xyplot(NRELMeteo)
print(p)
@ 
  \caption{1-min irradiation data from NREL-MIDC}
  \label{fig:NREL-MIDC}
\end{figure}

With this data, a \texttt{G0} object can be calculated. First, the direct and diffuse components of the data are used (\texttt{corr='none'}):

<<>>=
g0NREL <- calcG0(lat=lat, modeRad='bdI', dataRad=NRELMeteo, corr='none')
@ 

If these components were not available, a fd-kt hourly correlation is needed:

<<>>=
g0BRL <- calcG0(lat=lat, modeRad='bdI', dataRad=NRELMeteo, corr='BRL')
@ 

% \begin{figure}
%   \centering
% <<echo=TRUE, fig=TRUE>>=
% p <- xyplot(fd~kt, data=g0BRL, pch=19, alpha=0.3, cex=0.5)
% print(p)
% @ 
%   \caption{Diffuse fraction and clearness index correlation following the BRL model, with data from NREL-MIDC}
%   \label{fig:FdKtBRL}
% \end{figure}



\section{Irradiation and irradiance on the generator plane}

The solar irradiance incident on an inclined surface can be calculated
from the direct and diffuse irradiance on a horizontal surface, and
from the evolution of the angles of the Sun and the surface. The
transformation of the direct radiation is straightforward since only
geometric considerations are needed. However, the treatment of the
diffuse irradiance is more complex since it involves the modelling of
the atmosphere. There are several models for the estimation of diffuse
irradiance on an inclined surface. The one which combines simplicity
and acceptable results is the proposal of Hay and McKay. This model
divides the diffuse component in isotropic and anisotropic whose
values depends on a anisotropy index.  On the other hand, the
effective irradiance, the fraction of the incident irradiance that
reaches the cells inside a PV module, is calculated with the losses
due to the angle of incidence and dirtiness. This behaviour can be
simulated with a model proposed by Martin and Ruiz requiring
information about the angles of the surface and the level of dirtiness
\cite{Martin.Ruiz2001}.

The orientation, azimuth and incidence angle are calculated from the
results of \texttt{fSolI} or \texttt{calcSol} with the functions
\texttt{fTheta} and \texttt{fInclin}. These functions can calculate
the movement and irradiance for fixed systems, and two-axis and
horizontal N-S trackers. Besides, the movement of a horizontal NS
tracker due to the backtracking strategy \cite{Panico.Garvison.ea1991}
can be calculated with information about the tracker and the distance
between the trackers included in the system.

Both functions are integrated in \texttt{calcGef}, which construct an
object of class \texttt{Gef}. Once again, this class owns methods for
obtaining and displaying information.

For example, with the previous results, we can calculate the
irradiance and irradiation on a fixed surface. The figure
\ref{fig:GefG0vsTheta} shows the relation between the effective and
incident irradiance versus the cosine of the angle of incidence for
this system.

<<echo=TRUE>>=
gef<-calcGef(lat=37.2, modeRad='prev', dataRad=g0, beta=30)
print(gef)
@ 

%
\begin{figure}
\begin{centering}
<<echo=TRUE, fig=TRUE>>=
p<-xyplot(Gef/G~cosTheta|month, data=gef, type=c('p', 'smooth'), cex=0.4, alpha=0.5)
print(p)
@
\par\end{centering}
\caption{Relation between the effective and incident irradiance versus the cosine of the angle of incidence for a fixed system.\label{fig:GefG0vsTheta}}
\end{figure}

The next lines of code calculate the movement of a N-S horizontal axis
tracker with \emph{backtracking} (\texttt{modeShd='bt'}) and whose inclination angle is
limited to $\ang{60}$ (\texttt{betaLim=60}). The evolution of the inclination angle is
displayed in the figure \ref{fig:Backtracking}. The meaning of the
\texttt{distances} and \texttt{struct} arguments will be detailed in the
\ref{sec:shadows} section.

<<echo=FALSE>>=
G0dm=c(2766, 3491, 4494, 5912, 6989, 7742, 7919, 7027, 5369, 3562, 2814,
2179)
Ta=c(10, 14.1, 15.6, 17.2, 19.3, 21.2, 28.4, 29.9, 24.3, 18.2, 17.2, 15.2)
prom=list(G0dm=G0dm, Ta=Ta)
@ 
<<echo=TRUE>>=
structHoriz=list(L=4.83);
distHoriz=data.frame(Lew=structHoriz$L*4,H=0);

gefBT=calcGef(lat=37.2, dataRad=prom, sample='10 min',
  modeTrk='horiz',
  modeShd='bt', betaLim=60,
  distances=distHoriz,
  struct=structHoriz)
@ 

\begin{figure}
\begin{centering}
<<echo=TRUE, fig=TRUE>>=
p <- xyplot(r2d(Beta)~r2d(w),
     data=gefBT,
     type='l',
     xlab=expression(omega),
     ylab=expression(beta))
print(p)
@
\par\end{centering}

\caption{Evolution of the angle of inclination of a NS horizontal axis tracker with \emph{backtracking} and limitation of angle. \label{fig:Backtracking}}

\end{figure}


\chapter{Productivity of a Grid Connected PV System}

From the previous irradiance calculations, the function \texttt{fProd}
simulates the performance of a Grid Connected PV (GCPV) system paying attention
to some parameters of the system (characteristics of the PV module and the inverter, 
the electrical arrangement of the PV generator, and the losses of the system).

For example, the electrical power, voltage and current of a certain PV system is calculated below.

<<echo=TRUE>>=
inclin=data.frame(Gef=c(200,400,600,800,1000),Ta=25)
fProd(inclin)
@

First, \texttt{fProd} computes the Maximum Power Point (MPP) of the generator
(\texttt{Vmpp} and \texttt{Impp}) at the irradiance and ambient temperature conditions 
contained in \texttt{Inclin}. Next, it checks that this point is inside the MPP window
of the inverter, as defined by  \texttt{inverter\$Vmin} and \texttt{inverter\$Vmax}. 
If the MPP value is outside this range, the function asigns the limit value to the voltage,
and calculates the correspondent current value with a warning. 

Anyway, the inverter input voltage and current are  \texttt{Vdc} e
\texttt{Idc}. With the next piece of code, the \texttt{Vdc} value is
set to \texttt{Vmin} (the minimum value of the MPP window of the
inverter), $\SI{420}{\volt}$,  since \texttt{Vmpp} is below this value.

<<echo=TRUE>>=
inclin=data.frame(Gef=800,Ta=30)
gen1 = list(Nms = 10, Nmp = 11)
inv1=list(Ki = c(0.01,0.025,0.05), Pinv=25000, Vmin=420, Vmax=750, Gumb=20)
prod=fProd(inclin, generator=gen1, inverter=inv1)
print(prod)
@

For this configuration, the losses due to the voltage limitation are:

<<echo=TRUE>>=
with(prod,Vdc*Idc/(Vmpp*Impp))
@

The function \texttt{prodGCPV} integrates the calculation procedure of 
irradiation, irradiance and simulation of the GCPV system. It constructs an object
of class \texttt{ProdGCPV}. 

The next code computes the productivity of the previous GCPV system
working as fixed, NS horizontal axis tracking and two-axis tracking
systems.  The parameters of the generator, module, inverter and rest
of the system are those by default in \texttt{prodGCPV}. The
comparative of the intradaily power time series is shown at the figure
\ref{fig:ComparativaSeguimiento}. Later on, the \texttt{compare} and
\texttt{compareLosses} methods will be shown. They are useful for
comparisons of \emph{yearly} values.

<<echo=FALSE>>=
lat=37.2;
G0dm=c(2766, 3491, 4494, 5912, 6989, 7742, 7919, 7027, 5369, 3562, 2814,
2179)
Ta=c(10, 14.1, 15.6, 17.2, 19.3, 21.2, 28.4, 29.9, 24.3, 18.2, 17.2, 15.2)
prom=list(G0dm=G0dm, Ta=Ta)
@ 
<<echo=TRUE>>=
ProdFixed<-prodGCPV(lat=lat, dataRad=prom, keep.night=FALSE)
Prod2x<-prodGCPV(lat=lat, dataRad=prom, modeTrk='two', keep.night=FALSE)
ProdHoriz<-prodGCPV(lat=lat, dataRad=prom, modeTrk='horiz', keep.night=FALSE)
@

\begin{figure}
\begin{centering}
<<echo=TRUE, fig=TRUE, width=8>>=
ComparePac<-CBIND(two=as.zooI(Prod2x)$Pac, 
     horiz=as.zooI(ProdHoriz)$Pac, 
     fixed=as.zooI(ProdFixed)$Pac)
AngSol=as.zooI(as(ProdFixed, 'Sol'))
ComparePac=CBIND(AngSol, ComparePac)
mon=month(index(ComparePac))

p=xyplot(two+horiz+fixed~AzS|mon, data=ComparePac,
     type='l', auto.key=list(space='right', lines=TRUE, points=FALSE),ylab='Pac')
print(p)
@
\par\end{centering}

\caption{Comparative of intradaily power between tracker strategies.\label{fig:ComparativaSeguimiento}}
\end{figure}

\section{Using \texttt{mergesolaR}}
\label{sec:mergesolar}

The \texttt{mergesolaR} method is designed to merge \emph{daily} time series
of several \texttt{solaR} objects. 

For example, we can obtain the daily irradiation of the whole set of
meteorological stations of Madrid (Spain) and use this information to
calculate the productivity of a grid connected PV system. It is
possible to complete this process with the \texttt{lapply}
function. Therefore we obtain a list of \texttt{ProdGCPV} objects:
<<>>=
EstMadrid <- subset(RedEstaciones, NomProv=='Madrid')
nEstMadrid <- nrow(EstMadrid)
namesMadrid <- EstMadrid$NomEst

prodMadrid <- lapply(1:nEstMadrid,
                     function(x){try(prodGCPV(lat=41, modeRad='siar',
                                                  dataRad=list(prov=28, est=x,
                                                    start='01/01/2009', end='31/12/2010'))
                                     )})
names(prodMadrid) <- namesMadrid
okMadrid <- lapply(prodMadrid, class)!='try-error'
prodMadrid <- prodMadrid[okMadrid]
@ 

In order to prevent from the erroneous behaviour of some stations, the
code includes the use of \texttt{try}. Now it's time for
\texttt{mergesolaR}. Since we have a list of objects, \texttt{do.call}
can solve the problem:
<<>>=
YfMadrid <- do.call(mergesolaR, prodMadrid)
summary(YfMadrid)
@

The \texttt{mergesolaR} for a set of \texttt{ProdGCPV} objects merges
the daily time series of the \texttt{Yf} variable of each object. The
result is a multivariate \texttt{zoo} object where each column is the
daily productivity with the radiation data of each meteorological
station. It can be displayed (for example) with the
\texttt{horizonplot} function (figure \ref{fig:horizonMerge}). This
result will be revisited with the \texttt{Target Diagram} tool (figure
\ref{fig:TargetMerge}). 

\begin{figure}
  \centering
<<fig=TRUE, echo=TRUE>>=
print(
      horizonplot(YfMadrid-rowMeans(YfMadrid),
            origin=0,
            scales=list(y=list(relation='same')),
            colorkey=TRUE)
      )
@ 
  \caption{\texttt{Horizonplot} of the result of a \texttt{mergesolaR}
    call. Previously, the row mean is substracted from each column in
    order to show the deviation of each meteorological station from
    the daily mean of the set.}
  \label{fig:horizonMerge}
\end{figure}





\section{Shadows}
\label{sec:shadows}

The shadows on PV generators alter the performance of the PV generators and
reduce their productivity \cite{Perpinan2008}.
This package includes functions for the estimation of mutual shadows between generators 
from a same system. \texttt{fSombra2X}, \texttt{fSombraHoriz}, \texttt{fSombraEst},
calculate the shadows in two-axis, horizontal axis and fixed systems, respectively.
The function \texttt{fSombra6} is indicated for groups of 6 two-axis trackers. 
Finally, \texttt{fSombra} is a wrapper to the previous functions.

For example, the shadows factor of a tracker surrounded by five
trackers is calculated in the next code box. The dimensions of the
tracker structure and the configuration (rows and columns) of the
plant are defined by \texttt{struct}, while the distances between the
trackers are defined by \texttt{distances}.  The figure
\ref{fig:Sombreado2x} shows the evolution of the shadows factor during
the day (X axis) and year (Y axis).


<<echo=TRUE>>=
##I use method='spencer' and EoT='FALSE' in order to get w equispaced
##so levelplot method produces a suitable plot.
sol<-calcSol(lat=37.2, fBTd(mode='serie'), sample='10 min', 
             EoT=FALSE, keep.night=FALSE, method='spencer')
angGen<-fTheta(sol, beta=35);
Angles=CBIND(as.zooI(sol), angGen)

distances=data.frame(Lew=40,Lns=30,H=0)
struct=list(W=23.11, L=9.8, Nrow=2, Ncol=8)
ShdFactor<-fSombra6(Angles, distances, struct, prom=FALSE)
Angles$FS=ShdFactor
@                               


\begin{figure}
  \centering
<<echo=TRUE, fig=TRUE>>=
p <- levelplot(FS~w*day, data=Angles, par.settings=custom.theme(region=brewer.pal('YlOrBr', n=9)))
print(p)
@                               
\caption{Shadows in a PV plant with two-axis trackers.\label{fig:Sombreado2x}}
\end{figure}


Since the  \texttt{data.frame} \texttt{distances} does only have one row, 
the function \texttt{fSombra6} builds a symmetric grid around the point (0,0,0), 
which is the affected tracker. This grid can also be constructed with:

<<echo=TRUE>>=
distances=data.frame(Lew=c(-40,0,40,-40,40),Lns=c(30,30,30,0,0),H=0)
ShdFactor2<-fSombra6(Angles, distances, struct,prom=FALSE)

identical(coredata(ShdFactor), coredata(ShdFactor2))
@

Besides, \texttt{distances} can define a irregular grid around the affected tracker.
Since this tracker is situated at (0,0,0), \texttt{distances} must have five rows. 
When \texttt{prom=TRUE}, \texttt{fSombra6} provides a weighted averaged of the shadows
in the whole set of trackers, whose distribution in the PV plant is defined 
by \texttt{Nrow} and \texttt{Ncol}.

These functions are integrated in \texttt{calcShd}, \texttt{calcGef} and \texttt{prodGCPV}, 
as these examples show.

First, a two-axis tracking system.
<<echo=TRUE>>=
struct2x=list(W=23.11, L=9.8, Nrow=2, Ncol=8)
dist2x=data.frame(Lew=40, Lns=30, H=0)

prod2xShd<-prodGCPV(lat=lat, dataRad=prom, modeTrk='two', 
    modeShd='area', struct=struct2x, distances=dist2x)
@ 

Then, a N-S horizontal axis tracking system without backtracking, 
<<echo=TRUE>>=
#Horizontal N-S tracker
structHoriz=list(L=4.83);
distHoriz=data.frame(Lew=structHoriz$L*4,H=0);

#Without Backtracking
prodHorizShd<-prodGCPV(lat=lat, dataRad=prom, sample='10 min', 
    modeTrk='horiz',
    modeShd='area', betaLim=60,
    distances=distHoriz,
    struct=structHoriz)
@ 

and a N-S horizontal axis tracking system with backtracking, 
<<echo=TRUE>>=
#With Backtracking
prodHorizBT<-prodGCPV(lat=lat, dataRad=prom, sample='10 min', 
    modeTrk='horiz',
    modeShd='bt', betaLim=60,
    distances=distHoriz,
    struct=structHoriz)
@

Finally, we can compare the \emph{yearly} performance of these systems
with the method \texttt{compare} (fig. \ref{fig:compare}), and
calculate and compare their \emph{yearly} losses with the methods
\texttt{losses} and \texttt{compareLosses}
(fig. \ref{fig:compareLosses}), respectively.
 
\begin{figure}
  \centering
<<echo=TRUE, fig=TRUE>>=
comp <- compare(ProdFixed, Prod2x, ProdHoriz, prod2xShd, prodHorizShd, prodHorizBT)
head(comp)
@   
  \caption{Comparison of several \texttt{ProdGCPV} objects.}
  \label{fig:compare}
\end{figure}

\begin{figure}
  \centering
<<echo=TRUE, fig=TRUE>>=
compL <- compareLosses(ProdFixed, Prod2x, ProdHoriz, prod2xShd, prodHorizShd, prodHorizBT)
head(compL)
@
  \caption{Comparison of the losses of several \texttt{ProdGCPV} objects.}
  \label{fig:compareLosses}
\end{figure}


\section{Position of trackers in a PV plant}

The optimum distance between trackers or static structures of 
a PV grid connected plant depends on two main factors: 
the ground requirement ratio (defined as the ratio of the total ground
area to the PV generator area), 
and the productivity of the system including shadow losses. 
Therefore, the optimum separation may be the one which achieves the highest productivity 
with the lowest ground requirement ratio (GRR). 
However, this definition is not complete since the terrain characteristics 
and the costs of wiring or civil works could alter the decision. 

The function \texttt{optimShd} is a help for choosing this distance: 
it computes the productivity for a set of combinations of distances 
between the elements of the plant \cite{Perpinan2008}. The designer should adopt the decision
from these results with the adequate economical translations.

<<echo=FALSE>>=
lat=37.2;
G0dm=c(2766, 3491, 4494, 5912, 6989, 7742, 7919, 7027, 5369, 3562, 2814,
2179)
Ta=c(10, 14.1, 15.6, 17.2, 19.3, 21.2, 28.4, 29.9, 24.3, 18.2, 17.2, 15.2)
prom=list(G0dm=G0dm, Ta=Ta)
@

Let's analyse the configuration of a PV plant with NS horizontal axis trackers,
without  \emph{backtracking}, and a height of $\SI{4.83}{\meter}$. 
We are interested in a range of separations of 2 and 5 times this dimension. 
Besides, the analysis will be carried out with a limitation in the angle of inclination:

<<echo=TRUE>>=
structHoriz=list(L=4.83);
distHoriz=list(Lew=structHoriz$L*c(2,5));

Shd12Horiz<-optimShd(lat=lat, dataRad=prom,
                     modeTrk='horiz',
                     betaLim=60,
                     distances=distHoriz, res=2,
                     struct=structHoriz,
                     modeShd='area',
                     prog=FALSE)
@

The function \texttt{optimShd} constructs an object of class \texttt{Shade}
This class owns a S4 method of \texttt{plot} for displaying the results (figure \ref{fig:optimHoriz}).

%
\begin{figure}
\begin{centering}
<<echo=TRUE, fig=TRUE>>=
shadeplot(Shd12Horiz)
@
\par\end{centering}

\caption{Mutual shadows in a NS horizontal axis tracking PV system.\label{fig:optimHoriz}}

\end{figure}


Now, for a fixed system (figure \ref{fig:optimEst}):

<<echo=TRUE>>=
structFixed=list(L=5);
distFixed=list(D=structFixed$L*c(1,3));
Shd12Fixed<-optimShd(lat=lat, dataRad=prom,
                     modeTrk='fixed',
                     distances=distFixed, res=1,
                     struct=structFixed,
                     modeShd='area',
                     prog=FALSE)
@

%
\begin{figure}
\centering
<<echo=TRUE, fig=TRUE>>=
shadeplot(Shd12Fixed)
@

\caption{Mutual shadows in a PV plant with fixed structures.\label{fig:optimEst}}

\end{figure}


Last, we are interested in a two-axis tracker whose dimensions are
$\SI{23.11}{\meter}$ width and $\SI{9.8}{\meter}$ height. We will try
to design a PV plant with a grid of trackers of 2 rows and 8 columns.

<<echo=TRUE>>=
struct2x=list(W=23.11, L=9.8, Nrow=2, Ncol=8)
@

We will try the separations between $\SI{30}{\meter}$ and $\SI{50}{\meter}$ for the E-O direction
and between $\SI{20}{\meter}$ and $\SI{50}{\meter}$ for the N-S direction.

<<echo=TRUE>>=
dist2x=list(Lew=c(30,50),Lns=c(20,50))
@

\texttt{optimShd} constructs a sequence from the minimum to the
maximum value of \texttt{distances}, with \texttt{res} as the increment,
in meters, of the sequence. In this example, \texttt{res=5}. 
<<echo=TRUE>>=
ShdM2x<-optimShd(lat=lat, dataRad=prom, modeTrk='two', 
                 modeShd=c('area','prom'), distances=dist2x, struct=struct2x, 
                 res=5, prog=FALSE)
@
Besides, the \texttt{Shade} object includes the
local fitting of the sequence of \texttt{Yf} and \texttt{FS} values
(slots named \texttt{Yf.loess} and \texttt{FS.loess}). The
\texttt{predict} method is used with these \texttt{loess} slots inside
the \texttt{shadeplot} method of the \texttt{Shade} class (figure
\ref{fig:optim2x}).


\begin{figure}
\centering
<<echo=TRUE, fig=TRUE>>=
shadeplot(ShdM2x)
@

\caption{Mutual shadows in a two-axis tracking PV system for a combination of separations between trackers.\label{fig:optim2x}}

\end{figure}



\chapter{PV pumping systems}


\section{Simulation of centrifugal pumps}

The first step for the simulation of the performance of a PV pumping
system (PVPS) is the characterization of the pump under the
supposition of constant manometric height
\cite{Abella.Lorenzo.ea2003}. The function \texttt{fPump} computes the
performance of the different parts of a centrifugal pump fed by a
frequency converter following the affinity laws.

For example, we can characterize the performance of the SP8A44 pump (\url{http://net.grundfos.com/Appl/WebCAPS/InitCtrl?mode=1}) working with $H=\SI{40}{m}$. The information of this pump is stored in the dataset \texttt{pumpCoef}.

<<echo=TRUE>>=
data(pumpCoef) 
CoefSP8A44<-subset(pumpCoef, Qn==8&stages==44)

fSP8A44<-fPump(pump=CoefSP8A44,H=40)
@

The result of \texttt{fPump} is a set of functions which relate the electrical power and 
the flow, hydraulical and mechanical power, and frequency.
These functions allow the calculation of the performance for any electrical power 
inside the range of the pump (figures \ref{fig:EficienciaMotobomba} and \ref{fig:PotenciaMotobomba}):

<<echo=TRUE>>=
SP8A44=with(fSP8A44,{
		Pac=seq(lim[1],lim[2],by=100)
		Pb=fPb(Pac)
		etam=Pb/Pac
		Ph=fPh(Pac)
		etab=Ph/Pb
		f=fFreq(Pac)
		Q=fQ(Pac)
		result=data.frame(Q,Pac,Pb,Ph,etam,etab,f)})


SP8A44$etamb=with(SP8A44,etab*etam)

@

%
\begin{figure}
\centering
<<echo=TRUE, fig=TRUE>>=
lab=c(expression(eta[motor]), expression(eta[pump]), expression(eta[mp]))
p<-xyplot(etam+etab+etamb~Pac,data=SP8A44,type='l', ylab='Eficiencia')
print(p+glayer(panel.text(x[1], y[1], lab[group.number], pos=3)))
@


\caption{Efficiency of the motor and pump for several values of electrical power of a SP8A44 pump with $H=\SI{40}{m}$\label{fig:EficienciaMotobomba}}

\end{figure}


%
\begin{figure}
\centering
<<echo=TRUE, fig=TRUE>>=
lab=c(expression(P[pump]), expression(P[hyd]))
p<-xyplot(Pb+Ph~Pac,data=SP8A44,type='l', ylab='Power (W)', xlab='AC power (W)')
print(p+glayer(panel.text(x[length(x)], y[length(x)], lab[group.number], pos=3)))
@


\caption{\label{fig:PotenciaMotobomba}Mechanical and hydraulical power versus electrical power of a SP8A44 pump with $H=\SI{40}{m}$.}

\end{figure}



\section{Nomograms of PVPS}

The international standard IEC 61725 is of common usage in public licitations of PVPS. 
This standard proposes a equation of the irradiance profile with several parameters 
such as the length of the day, the daily irradiation and the maximum value of the irradiance.
With this profile, the performance of a PVPS can be calculated for several manometric heights and
nominal PV power values. A nomogram can display the set of combinations. 
This graphical tool can help to choose the best combination of pump and PV generator 
for certain conditions of irradiation and height \cite{Abella.Lorenzo.ea2003}. 

This kind of graphics is provided by the function
\texttt{NmgPVPS}. For example, the figure \ref{fig:Nomograma} is a nomogram
for the SP8A44 pump working in a range of heights from 50 to 80
meters, with different PV generators. The peculiar shape of the
curve of 50 meters shows that this pump does not work correctly with
this height.

\texttt{}%
\begin{figure}
\begin{centering}
<<echo=TRUE, fig=TRUE>>=
Pg=seq(3000,5500,by=500);
H=seq(50,80,by=5);
NmgSP8A44<-NmgPVPS(pump=CoefSP8A44, Pg=Pg, H=H, Gd=6000,
     title='Selection of Pumps', theme=custom.theme())
print(NmgSP8A44$plot)
@
\par\end{centering}

\caption{Nomogram for the SP8A44 pump working in a range of heights from 50 to 80 meters, 
with different PV generators.\label{fig:Nomograma} }

\end{figure}


\section{Productivity of PVPS}

A different approach is to simulate the performance of the PVPS
following the same procedure as the one described for the GCPV systems.
The function \texttt{prodPVPS} is the equivalent to the function \texttt{prodGCPV}.
The inputs are very similar between them, although there are some changes 
due to the different composition of the system. This function does not
allow for the calculation of shadows.

Once again with the SP8A44 pump, we compute the flow to be produced by this pump
with a PV generator of $\SI{5500}{\wattpeak}$ and a manometric height of 50 meters. 
The relation between flow and effective irradiance is displayed in the figure \ref{fig:CaudalIrradiancia}. 

<<echo=TRUE>>=
prodSP8A44<-prodPVPS(lat=41,
             modeRad='siar',
             dataRad=list(prov=28,est=3,
               start='01/01/2009', end='31/12/2009'),
             pump=CoefSP8A44, Pg=5500, H=50)
as.zooY(prodSP8A44)
@



\begin{figure}
\begin{centering}
<<echo=TRUE, fig=TRUE>>=
p=xyplot(Q~Gef|month, data=prodSP8A44, 
  cex=0.5, type=c('p', 'smooth'), 
  col.symbol='gray', col.line='black')
print(p)
@
\par\end{centering}

\caption{Flow versus irradiance of a PVPS with a SP8A44 pump and a PV generator with a nominal power of $\SI{5500}{\wattpeak}$ and a manometric height of 50 meters.\label{fig:CaudalIrradiancia}}

\end{figure}

Let's try to obtain more water with this pump using a larger PV
generator of $\SI{7000}{\wattpeak}$.  However, we can check that this
is not a correct decision since the productivity has decreased.  The
figure \ref{fig:CaudalIrradianciaLimite} shows that during the central
months of the year, during the maximum irradiance periods, the pump
reaches its limits of flow and frequency, and so the frequency
converter stops the system. Finally, the figure
\ref{fig:ComparisonPVPS} shows the evolution of the daily productivity
of these two configurations.

<<echo=TRUE>>=
prodSP8A44Lim<-prodPVPS(lat, modeRad='prev', dataRad=prodSP8A44,
                        pump=CoefSP8A44, H=50, Pg=7000)  
as.zooY(prodSP8A44Lim)
@
%
\begin{figure}
\begin{centering}
<<echo=FALSE, fig=TRUE>>=
p=xyplot(Q~Gef|month, data=prodSP8A44Lim, 
  cex=0.5, type=c('p', 'smooth'), 
  col.symbol='gray', col.line='black')
print(p)
@
\par\end{centering}

\caption{Water flow versus irradiance of a PVPS system with a SP8A44 pump and a generator of 
  $\SI{7000}{\wattpeak}$ with a manometric height of 50 meters.\label{fig:CaudalIrradianciaLimite}}

\end{figure}

\begin{figure}
  \centering
<<echo=TRUE, fig=TRUE>>=
compPVPS <- mergesolaR(prodSP8A44, prodSP8A44Lim)
print(xyplot(compPVPS, superpose=TRUE, ylab='Yf'))
@ 
  \caption{Comparison of the daily productivity of two pumping PV systems.\label{fig:ComparisonPVPS}}
\end{figure}

\chapter{Statistical analysis of PV plants}
\label{sec:AnalisisEstadistico}
In a PV plant, the individual systems are theoretically
identical and their performance along the time should be the same. Due to their
practical differences --power tolerance, dispersion losses, dust--,
the individual performance of each system will deviate from the average
behaviour. However, when a system is performing correctly, these deviations
are constrained inside a range and should not be regarded as a sign
of malfunctioning. 

If these common deviations are assumed as a random process, a statistical
analysis of the performance of the whole set of systems can identify
a faulty system as the one that departs significantly from the mean
behaviour.

The functions \texttt{analyzeData} and \texttt{Target Diagram} compare
the daily performance of each system with a reference (for example,
the median of the whole set) during a time period of N days preceding
the current day. They calculate a set of statistics of the performance
of the PV plant as a whole, and another set of the comparison with the
reference. This statistical analysis can be summarised with a
graphical tool named "Target Diagram", which plots together the root
mean square difference, the average difference and the standard
deviation of the difference. Besides, this diagram includes the sign
of the difference of the standard deviations of the system and the
reference \cite{Perpinan2009}.

The next example uses a dataset of productivity from a PV plant
composed of 22 systems (\texttt{data(prodEx)}).  It is clear that the
system no.20 is not working correctly during these periods
(horizonplot of figure \ref{fig:StatisticalAnalysis} and target
diagram of figure \ref{fig:TargetDiagram}).

<<echo=TRUE>>=
data(prodEx)
prodStat<-analyzeData(prodEx)
@

%
\begin{figure}
  \centering
<<echo=TRUE, fig=TRUE>>=
dif<-prodEx-prodStat$stat$Median;

day=as.Date('2008-8-29')
p <- horizonplot(window(dif, start=day-60, end=day),
                 origin=0, layout=c(1, 22), 
                 colorkey=TRUE, colorkey.digits=1,
                 scales=list(y=list(relation="same")))

print(p)
@ 

\caption{Horizonplot of the differences of productivity of a set of 22 PV systems.}
\label{fig:StatisticalAnalysis}
\end{figure}

%
\begin{figure}
  \centering
 
<<echo=TRUE, fig=TRUE>>=
ndays=c(5, 10, 15, 20)
palette=brewer.pal(n=length(ndays), name='Set1')
TDColor<-TargetDiagram(prodEx, end=day, ndays=ndays,
color=palette)
print(TDColor$plot)
@ 
  \caption{``Target Diagram'' of the statistical analysis of a set of 22 systems during various time periods.}
  \label{fig:TargetDiagram}
\end{figure}

Let's remember the example devoted to \texttt{mergesolaR}, with the
result displayed in the figure \ref{fig:horizonMerge}. The function
\texttt{TargetDiagram} is an alternative tool to show the behaviour of
the set of meteorological stations (figure \ref{fig:TargetMerge}).

\begin{figure}
  \centering
<<echo=TRUE, fig=TRUE, width=8>>=
TDMadrid <- TargetDiagram(YfMadrid, 
                    end=as.POSIXct('2010-12-31'), 
                    ndays=c(10, 20, 30, 40, 50, 60), cex=0.5)
print(TDMadrid$plot)
@ 
  \caption{Target Diagram of the result of the mergesolaR example}
  \label{fig:TargetMerge}
\end{figure}

\appendix
\chapter*{Changes}

\section*{solaR 0.24}
\begin{itemize}
\item The url of \url{http://www.mapa.es/siar} has been changed to
  \url{http://www.marm.es/siar}. The \texttt{readMAPA} function has
  been fixed. Moreover, this function has been renamed to
  \texttt{readSIAR}, although it is still available as \texttt{readMAPA}.
\item The arguments \texttt{prom}, \texttt{mapa}, \texttt{bd} and
  \texttt{bdI} of \texttt{calcG0}, \texttt{calcGef}, \texttt{prodGCPV}
  and \texttt{prodPVPS} have been unified in the new
  \texttt{dataRad}. Those arguments are now deprecated although still available.
\item The mode \texttt{modeRad='mapa'} is now deprecated and renamed
  to \texttt{modeRad='siar'}.
\end{itemize}

\section*{solaR 0.23}
\begin{itemize}
\item A bug in \texttt{fInclin} when \texttt{horizBright=FALSE} has been corrected. See help(calcGef).
\item Four sun-geometry methods are now available in \texttt{fSolD, fSolI, calcSol}.
\item \texttt{ThzS, AlS and AzS} are not set to NA outside sunrise-sunset period (suggested by Joerg Trentmann)
\item \texttt{fInclin} and \texttt{calcGef} gain a new argument, \texttt{HCPV} (suggested by Álvaro Segovia)
\item A bug in \texttt{fSolI} and \texttt{fTheta} when \texttt{lat=0} has been corrected (suggested by Alberto Soria)
\item \texttt{h2d} and \texttt{d2h} small functions added.
\end{itemize}

\section*{solaR 0.22}
\begin{itemize}
\item A new \texttt{mergesolaR} method has been defined for merging
  \texttt{solaR} objects.
\item The calculation of the sunset time has been improved.
\item The voltage dependency of the efficiency curve of the inverter
  is now included in \texttt{fProd} and \texttt{calcGCPV}.
\item The default values of the \texttt{module}, \texttt{generator}
  and \texttt{inverter} of both \texttt{fProd}, \texttt{calcGCPV} and
  \texttt{optimShd} is now documented.
\item The help page of optimShd now explains correctly the concept of GRR.
\item The plot method for Shade has been renamed to shadeplot.
\item The \texttt{as.data.frame} method of the \texttt{Shade} class is
  now exported.
\end{itemize}

\section*{solaR 0.21}

solaR is now able to calculate from both daily and sub-daily irradiation values. Besides,

\begin{itemize}
\item \texttt{calcSol} and \texttt{fSolI} gain a \texttt{"BTi"} argument for intradaily time bases.

\item \texttt{fCompI} gains a \texttt{"G0I"} argument for intradaily irradiation series.

\item \texttt{fCompI} gains both \texttt{"corr"} and \texttt{"f"}.

\item \texttt{calcG0, calcGef, prodGCPV} and \texttt{prodPVPS} gain a new \texttt{"bdI"} argument
  for intradaily irradiation, and the \texttt{"corr", "f"} arguments.

\item The \texttt{"bd"} (and the new \texttt{"bdI"}) argument of \texttt{"calcG0"} can be now a
  \texttt{"Meteo"} object.  The \texttt{"file"} component of this argument can be now a
  \texttt{"zoo"} object.

\item New methods (\texttt{"losses"}, \texttt{"compareLosses"} and \texttt{"compare"}) are available for \texttt{"Gef"} and \texttt{"ProdGCPV"} classes.

\item The \texttt{"corr"} argument of \texttt{"fCompD"} (and \texttt{"fCompI"}) can be now
  \texttt{"corr=none"}.

\item The correlations between the diffuse fraction and the clearness
  index are now coded outside \texttt{"fCompD"} as separate functions. 
  Several new correlations have been included, both for monthly/daily values
  and for intradaily values.

\item New small functions for \texttt{difftime} objects have been included.

\end{itemize}


\section*{solaR 0.20}

\begin{itemize}
\item The package is now almost entirely designed with S4 classes and methods.

\item The time series object are constructed with the 'zoo' package.

\item Most of the functions and arguments have been renamed in order to 
ease the understanding by international users.

\item Two new functions have been included for the statistical analysis of a PV plant composed of several systems.

\item The package dependencies have been optimized.

\item Several new small functions for date-time calculations are now available.


\end{itemize}

\backmatter

\bibliography{solaR}
\bibliographystyle{jss}

\end{document}
