#Written by Vera Lindsay, Feb 2005

#Begining of script:

wind.plot <- function(Station,FILE,type=1){
  print(type)

#Read in the data file: 
  if(type == 1){
dat<-read.table(file=FILE,header=F, na.strings="NA", 
col.names=c("stn","year","mo","day","hr","tmp","rh","wspd","wdir","ppt"))
}else if (type == 2){
dat<-read.table(file=FILE,header=F, na.strings="NA", 
col.names=c("year","mo","day","hr","wspd","wdir"))
}else if (type == 3){
  dat<-read.table(file=FILE,header=F, na.strings="NA", 
col.names=c("year","mo","day","hr","tmp","wdir","wspd"))
}else{
 stop("Unknowen station type, please specify 1 for mof,2 for ec or 3 for mwlep.")
}

#dat[1:100,1:10]
val<-split(dat,dat$day)
day1<-val$"1"
day2<-val$"2"
day3<-val$"3"
print("Done secton 1")
#Calculat means and create graph for each day and entire data set:

# Initiate loop and set variables:
g<-1
for(g in 1:4){
if (g ==1){
day<-day1
Day<-"Day_1"
}else if (g ==2){
day<-day2
Day<-"Day_2"
}else if (g ==3){
day<-day3
Day<-"Day_3"
}else if (g == 4) {
day<-dat
Day<-"All_Days"
}

T<-paste(Day,Station,sep="_")
EXT<-".jpg"
File<-paste("graphs/2.5x5/",T,EXT,sep="")
File2 <- paste("formated_stn/",T,".means",sep="")
sdat<-day$wspd
ddat<-day$wdir
hr<-day$hr

# Calculate mean wind direction and speed for each day:
# this is a copy from wind_vec-means.r
i<-1
C<-3.1415/180
if (type == 1){
C2<-(-1000/3600)
}else{
  C2<--1
}
C3<-180/3.1415
#print("Done setting variables")

#split into componant vectors:
v<-sdat*C2*cos(ddat*10*C)
u<-sdat*C2*sin(ddat*10*C)

#Calculate hourly means for v and u:
value<-split(v,hr)
vMean<-sapply(value,mean,na.rm=TRUE)

value<-split(u,hr)
uMean<-sapply(value,mean,na.rm=TRUE)

#Scalor average (speed)
s<-sdat*C2*-1
value<-split(s,hr)
sMean<-sapply(value,mean,na.rm=TRUE)


#vector average (velocity)
vVec<-sqrt((uMean*uMean)+(vMean*vMean))
cons<-(vVec/sMean)
#print("done calculating Vector ave Vol")
#Direction:
#remove calms from u and v:
w<-1
for(w in 1:length(v)){
if (!is.na(sdat[w]) &sdat[w]==0){
v[w]<-NA
u[w]<-NA
}
}
#print ("done removing 0 from U and V")
#re-calculate means for u and v:

value<-split(v,hr)
v2Mean<-sapply(value,mean, na.rm=TRUE)

value<-split(u,hr)
u2Mean<-sapply(value,mean,na.rm=TRUE)

print("Calculating mean direction d")
C3<-180/pi
i<-1
d=matrix(data=NA,nrow=length(v2Mean),ncol=1)
for(i in 1:length(v2Mean))
if(!is.na(u2Mean[i]) & !is.na(v2Mean[i]) &(u2Mean[i]>0) && (v2Mean[i]>=0)){
d[i,]<-180 +(90-C3*atan(v2Mean[i]/u2Mean[i]))
}else if(!is.na(u2Mean[i]) & !is.na(v2Mean[i]) & (u2Mean[i]<0) && (v2Mean[i]>=0)){
d[i,]<-(270+C3*atan(-v2Mean[i]/u2Mean[i]))-180
}else if(!is.na(u2Mean[i]) & !is.na(v2Mean[i]) & (u2Mean[i]<0) && (v2Mean[i]<0)){
d[i,]<-(279-(C3*atan(v2Mean[i]/u2Mean[i]))-180)
}else{
d[i,]<-180+(90+C3*atan(-v2Mean[i]/u2Mean[i]))
}

writing values out to a file:
temp<-matrix(data=NA,nrow=24,ncol=5)
i<-1
for(i in 1:24){
temp[i,]<-c(i,sMean[i],vVec[i],cons[i],d[i])
}
write.table(temp,file=File2,col.names=c("hr","scalorM","vectorM","cons","dirMean"), row.names=F)
#End of wind_vec_means.r

#create plot for wind speed, direction and consistancy:
print("starting Graphs")

#In order to export eps files change EXE to eps(line 64),
#comment the line below and dev.off and uncomment the x11 line adnd de.copy...

jpeg(file=File,pointsize=10,width=800,height=500)
#x11(width=5,height=2.5)
par(mar=c(3,3,2,3))
par(mgp=c(3,0.5,0))
par(cex=.9)
par(tcl=-0.3)
plot(sMean,ylim=c(0,4),xlim=c(0,24),ann=F,pch=4,axes=F)
points(sMean,type="l")
mgp=c(3,0.5,0)
mtext("Wind Speed (m/s) and Consistency",font=2,side=2,line=2)
mtext("Hour",font=2,side=1,line=1.5,)
mtext(T,font=2,side=3,line=1)
axis(2,las=2)
axis(1,c(1,4,8,12,16,20,24))
box()
points(cons,col=3,pch=20)
points(cons,type="l",lty=2,col=3)
par(new=T)
par(mgp=c(3,0.5,0))
plot(d,axes=F,col=4,ann=F,pch=4,ylim=c(0,360))
axis(4,c(0,90,180,270,360),las=2)
mtext("Wind Direction (degrees)",side=4,font=2,line=2)
legend(3,365,legend=c("Avg Wspd","Avg Wdir","Consistency"),lty=c(1,NA,2),pch=c(4,4,20),col=c(1,4,3
),cex=0.9)

#dev.copy2eps(file=File)
dev.off()
}
}