Difference between revisions of "Kopra"

# Datat -> N:\Huippuyksikko\Tutkimus\R19_Kopra\Data\Model_Converted


# PM2.5

# Concentration_Background
# Kaikki


pm25kaikki <- function(koko){
# 50x50 ja 10x10 ruudut & PM2.5


AlkuAika <- Sys.time()



SummaPixCi <- NULL
polku <- "N:\\Huippuyksikko\\Tutkimus\\R19_Kopra\\Data\\Model_Converted\\"
talletuspolku <- "c:\\omat\\dataa\\Concentration_Backround\\PM2.5\\"

if(koko==50){
PopulaatioData <- NULL
PopulaatioData <-read.table(paste(polku,"Population_Europe\\50km\\ciesin_center_63lon_50km_2",sep=""))
}

if(koko==10){
PopulaatioData <- NULL
PopulaatioData <-read.table(paste(polku,"Population_Europe\\10km\\ciesin_center_63lon_10km_2",sep=""))
}




for(i in 1:12){

if(koko==50){
tiedosto <- paste("Concentration_Background\\2.5\\joinfilelist_month",i,".joined.out",sep="")
}
if(koko==10){
tiedosto <- paste("Concentration_Background\\2.5\\joinfilelist_month",i,".joined.out.1",sep="")
}

PitoisuusData <- NULL
PitoisuusData  <- read.table(paste(polku, tiedosto ,sep=""))



PitoisuusMatriisi <- NULL

for(j in 4:dim(PitoisuusData)[2])
PitoisuusMatriisi <- cbind(PitoisuusMatriisi, PitoisuusData[,j])


apu <- NULL
apu <- matrix(1,dim(PitoisuusMatriisi)[2],1)



PitoisuusVektori <- PitoisuusMatriisi%*%apu



SummaPixCi[i] <- (t(PitoisuusVektori)%*%PopulaatioData[,3])
}


# BR <- 20 m^3/päivä
# BR <- 20/(24 * 60 * 60) m^3/s

# Q <- 3479.5265 Gg/vuosi
# Q <- 3479.5265 * 10^6/(365*24*60*60) kg/s


BR <- 20/(24 * 60 * 60)

# Luetaan joka maalle oma Q_i tiedostosta 
# C:\omat\dataa\MaaKohtainenQ.csv ja lasketaan yhteen.

MaaKohtainenQ <- read.csv2("C:\\omat\\dataa\\MaaKohtainenQ.csv", header=FALSE)

Q <- sum(MaaKohtainenQ[,5])



IFConcentrationBackround <- NULL

for(i in 1:12)
IFConcentrationBackround[i] <- SummaPixCi[i]/(Q * 10^6/(365*24*60*60))* BR

if(koko==50){
TalletusKohde <- paste(talletuspolku,"IFConcentration_Backround50km.txt",sep="")
}

if(koko==10){
TalletusKohde <- paste(talletuspolku,"IFConcentration_Backround10km.txt",sep="")
}

write.table(IFConcentrationBackround, file = TalletusKohde, row.names = FALSE, col.names = FALSE)

LoppuAika <- Sys.time()

LoppuAika - AlkuAika



#rm(MaaKohtainenQ,Q, SummaPixCi, apu, polku, talletuspolku, PopulaatioData, AlkuAika, LoppuAika, tiedosto, 
#PitoisuusData, PitoisuusMatriisi,  BR, IFConcentrationBackround, TalletusKohde)
}

# Datat -> N:\Huippuyksikko\Tutkimus\R19_Kopra\Data\Model_Converted


# PM2.5-10

# Concentration_Background
# Kaikki


pm2510kaikki <- function(koko){
# 50x50 ja 10x10 ruudut & PM2.5-10


AlkuAika <- Sys.time()



SummaPixCi <- NULL
polku <- "N:\\Huippuyksikko\\Tutkimus\\R19_Kopra\\Data\\Model_Converted\\"
talletuspolku <- "c:\\omat\\dataa\\Concentration_Backround\\PM2.5-10\\"

if(koko==50){
PopulaatioData <- NULL
PopulaatioData <-read.table(paste(polku,"Population_Europe\\50km\\ciesin_center_63lon_50km_2",sep=""))
}

if(koko==10){
PopulaatioData <- NULL
PopulaatioData <-read.table(paste(polku,"Population_Europe\\10km\\ciesin_center_63lon_10km_2",sep=""))
}




for(i in 1:12){

if(koko==50){
tiedosto <- paste("Concentration_Background\\2.5-10\\joinfilelist_month",i,".joined.out",sep="")
}
if(koko==10){
tiedosto <- paste("Concentration_Background\\2.5-10\\joinfilelist_month",i,".joined.out.1",sep="")
}

PitoisuusData <- NULL
PitoisuusData  <- read.table(paste(polku, tiedosto ,sep=""))



PitoisuusMatriisi <- NULL

for(j in 4:dim(PitoisuusData)[2])
PitoisuusMatriisi <- cbind(PitoisuusMatriisi, PitoisuusData[,j])


apu <- NULL
apu <- matrix(1,dim(PitoisuusMatriisi)[2],1)



PitoisuusVektori <- PitoisuusMatriisi%*%apu



SummaPixCi[i] <- (t(PitoisuusVektori)%*%PopulaatioData[,3])
}


# BR <- 20 m^3/päivä
# BR <- 20/(24 * 60 * 60) m^3/s

# Q <- 3479.5265 Gg/vuosi
# Q <- 3479.5265 * 10^6/(365*24*60*60) kg/s


BR <- 20/(24 * 60 * 60)

# Luetaan joka maalle oma Q_i tiedostosta 
# C:\omat\dataa\MaaKohtainenQ.csv ja lasketaan yhteen.

MaaKohtainenQ <- read.csv2("C:\\omat\\dataa\\MaaKohtainenQ.csv", header=FALSE)

Q <- sum(MaaKohtainenQ[,6])



IFConcentrationBackround <- NULL

for(i in 1:12)
IFConcentrationBackround[i] <- SummaPixCi[i]/(Q * 10^6/(365*24*60*60))* BR

if(koko==50){
TalletusKohde <- paste(talletuspolku,"IFConcentration_Backround50km.txt",sep="")
}

if(koko==10){
TalletusKohde <- paste(talletuspolku,"IFConcentration_Backround10km.txt",sep="")
}

write.table(IFConcentrationBackround, file = TalletusKohde, row.names = FALSE, col.names = FALSE)

LoppuAika <- Sys.time()

LoppuAika - AlkuAika



#rm(MaaKohtainenQ,Q, SummaPixCi, apu, polku, talletuspolku, PopulaatioData, AlkuAika, LoppuAika, tiedosto, 
#PitoisuusData, PitoisuusMatriisi,  BR, IFConcentrationBackround, TalletusKohde)
}

# Datat -> N:\Huippuyksikko\Tutkimus\R19_Kopra\Data\Model_Converted

# Luetaan ensin joka maalle oma Q_i tiedostosta 
# C:\omat\dataa\MaaKohtainenQ.txt

#PM2.5

MaaKohtainenQ <- NULL
MaaKohtainenQ <- read.csv2("C:\\omat\\dataa\\MaaKohtainenQ.csv", header=FALSE)

	
Vaesto <- "Eurooppa"
HilaKoko <- 10
polku <- "N:\\Huippuyksikko\\Tutkimus\\R19_Kopra\\Data\\Model_Converted\\"

if(Vaesto == "Eurooppa"){
#Vaesto = Eurooppa
PopulaatioData <- NULL
PopulaatioData<-read.table(paste(polku,"Population_Europe\\",HilaKoko,"km\\ciesin_center_63lon_", HilaKoko,"km_2",sep=""))
}
if(Vaesto == "Suomi"){
#Vaesto = Suomi
PopulaatioData <- NULL
PopulaatioData<-read.table(paste(polku,"Population_Finland\\","fin_v00_",HilaKoko,"km.txt",sep=""))
}

apuPD <- NULL
apuPD <- PopulaatioData[PopulaatioData[,3]!=0,]

if(HilaKoko == 10) paate <- "out.1"
if(HilaKoko == 50) paate <- "out"

#Maakohtainen_Concentration_Background



AlkuAika <- Sys.time()

SummaPixCi <- NULL
talletuspolku <- "c:\\omat\\dataa\\Concentration_Backround\\PM2.5\\"

apu<-NULL

for(i in 1:12){
tiedosto<-paste("Concentration_Background\\2.5\\joinfilelist_month",i,".joined.",paate,sep="")

apuPitData <- NULL
apuPitData <- read.table(paste(polku, tiedosto,sep=""))

PitoisuusData <- NULL
PitoisuusData <- apuPitData[PopulaatioData[,3]!=0,]


# Yhdistetään Saksat Saksaksi ja Venäjän alueet Venäjäksi. Sekä poistetaan LO, LA ja All

PitoisuusMatriisi <- cbind(PitoisuusData[,4:18],PitoisuusData[,19] + PitoisuusData[,20],
PitoisuusData[,21:40], PitoisuusData[,41]+ PitoisuusData[,42]+PitoisuusData[,43]+PitoisuusData[,44],
PitoisuusData[,45:49])

apu <- (t(PitoisuusMatriisi)%*%apuPD[,3])

SummaPixCi <- rbind(SummaPixCi,t(apu))
}

# BR <- 20 m^3/päivä
# BR <- 20/(24 * 60 * 60) m^3/s
# Q <- 3479.5265 Gg/vuosi
# Q <- 3479.5265 * 10^6/(365*24*60*60) kg/s

BR <- 20/(24 * 60 * 60)

maakohtainenIFConcentrationBackround <- NULL

#PM2.5 on Q taulussa sarakkeessa 5

for(i in 1:12)
maakohtainenIFConcentrationBackround <- 
rbind(maakohtainenIFConcentrationBackround,(SummaPixCi[i,]/(t(MaaKohtainenQ[,5]) * 10^6/(365*24*60*60))) * BR)

TalletusKohde <- paste(talletuspolku,"maakohtainen_IFConcentration_Backround_",HilaKoko,"km_",Vaesto,".txt",sep="")
write.table(maakohtainenIFConcentrationBackround, file = TalletusKohde, row.names = FALSE, col.names = FALSE)

LoppuAika <- Sys.time()

LoppuAika - AlkuAika



rm(MaaKohtainenQ, SummaPixCi50km, apu, polku, talletuspolku, PopulaatioData, AlkuAika, LoppuAika, tiedosto, 
PitoisuusData, PitoisuusMatriisi,  BR, maakohtainenIFConcentrationBackround50km, TalletusKohde)




#----------------------------------------------------------------------
# Piirto PM2.5 maat erikseen
Vaesto <- "Suomi"
yAkseliYlaraja <- 4*10^(-7)
IFmat10<-
read.table(paste("c:\\omat\\dataa\\concentration_backround\\PM2.5\\maakohtainen_IFConcentration_Backround_10km_",
Vaesto,".txt",sep=""))
IFmat50<-
read.table(paste("c:\\omat\\dataa\\concentration_backround\\PM2.5\\maakohtainen_IFConcentration_Backround_50km_",
Vaesto,".txt",sep=""))

maat <- read.table("c:\\omat\\dataa\\maat.txt")

postscript(file=
paste("c:\\omat\\dataa\\concentration_backround\\PM2.5\\maaterikseen_",Vaesto,".ps",sep=""),
width=7,height=15,horizontal=F)

# par(mfrow = c(2,1))

for(i in 1:42){
maa <- paste(maat[i+1,1], maat[i+1,2],sep=", ")
plot(c(1:12),IFmat10[,i],type="b", xlim=c(1,12), ylim=c(0, yAkseliYlaraja),xlab="kuukausi",
ylab="IF 10km (punainen) ja 50km (musta)", col=rgb(1,0,0), main=maa)
par(new=T)
plot(c(1:12),IFmat50[,i],type="b", xlim=c(1,12), ylim=c(0, yAkseliYlaraja),xlab="kuukausi",
ylab="IF 10km (punainen) ja 50km (musta)", col=rgb(0,0,0), main=maa)
}

dev.off()

#---------------------------------------------------------





#---------------------




pm2510 <- function(koko){
# 50x50 ja 10x10 ruudut & PM2.5-10

AlkuAika <- Sys.time()

# Luetaan ensin joka maalle oma Q_i tiedostosta 
# C:\omat\dataa\MaaKohtainenQ.txt

MaaKohtainenQ <- read.csv2("C:\\omat\\dataa\\MaaKohtainenQ.csv", header=FALSE)




SummaPixCi <- NULL
polku <- "N:\\Huippuyksikko\\Tutkimus\\R19_Kopra\\Data\\Model_Converted\\"
talletuspolku <- "c:\\omat\\dataa\\Concentration_Backround\\PM2.5-10\\"

if(koko==50){
PopulaatioData <- NULL
PopulaatioData <-read.table(paste(polku,"Population_Europe\\50km\\ciesin_center_63lon_50km_2",sep=""))
}

if(koko==10){
PopulaatioData <- NULL
PopulaatioData <-read.table(paste(polku,"Population_Europe\\10km\\ciesin_center_63lon_10km_2",sep=""))
}


apu<-NULL


for(i in 1:12){

if(koko==50){
tiedosto <- paste("Concentration_Background\\2.5-10\\joinfilelist_month",i,".joined.out",sep="")
}
if(koko==10){
tiedosto <- paste("Concentration_Background\\2.5-10\\joinfilelist_month",i,".joined.out.1",sep="")
}

PitoisuusData <- NULL
PitoisuusData  <- read.table(paste(polku, tiedosto ,sep=""))

# Yhdistetään Saksat Saksaksi ja Venäjän alueet Venäjäksi. Sekä poistetaan LO, LA ja All

PitoisuusMatriisi <- cbind(PitoisuusData[,4:18],PitoisuusData[,19] + PitoisuusData[,20],
PitoisuusData[,21:40], PitoisuusData[,41]+ PitoisuusData[,42]+PitoisuusData[,43]+PitoisuusData[,44],
PitoisuusData[,45:49])


apu <- (t(PitoisuusMatriisi)%*%PopulaatioData[,3])

SummaPixCi <- rbind(SummaPixCi,t(apu))

}


# BR <- 20 m^3/päivä
# BR <- 20/(24 * 60 * 60) m^3/s

# Q <- 3479.5265 Gg/vuosi
# Q <- 3479.5265 * 10^6/(365*24*60*60) kg/s


BR <- 20/(24 * 60 * 60)

maakohtainenIFConcentrationBackround <- NULL

for(i in 1:12)
maakohtainenIFConcentrationBackround <- rbind(maakohtainenIFConcentrationBackround,(SummaPixCi[i,]/(t(MaaKohtainenQ[,6]) * 10^6/(365*24*60*60))) * BR)

if(koko==50){
TalletusKohde <- paste(talletuspolku,"maakohtainen_IFConcentration_Backround50km.txt",sep="")
}

if(koko==10){
TalletusKohde <- paste(talletuspolku,"maakohtainen_IFConcentration_Backround10km.txt",sep="")
}

write.table(maakohtainenIFConcentrationBackround, file = TalletusKohde, row.names = FALSE, col.names = FALSE)

LoppuAika <- Sys.time()

LoppuAika - AlkuAika



#rm(MaaKohtainenQ, SummaPixCi, apu, polku, talletuspolku, PopulaatioData, AlkuAika, LoppuAika, tiedosto, 
#PitoisuusData, PitoisuusMatriisi,  BR, maakohtainenIFConcentrationBackround, TalletusKohde)
}

#---------------------

# Datat -> N:\Huippuyksikko\Tutkimus\R19_Kopra\Data\Model_Converted


#PM2.5 Erikseen DOM, OTH, TRA, AGR,...

#Concentration_Finland

AlkuAika <- Sys.time()

PaastoLahde <- "AGR"
HilaKoko <- "10"
Vaesto <- "Eurooppa"
polku <- "N:\\Huippuyksikko\\Tutkimus\\R19_Kopra\\Data\\Model_Converted\\"
talletuspolku <- paste("c:\\omat\\dataa\\Concentration_Finland\\",PaastoLahde,"\\", sep="")

Qtaulu <- NULL
Qtaulu <- read.table(paste("N:\\Huippuyksikko\\Tutkimus\\R19_Kopra\\Data\\Model_Origin\\Emission_Silam_massbudget\\",
"Suomi","PM25Q.txt",sep=""))

if(Vaesto == "Eurooppa"){
#Vaesto = Eurooppa
PopulaatioData <- NULL
PopulaatioData<-read.table(paste(polku,"Population_Europe\\",HilaKoko,"km\\ciesin_center_63lon_", HilaKoko,"km_2",sep=""))
}
if(Vaesto == "Suomi"){
#Vaesto = Suomi
PopulaatioData <- NULL
PopulaatioData<-read.table(paste(polku,"Population_Finland\\","fin_v00_",HilaKoko,"km.txt",sep=""))
}

SummaPixCi <- NULL
apu <- NULL
Pitoisuus <- NULL
#Paasto <- NULL
IFConcentration <-NULL
TalletusKohde <- NULL
Q <- NULL


#tiedosto <- "Emission_grads\\fmsyke_data_joined_original_hila"
#Paasto <- read.table(paste(polku, tiedosto, sep=""))
#AGR:ssä vain 2 saraketta
#Q <- Paasto[,14] + Paasto[,20]
#Q <- Paasto[,12] + Paasto[,18] + Paasto[,24]

tiedosto <- paste("Concentration_Finland\\", HilaKoko,"km\\pm_0.1\\joinfilelist_PM_", PaastoLahde,".joined.out", sep="")
pm01 <- read.table(paste(polku, tiedosto, sep=""))

tiedosto <- paste("Concentration_Finland\\", HilaKoko,"km\\pm_0.1_1\\joinfilelist_PM_", PaastoLahde,".joined.out", sep="")
pm011 <- read.table(paste(polku, tiedosto, sep=""))

tiedosto <- paste("Concentration_Finland\\", HilaKoko,"km\\pm_1_2.5\\joinfilelist_PM_", PaastoLahde,".joined.out", sep="")
pm125 <- read.table(paste(polku, tiedosto, sep=""))


apu <- pm01 + pm011 + pm125

Pitoisuus <- cbind(Pitoisuus,apu[,3])
for( i in 2:9) Pitoisuus <- cbind(Pitoisuus,apu[,(i + 2 + 3)])
for( i in 10:12) Pitoisuus <- cbind(Pitoisuus,apu[,(i + 2  - 8)])

for(i in 1:12){
SummaPixCi[i] <- sum(PopulaatioData[,3]*Pitoisuus[,i])
}

# BR <- 20 m^3/päivä
# BR <- 20/(24 * 60 * 60) m^3/s
BR <- 20/(24 * 60 * 60)

#Q Kiloa sekunnissa
Q <- sum(Qtaulu[PaastoLahde==Qtaulu[,1],2:4])*(365/372)/(365*24*60*60)

IFConcentration <- SummaPixCi * BR/Q

TalletusKohde <- paste(talletuspolku,Vaesto,"_",HilaKoko,"km.txt",sep="")
write.table(IFConcentration, file = TalletusKohde, row.names = FALSE, col.names = FALSE)

LoppuAika <- Sys.time()

LoppuAika - AlkuAika

#rm(apu, SummaPixCi50km, polku, talletuspolku, PopulaatioData, AlkuAika, LoppuAika, tiedosto, 
#PitoisuusData, Q, BR, IFConcentration_FinlandAGR50km, TalletusKohde,PitoisuusAGR)


#-----------------------------------
# Kuvien piirto

#PaastoLahde <- "TRA"
piirtopolku <- paste("c:\\omat\\dataa\\Concentration_Finland\\",PaastoLahde,"\\", sep="")
postscript(file=paste(piirtopolku,"IF_kuva_",Vaesto,".ps", sep=""), width=10,height=10)#,horizontal=F)
IFmat10<-read.table(paste(piirtopolku,Vaesto,"_10km.txt", sep=""))
IFmat50<-read.table(paste(piirtopolku,Vaesto,"_50km.txt", sep=""))
plot(c(1:12),IFmat10[,1],type="b", xlim=c(1,12), ylim=c(0,1.2*10^(-6)),xlab="kuukausi", ylab="IF 10km ja 50 km (50km on musta)", col=rgb(1,0,0))
par(new=T)
plot(c(1:12),IFmat50[,1],type="b", xlim=c(1,12), ylim=c(0,1.2*10^(-6)),xlab="kuukausi", ylab="IF 10km ja 50 km (50km on musta)", col=rgb(0,0,0))
# musta on 0,0,0
dev.off()

# Datat -> N:\Huippuyksikko\Tutkimus\R19_Kopra\Data\Model_Converted


#--------------


# Concentration_Point

# 50x50 ruudut

AlkuAika <- Sys.time()

SummaPixCi50km <- matrix(nrow=12,ncol=18)

polku <- "N:\\Huippuyksikko\\Tutkimus\\R19_Kopra\\Data\\Model_Converted\\"
talletuspolku <- "c:\\omat\\dataa\\concentration_point\\"

PopulaatioData <- NULL
PopulaatioData<-read.table(paste(polku,"Population_Europe\\50km\\ciesin_center_63lon_50km_2",sep=""))


for(i in 1:12){


tiedosto1<-paste("joinfilelist_",i,"_PM_0_1.joined.out",sep="")
tiedosto2<-paste("joinfilelist_",i,"_PM_0_1_1.joined.out",sep="")
tiedosto3<-paste("joinfilelist_",i,"_PM_1_2_5.joined.out",sep="")

PitoisuusData1 <- NULL
PitoisuusData1 <-read.table(paste(polku,"Concentration_Point\\",tiedosto1,sep=""))

PitoisuusData2 <- NULL
PitoisuusData2 <-read.table(paste(polku,"Concentration_Point\\",tiedosto2,sep=""))

PitoisuusData3 <- NULL
PitoisuusData3 <-read.table(paste(polku,"Concentration_Point\\",tiedosto3,sep=""))

PM2.5PitoisuusData <- NULL
PM2.5PitoisuusData <- PitoisuusData1 + PitoisuusData2 + PitoisuusData3


for(j in 1:18)
SummaPixCi50km[i,j] <- sum( PM2.5PitoisuusData[,(j+2)] * PopulaatioData[,3] )

}



# BR <- 20 m^3/päivä
# BR <- 20/(24 * 60 * 60) m^3/s

BR <- 20/(24 * 60 * 60)


# Q <- 3 * 1/(365 * 24 * 60 * 60)kg/s

Q <- 3/(365 * 24 * 60 * 60)

IFConcentrationPoint50km <- SummaPixCi50km * BR/Q



TalletusKohde <- paste(talletuspolku,"IFConcentrationPoint50km.txt",sep="")
write.table(IFConcentrationPoint50km, file = TalletusKohde, row.names = FALSE, col.names = FALSE)



LoppuAika <- Sys.time()

LoppuAika - AlkuAika


rm(SummaPixCi50km, polku, talletuspolku, PopulaatioData, AlkuAika, LoppuAika, tiedosto1, tiedosto2, tiedosto3, PitoisuusData1, PitoisuusData2, PitoisuusData3, PM2.5PitoisuusData, Q, BR, IFConcentrationPoint50km, TalletusKohde)



#-----------------------------




# 10x10 ruudut

AlkuAika <- Sys.time()

SummaPixCi10km <- matrix(nrow=12,ncol=18)

polku <- "N:\\Huippuyksikko\\Tutkimus\\R19_Kopra\\Data\\Model_Converted\\"
talletuspolku <- "c:\\omat\\dataa\\concentration_point\\"

PopulaatioData <- NULL
PopulaatioData<-read.table(paste(polku,"Population_Europe\\10km\\ciesin_center_63lon_10km_2",sep=""))


for(i in 1:12){


tiedosto1<-paste("joinfilelist_",i,"_PM_0_1.joined.out.1",sep="")
tiedosto2<-paste("joinfilelist_",i,"_PM_0_1_1.joined.out.1",sep="")
tiedosto3<-paste("joinfilelist_",i,"_PM_1_2_5.joined.out.1",sep="")

PitoisuusData1 <- NULL
PitoisuusData1 <-read.table(paste(polku,"Concentration_Point\\",tiedosto1,sep=""))

PitoisuusData2 <- NULL
PitoisuusData2 <-read.table(paste(polku,"Concentration_Point\\",tiedosto2,sep=""))

PitoisuusData3 <- NULL
PitoisuusData3 <-read.table(paste(polku,"Concentration_Point\\",tiedosto3,sep=""))

PM2.5PitoisuusData <- NULL
PM2.5PitoisuusData <- PitoisuusData1 + PitoisuusData2 + PitoisuusData3


for(j in 1:18)
SummaPixCi10km[i,j] <- sum( PM2.5PitoisuusData[,(j+2)] * PopulaatioData[,3] )


}



# BR <- 20 m^3/päivä
# BR <- 20/(24 * 60 * 60) m^3/s

BR <- 20/(24 * 60 * 60)


# Q <- 3 * 1/(365 * 24 * 60 * 60)kg/s

Q <- 3/(365 * 24 * 60 * 60)

IFConcentrationPoint10km <- SummaPixCi10km * BR/Q



TalletusKohde <- paste(talletuspolku,"IFConcentrationPoint10km.txt",sep="")
write.table(IFConcentrationPoint10km, file = TalletusKohde, row.names = FALSE, col.names = FALSE)



LoppuAika <- Sys.time()

LoppuAika - AlkuAika


rm(SummaPixCi10km, polku, talletuspolku, PopulaatioData, AlkuAika, LoppuAika, tiedosto1, tiedosto2, tiedosto3, PitoisuusData1, PitoisuusData2, PitoisuusData3, PM2.5PitoisuusData, Q, BR, IFConcentrationPoint10km, TalletusKohde)

#-----------------------------------------------------------------------
#Maan pinnalla sijaitsevien pisteiden etäisyyden laskeva funktio
#-----------------------------------------------------------------------

#maan säde n. 6371 km 

Etaisyys <- function(sade, p1LO, p1LA, p2LO, p2LA)
{
piste1LA <- p1LA * (2*pi/360)
piste1LO <- p1LO * (2*pi/360)
piste2LA <- p2LA * (2*pi/360)
piste2LO <- p2LO * (2*pi/360)
#pallokoorinaatistoa käyttäen
SuoraValiMatka <- sade * sqrt(  (cos(piste1LO)*cos(piste1LA) - cos(piste2LO)*cos(piste2LA))^2
                              + (sin(piste1LO)*cos(piste1LA) - sin(piste2LO)*cos(piste2LA))^2
                              + (sin(piste1LA) - sin(piste2LA))^2 )
#cosini lause
Kulma <- acos(1 - SuoraValiMatka^2/(2*sade^2))
#kaaren pituus on kulma(radiaaneina)*säde
sade * Kulma
}



#----------------------------------------------------------------------------
# PopulaatioData siirtävä funktio
#-------------------------------------------------------------------------

LohkoPituus <- 354
Lohkoja <-  dim(PopulaatioData)[1]/LohkoPituus
# PD <- PDSiirto(0,0,2,0,PopulaatioData, Lohkoja, LohkoPituus)

PDSiirto <- function(poh, ita, ete, lan, PopulaatioData, Lohkoja, LohkoPituus){
UusiPD <- NULL
UusiPD <- PopulaatioData

#Siirto etelään
if(ete > 0){
#Siirretään dataa
for( i in 1:(Lohkoja-ete)){
UusiPD[((i-1)*LohkoPituus + 1):(i*LohkoPituus),3] <- 
PopulaatioData[((i + ete - 1)*LohkoPituus + 1):((i+ete)*LohkoPituus),3]
}
#Täytetään lopusta nollilla
for(i in 1:ete){
UusiPD[((Lohkoja-ete)*LohkoPituus+1):((Lohkoja - ete +1)*LohkoPituus),3] <- t(t(rep(0,LohkoPituus)))
}
}

if(lan > 0){
#Siirto länteen
#Siirretään dataa
for( i in 1:Lohkoja){
apu1 <- NULL
apu1 <- PopulaatioData[((i-1)*LohkoPituus + 1):(i*LohkoPituus),3]
UusiPD[((i-1)*LohkoPituus + 1):(i*LohkoPituus),3] <- 
t(t(c(apu1[(1+lan):LohkoPituus],rep(0,lan))))
}
}

if(poh > 0){
#Siirto pohjoiseen
#Siirretään dataa
for( i in (poh + 1):(Lohkoja)){
UusiPD[((i-1)*LohkoPituus + 1):(i*LohkoPituus),3] <- 
PopulaatioData[((i-poh)*LohkoPituus + 1):((i-poh+1)*LohkoPituus),3]
}
#Täytetään lopusta nollilla
for(i in 1:poh){
UusiPD[((i-1)*LohkoPituus+1):((i)*LohkoPituus),3] <- t(t(rep(0,LohkoPituus)))
}
}

if(ita > 0){
#Siirto itään
#Siirretään dataa
for( i in 1:Lohkoja){
apu1 <- NULL
apu1 <- PopulaatioData[((i-1)*LohkoPituus + 1):(i*LohkoPituus),3]
UusiPD[((i-1)*LohkoPituus + 1):(i*LohkoPituus),3] <-
t(t(c(rep(0,ita),apu1[1:(LohkoPituus - ita)])))
}
}

UusiPD
}


#-----------------------------------------------------------
#Lähimpien Ruutujen etsintä
#-----------------------------------------------------------
# Annetaan pisteen Lo ja La koordinaatit, piste itse, Populaatio data ja Pienhiukkasdata. 
# Viimeisin annetaan, jotta ruutujen määrää voidaan vähentää. 
#-----------------------------------------------------------
LahimmatRuudut <- function(LO,LA, PD, PM, Piste){
LOLisays <- 0.1981
LOLisaysArvio <- 0.199 
LALisays <- 0.08995
LALisaysArvio <- 0.09

apuPD <- NULL
apuPD <- PD[PD[,3]!=0 & (PM[,(1 + 2*Piste)]!=0 | PM[,(2 + 2*Piste)]!=0),]

apu1 <- NULL

#Sade 50 km
for(k in 1:dim(apuPD)[1]){
if( (apuPD[k,1] >= LO   &  apuPD[k,2] >= LA) | (apuPD[k,1] < LO   &  apuPD[k,2] < LA ) ){
apu1[k] <- 
(Etaisyys(6371,LO, LA, (apuPD[k,1] + LOLisays), (apuPD[k,2])) <= 50 &
Etaisyys(6371,LO, LA, (apuPD[k,1]), (apuPD[k,2] + LALisays)) <= 50)
}

if((apuPD[k,1] < LO   &  apuPD[k,2] >= LA) | (apuPD[k,1] >= LO   &  apuPD[k,2] < LA ) ){
apu1[k] <- 
(Etaisyys(6371,LO, LA, (apuPD[k,1]), (apuPD[k,2])) <= 50 &
Etaisyys(6371,LO, LA, (apuPD[k,1] + LOLisays), (apuPD[k,2] + LALisays)) <= 50)
}
}
apuPD <- apuPD[apu1,]


apu2 <- NULL
apu2 <- rep(FALSE,(dim(apuPD)[1])) 
for(k in 1:dim(apuPD)[1]){
if( (apuPD[k,1] <= LO   &  apuPD[k,2] <= LA) & 
(LO < (apuPD[k,1] + LOLisays) &  LA < (apuPD[k,2] + LALisays)) )
{
apu2[k] <- (1==1)
Ruutu <- k
}
}

for(k in 1:dim(apuPD)[1]){
if
(
 ((apuPD[Ruutu,1]-LOLisaysArvio)<= apuPD[k,1]) & (apuPD[k,1] <= (apuPD[Ruutu,1]+LOLisaysArvio)) &
 ((apuPD[Ruutu,2]-LALisaysArvio) <= apuPD[k,2]) & (apuPD[k,2] <= (apuPD[Ruutu,2]+LALisaysArvio))
)
{
apu2[k] <- (1==1)
if(length(apu2[apu2]) == 9){ k <- dim(apuPD)[1] + 1}
}
}
apu2
}



#------------------------------------------------------
# IF-datan piirtävä koodi
#------------------------------------------------
#polku <- "N:\\Huippuyksikko\\Tutkimus\\R19_Kopra\\Mallit\\IntakeFraction\\Distance\\"
polku <- "c:\\omat\\dataa\\concentration_point\\"
talletuspolku <- "N:\\Huippuyksikko\\Tutkimus\\R19_Kopra\\Mallit\\IntakeFraction\\Distance\\"
SadeLisays <- 50
PieninSade <- 50
MaxSade <- 1500

for(i in 1:1){	
Piste <- i 
#data <- paste(polku,"IF_Piste_",Piste,"_Distance.txt",sep="")
data <- paste(polku,"IF_Piste_",Piste,"_Distance_Siirretty_kk3.txt",sep="")
#data <- paste(polku,"IF_Piste_",Piste,"_Distance_alle_50.txt",sep="")

#kohde <- paste(talletuspolku, "Piste_",Piste,".ps",sep="")
kohde <- paste(talletuspolku, "Piste_",Piste,"_Siirretty_kk3.ps",sep="")
#kohde <- paste(talletuspolku, "Piste_",Piste,"_alle_50.ps",sep="")

postscript(file=kohde, width=10,height=10,horizontal=F)

IFDist<-read.table(data)

xAkseli <- seq(PieninSade, MaxSade, by=SadeLisays)

for(kk in 1:1){
yAkseliHI <- NULL
yAkseliLO <- NULL

for(y in 1:(dim(IFDist)[2]/2)){
yAkseliHI[y] <- IFDist[kk,(2*y - 1)]
yAkseliLO[y] <- IFDist[kk,(2*y)]
}

plot(xAkseli,yAkseliHI, xlim=c(1,(MaxSade+ 2* SadeLisays)), ylim=c(0, 2.5*10^(-6)),
xlab="Säde", ylab="IF (high = punainen, low = musta)", col=rgb(1,0,0), main = paste("Piste ", Piste, " kuukausi ", kk, sep=""))
par(new=T)
plot(xAkseli,yAkseliLO, xlim=c(1,(MaxSade+ 2* SadeLisays)), ylim=c(0, 2.5*10^(-6)),
xlab="Säde", ylab="IF (high = punainen, low = musta)", col=rgb(0,0,0), main = paste("Piste ", Piste, " kuukausi ", kk, sep=""))
}

# musta on 0,0,0
dev.off()

}




#-------------------------------------------------------------------------------------------------------------
#Perustaso

DataRuudut <- 10
# HUOM! tiedostoissa pääte .out on 50km ja .out.1 on 10km

if (DataRuudut == 10) paate <- ".out.1"
if (DataRuudut == 50) paate <- ".out"

SadeLisays <- 50
PieninSade <- 50
MaxSade <- 50


LOLisays <- 0.1981
LOLisaysArvio <- 0.199 
LALisays <- 0.08995
LALisaysArvio <- 0.09

PisteKK <- c(25,60,22,60.5,22,63,26,62,30,63,24.5,65,29.5,65,24,67.5,28,69)

polku <- "N:\\Huippuyksikko\\Tutkimus\\R19_Kopra\\Data\\Model_Converted\\"
talletuspolku <- "N:\\Huippuyksikko\\Tutkimus\\R19_Kopra\\Mallit\\IntakeFraction\\Distance\\"


PopulaatioData <- NULL
PopulaatioData <- 
read.table(paste(polku,"Population_Europe\\",DataRuudut,"km\\ciesin_center_63lon_",DataRuudut,"km_2",sep=""))




for( Piste in 1:1){
LO <- PisteKK[(2*Piste-1)]
LA <- PisteKK[(2*Piste)]

SummaPixCi <- NULL

# Seuraavassa loopissa kk on kuukausi
for(kk in 1:12){
tiedosto1<-paste("joinfilelist_",kk,"_PM_0_1.joined", paate,sep="")
tiedosto2<-paste("joinfilelist_",kk,"_PM_0_1_1.joined", paate,sep="")
tiedosto3<-paste("joinfilelist_",kk,"_PM_1_2_5.joined", paate,sep="")

PitoisuusData1 <- NULL
PitoisuusData1 <-read.table(paste(polku,"Concentration_Point\\",tiedosto1,sep=""))

PitoisuusData2 <- NULL
PitoisuusData2 <-read.table(paste(polku,"Concentration_Point\\",tiedosto2,sep=""))

PitoisuusData3 <- NULL
PitoisuusData3 <-read.table(paste(polku,"Concentration_Point\\",tiedosto3,sep=""))

PM2.5PitoisuusData<- NULL
PM2.5PitoisuusData <- PitoisuusData1 + PitoisuusData2 + PitoisuusData3

apuPM <- NULL
apuPM <- PM2.5PitoisuusData[(PM2.5PitoisuusData[,(1 + 2*Piste)]!=0 | PM2.5PitoisuusData[,(2 + 2*Piste)]!=0)
& PopulaatioData[,3]!=0,]

apuPD <- NULL
apuPD <- PopulaatioData[PopulaatioData[,3]!=0 &
(PM2.5PitoisuusData[,(1 + 2*Piste)]!=0 | PM2.5PitoisuusData[,(2 + 2*Piste)]!=0),]


apuHI <- NULL
apuLO <- NULL
apuSumma <- NULL


#apu1 <- rep(FALSE, dim(apuPD)[1])


for(Sade in seq(MaxSade, PieninSade, by= -SadeLisays )){
apu1 <- NULL

for(k in 1:dim(apuPD)[1]){
if( (apuPD[k,1] >= LO   &  apuPD[k,2] >= LA) | (apuPD[k,1] < LO   &  apuPD[k,2] < LA ) ){
apu1[k] <- 
(Etaisyys(6371,LO, LA, (apuPD[k,1] + LOLisays), (apuPD[k,2])) < Sade &
Etaisyys(6371,LO, LA, (apuPD[k,1]), (apuPD[k,2] + LALisays)) < Sade)
}

if((apuPD[k,1] < LO   &  apuPD[k,2] >= LA) | (apuPD[k,1] >= LO   &  apuPD[k,2] < LA ) ){
apu1[k] <- 
(Etaisyys(6371,LO, LA, (apuPD[k,1]), (apuPD[k,2])) < Sade &
Etaisyys(6371,LO, LA, (apuPD[k,1] + LOLisays), (apuPD[k,2] + LALisays)) < Sade)
}

}
apuPM <- apuPM[apu1,]
apuPD <- apuPD[apu1,]
}


#-------------------------------------------------------
#Viimeisten ruutujen etsintä
apu2 <- NULL
apu2 <- LahimmatRuudut(LO,LA, apuPD, apuPM, Piste)


apuPM <- apuPM[apu2,]
apuPD <- apuPD[apu2,]
#------------------------------------------------------------------

apuHI[1] <- sum( apuPM[,(1 + 2*Piste)] * apuPD[,3] )
apuLO[1] <- sum( apuPM[,(2 + 2*Piste)] * apuPD[,3] )

apuSumma[1] <- apuHI[1]
apuSumma[2] <- apuLO[1]


SummaPixCi <- rbind(SummaPixCi, t(apuSumma))
}

# BR <- 20 m^3/päivä
# BR <- 20/(24 * 60 * 60) m^3/s
BR <- 20/(24 * 60 * 60)

# Q <- 3 * 1/(365 * 24 * 60 * 60)kg/s
Q <- 3/(365 * 24 * 60 * 60)

IF <- SummaPixCi * BR/Q

TalletusKohde <- paste(talletuspolku,"IF_Piste_",Piste,"_Perus_Distance.txt",sep="")
write.table(IF, file = TalletusKohde, row.names = FALSE, col.names = FALSE)
}

#-----------------------------------------------------------------------
#Maan pinnalla sijaitsevien pisteiden etäisyyden laskeva funktio
#-----------------------------------------------------------------------

#maan säde n. 6371 km 
#jos laskutapa on 1, niin lasketaan 2 pisteen välimatka pallon pintaa pitkin
# jos laskutapa =2, niin lasketaan pisteiden (p1LO, p1LA) ja (p2LO, p2LA) euklidinen etäisyys,
# tällöin säteellä ei ole merkitystä 

Etaisyys <- function(sade, p1LO, p1LA, p2LO, p2LA, laskutapa)
{

if(laskutapa == 1){
piste1LA <- p1LA * (2*pi/360)
piste1LO <- p1LO * (2*pi/360)
piste2LA <- p2LA * (2*pi/360)
piste2LO <- p2LO * (2*pi/360)
#pallokoorinaatistoa käyttäen
SuoraValiMatka <- sade * sqrt(  (cos(piste1LO)*cos(piste1LA) - cos(piste2LO)*cos(piste2LA))^2
                              + (sin(piste1LO)*cos(piste1LA) - sin(piste2LO)*cos(piste2LA))^2
                              + (sin(piste1LA) - sin(piste2LA))^2 )
#cosini lause
Kulma <- acos(1 - SuoraValiMatka^2/(2*sade^2))
#kaaren pituus on kulma(radiaaneina)*säde
tulos <- sade * Kulma
}
if(laskutapa == 2)
{
tulos <- sqrt((p1LO - p2LO)^2 + (p1LA - p2LA)^2)
}
tulos
}



polku <- "N:\\Huippuyksikko\\Tutkimus\\R19_Kopra\\Data\\Model_Converted\\"
talletuspolku <- "N:\\Huippuyksikko\\Tutkimus\\R19_Kopra\\Mallit\\IntakeFraction\\Distance\\"


PopulaatioData <- NULL
PopulaatioData <- 
read.table(paste(polku,"Population_Europe\\",DataRuudut,"km\\ciesin_center_63lon_",DataRuudut,"km_2",sep=""))



#----------------------------------------------------------------------------
# PopulaatioData siirtävä funktio
#-------------------------------------------------------------------------

LohkoPituus <- 354
Lohkoja <-  dim(PopulaatioData)[1]/LohkoPituus
# sPD <- PDSiirto(0,0,2,0,PopulaatioData, Lohkoja, LohkoPituus)

PDSiirto <- function(poh, ita, ete, lan, PopulaatioData, Lohkoja, LohkoPituus){
UusiPD <- NULL
UusiPD <- PopulaatioData

#Siirto etelään
if(ete > 0){
#Siirretään dataa
for( i in 1:(Lohkoja-ete)){
UusiPD[((i-1)*LohkoPituus + 1):(i*LohkoPituus),3] <- 
PopulaatioData[((i + ete - 1)*LohkoPituus + 1):((i+ete)*LohkoPituus),3]
}
#Täytetään lopusta nollilla
for(i in 1:ete){
UusiPD[((Lohkoja-ete)*LohkoPituus+1):((Lohkoja - ete +1)*LohkoPituus),3] <- t(t(rep(0,LohkoPituus)))
}
}

if(lan > 0){
#Siirto länteen
#Siirretään dataa
for( i in 1:Lohkoja){
apu1 <- NULL
apu1 <- PopulaatioData[((i-1)*LohkoPituus + 1):(i*LohkoPituus),3]
UusiPD[((i-1)*LohkoPituus + 1):(i*LohkoPituus),3] <- 
t(t(c(apu1[(1+lan):LohkoPituus],rep(0,lan))))
}
}

if(poh > 0){
#Siirto pohjoiseen
#Siirretään dataa
for( i in (poh + 1):(Lohkoja)){
UusiPD[((i-1)*LohkoPituus + 1):(i*LohkoPituus),3] <- 
PopulaatioData[((i-poh)*LohkoPituus + 1):((i-poh+1)*LohkoPituus),3]
}
#Täytetään lopusta nollilla
for(i in 1:poh){
UusiPD[((i-1)*LohkoPituus+1):((i)*LohkoPituus),3] <- t(t(rep(0,LohkoPituus)))
}
}

if(ita > 0){
#Siirto itään
#Siirretään dataa
for( i in 1:Lohkoja){
apu1 <- NULL
apu1 <- PopulaatioData[((i-1)*LohkoPituus + 1):(i*LohkoPituus),3]
UusiPD[((i-1)*LohkoPituus + 1):(i*LohkoPituus),3] <-
t(t(c(rep(0,ita),apu1[1:(LohkoPituus - ita)])))
}
}

UusiPD
}


PData <- LahimmatRuudut(25,60.08,PopulaatioData,1, 58, 65)
apuri <- PopulaatioData[58 <= PopulaatioData[,2] & PopulaatioData[,2] <= 65, ]
apuri[PData,]
PData
#-----------------------------------------------------------
#Lähimpien Ruutujen etsintä
#-----------------------------------------------------------
# Annetaan pisteen Lo ja La koordinaatit, piste itse, Populaatio data ja Pienhiukkasdata. 
# Viimeisin annetaan, jotta ruutujen määrää voidaan vähentää. 
# Voidaan myös rajata antamalla alaraja ja yläraja Latitudelle
#-----------------------------------------------------------
LahimmatRuudut <-function(LO,LA, PD, Sade,Sateella, LAalaraja, LAylaraja){
#function(LO,LA, PD, PM, Piste){
LOLisays <- 0.1981
LOLisaysArvio <- 0.199 
LALisays <- 0.08995
LALisaysArvio <- 0.09

apuPD <- NULL
#apuPD <- PD[PD[,3]!=0 & (PM[,(1 + 2*Piste)]!=0 | PM[,(2 + 2*Piste)]!=0),]
#apuPD <- PD[PD[,3]!=0 & LAalaraja <= PD[,2] & PD[,2] <= LAylaraja,]
apuPD <- PD[LAalaraja <= PD[,2] & PD[,2] <= LAylaraja,]

apu2 <- NULL

if(Sateella==1){
for(k in 1:dim(apuPD)[1]){
if( (apuPD[k,1] >= LO   &  apuPD[k,2] >= LA) | (apuPD[k,1] < LO   &  apuPD[k,2] < LA ) ){
apu2[k] <- 
(Etaisyys(6371,LO, LA, (apuPD[k,1] + LOLisays), (apuPD[k,2]),1) <= Sade &
Etaisyys(6371,LO, LA, (apuPD[k,1]), (apuPD[k,2] + LALisays),1) <= Sade)
}

if((apuPD[k,1] < LO   &  apuPD[k,2] >= LA) | (apuPD[k,1] >= LO   &  apuPD[k,2] < LA ) ){
apu2[k] <- 
(Etaisyys(6371,LO, LA, (apuPD[k,1]), (apuPD[k,2]),1) <= Sade &
Etaisyys(6371,LO, LA, (apuPD[k,1] + LOLisays), (apuPD[k,2] + LALisays),1) <= Sade)
}
}
#apuPD <- apuPD[apu2,]
}


if(Sateella==0)
{
apu2 <- rep(FALSE,(dim(apuPD)[1])) 
for(k in 1:dim(apuPD)[1]){
if( (apuPD[k,1] <= LO   &  apuPD[k,2] <= LA) & 
(LO < (apuPD[k,1] + LOLisays) &  LA < (apuPD[k,2] + LALisays)) )
{
apu2[k] <- (1==1)
Ruutu <- k
}
}

for(k in 1:dim(apuPD)[1]){
if
(
 ((apuPD[Ruutu,1]-LOLisaysArvio)<= apuPD[k,1]) & (apuPD[k,1] <= (apuPD[Ruutu,1]+LOLisaysArvio)) &
 ((apuPD[Ruutu,2]-LALisaysArvio) <= apuPD[k,2]) & (apuPD[k,2] <= (apuPD[Ruutu,2]+LALisaysArvio))
)
{
apu2[k] <- (1==1)
if(length(apu2[apu2]) == 9){ k <- dim(apuPD)[1] + 1}
}
}
}
apuPD[apu2,]
}



#------------------------------------------------------
# IF-datan piirtävä koodi
#------------------------------------------------
#polku <- "N:\\Huippuyksikko\\Tutkimus\\R19_Kopra\\Mallit\\IntakeFraction\\Distance\\"
#polku <- "c:\\omat\\dataa\\concentration_point\\"
polku <- "N:\\Huippuyksikko\\Tutkimus\\R19_Kopra\\Mallit\\IntakeFraction\\Distance\\"
talletuspolku <- "N:\\Huippuyksikko\\Tutkimus\\R19_Kopra\\Mallit\\IntakeFraction\\Distance\\"
SadeLisays <- 50
PieninSade <- 50
MaxSade <- 1700

for(i in 1:1){	
Piste <- i 
data <- paste(polku,"IF_Piste_",Piste,"_Distance.txt",sep="")
#data <- paste(polku,"IF_Piste_",Piste,"_Distance_alle_50.txt",sep="")
#data <- paste(polku,"IF_Piste_",Piste,"_Distance_Siirretty.txt",sep="")

#kohde <- paste(talletuspolku, "Piste_",Piste,".ps",sep="")
kohde <- paste(talletuspolku, "Piste_",Piste,"_Siirretty.ps",sep="")
#kohde <- paste(talletuspolku, "Piste_",Piste,"_alle_50.ps",sep="")

postscript(file=kohde, width=10,height=10,horizontal=F)

IFDist<-read.table(data)

xAkseli <- seq(PieninSade, MaxSade, by=SadeLisays)

for(kk in 1:12){
yAkseliHI <- NULL
yAkseliLO <- NULL

for(y in 1:(dim(IFDist)[2]/2)){
yAkseliHI[y] <- IFDist[kk,(2*y - 1)]
yAkseliLO[y] <- IFDist[kk,(2*y)]
}

plot(xAkseli,yAkseliHI, xlim=c(1,(MaxSade+ 2* SadeLisays)), ylim=c(0, 5*10^(-6)),
xlab="Säde", ylab="IF (high = punainen, low = musta)", col=rgb(1,0,0), main = paste("Piste ", Piste, " kuukausi ", kk, sep=""))
par(new=T)
plot(xAkseli,yAkseliLO, xlim=c(1,(MaxSade+ 2* SadeLisays)), ylim=c(0, 5*10^(-6)),
xlab="Säde", ylab="IF (high = punainen, low = musta)", col=rgb(0,0,0), main = paste("Piste ", Piste, " kuukausi ", kk, sep=""))
}

# musta on 0,0,0
dev.off()

}




#-------------------------------------------------------------------------------------------------------------


DataRuudut <- 10
# HUOM! tiedostoissa pääte .out on 50km ja .out.1 on 10km

if (DataRuudut == 10) paate <- ".out.1"
if (DataRuudut == 50) paate <- ".out"

SadeLisays <- 50
PieninSade <- 50
MaxSade <- 1700


LOLisays <- 0.1981
LOLisaysArvio <- 0.199 
LALisays <- 0.08995
LALisaysArvio <- 0.09

PisteKK <- c(25,60,22,60.5,22,63,26,62,30,63,24.5,65,29.5,65,24,67.5,28,69)

polku <- "N:\\Huippuyksikko\\Tutkimus\\R19_Kopra\\Data\\Model_Converted\\"
talletuspolku <- "N:\\Huippuyksikko\\Tutkimus\\R19_Kopra\\Mallit\\IntakeFraction\\Distance\\"


PopulaatioData <- NULL
PopulaatioData <- 
read.table(paste(polku,"Population_Europe\\",DataRuudut,"km\\ciesin_center_63lon_",DataRuudut,"km_2",sep=""))

PopulaatioData <- sPD



for( Piste in 1:1){
 LO <- PisteKK[(2*Piste-1)]
 LA <- PisteKK[(2*Piste)]

 SummaPixCi <- NULL

# Seuraavassa loopissa kk on kuukausi
 for(kk in 1:12){
  tiedosto1<-paste("joinfilelist_",kk,"_PM_0_1.joined", paate,sep="")
  tiedosto2<-paste("joinfilelist_",kk,"_PM_0_1_1.joined", paate,sep="")
  tiedosto3<-paste("joinfilelist_",kk,"_PM_1_2_5.joined", paate,sep="")

  PitoisuusData1 <- NULL
  PitoisuusData1 <-read.table(paste(polku,"Concentration_Point\\",tiedosto1,sep=""))

  PitoisuusData2 <- NULL
  PitoisuusData2 <-read.table(paste(polku,"Concentration_Point\\",tiedosto2,sep=""))

  PitoisuusData3 <- NULL
  PitoisuusData3 <-read.table(paste(polku,"Concentration_Point\\",tiedosto3,sep=""))

  PM2.5PitoisuusData<- NULL
  PM2.5PitoisuusData <- PitoisuusData1 + PitoisuusData2 + PitoisuusData3

  apuPM <- NULL
  apuPM <- PM2.5PitoisuusData[(PM2.5PitoisuusData[,(1 + 2*Piste)]!=0 | PM2.5PitoisuusData[,(2 + 2*Piste)]!=0)
  & PopulaatioData[,3]!=0,]

  apuPD <- NULL
  apuPD <- PopulaatioData[PopulaatioData[,3]!=0 &
  (PM2.5PitoisuusData[,(1 + 2*Piste)]!=0 | PM2.5PitoisuusData[,(2 + 2*Piste)]!=0),]

  apuHI <- NULL
  apuLO <- NULL
  apuSumma <- NULL

#apu1 <- rep(FALSE, dim(apuPD)[1])

  sadeindeksi <- 0

  for(Sade in seq(MaxSade, PieninSade, by= -SadeLisays )){
   apu1 <- NULL
   sadeindeksi <- sadeindeksi + 1

   for(k in 1:dim(apuPD)[1]){
    if( (apuPD[k,1] >= LO   &  apuPD[k,2] >= LA) | (apuPD[k,1] < LO   &  apuPD[k,2] < LA ) ){
     apu1[k] <- 
     (Etaisyys(6371,LO, LA, (apuPD[k,1] + LOLisays), (apuPD[k,2]),1) < Sade &
     Etaisyys(6371,LO, LA, (apuPD[k,1]), (apuPD[k,2] + LALisays),1) < Sade)
    }

    if((apuPD[k,1] < LO   &  apuPD[k,2] >= LA) | (apuPD[k,1] >= LO   &  apuPD[k,2] < LA ) ){
     apu1[k] <- 
     (Etaisyys(6371,LO, LA, (apuPD[k,1]), (apuPD[k,2]),1) < Sade &
     Etaisyys(6371,LO, LA, (apuPD[k,1] + LOLisays), (apuPD[k,2] + LALisays),1) < Sade)
    }
   }
   apuPM <- apuPM[apu1,]
   apuPD <- apuPD[apu1,]

   apuHI[sadeindeksi] <- sum( apuPM[,(1 + 2*Piste)] * apuPD[,3] )
   apuLO[sadeindeksi] <- sum( apuPM[,(2 + 2*Piste)] * apuPD[,3] )

   aputaulu <- c(Sade, dim(apuPD)[1])

   TalletusKohde <- paste(talletuspolku,"Piste_",Piste,"\\IF_Piste_",Piste,"_kk_",kk,"_Sade_",Sade,"_Distance.txt",sep="")
   write.table(t(aputaulu), file = TalletusKohde, row.names = FALSE, col.names = FALSE)
  }
 
#-------------------------------------------------------
#Viimeisten ruutujen etsintä
#apu2 <- NULL
#apu2 <- LahimmatRuudut(LO,LA, apuPD, apuPM, Piste)
#apuPM <- apuPM[apu2,]
#apuPD <- apuPD[apu2,]
#apuHI[1] <- sum( apuPM[,(1 + 2*Piste)] * apuPD[,3] )
#apuLO[1] <- sum( apuPM[,(2 + 2*Piste)] * apuPD[,3] )
#apuSumma[1] <- apuHI[1]
#apuSumma[2] <- apuLO[1]
#------------------------------------------------------------------

  ylaraja <- length(apuHI)
  for(i in 1: ylaraja){
   apuSumma[(2*i - 1)] <- apuHI[(ylaraja - i + 1)]
   apuSumma[2*i] <- apuLO[(ylaraja - i + 1)]
  }

  SummaPixCi <- rbind(SummaPixCi, t(apuSumma))
 }

# BR <- 20 m^3/päivä
# BR <- 20/(24 * 60 * 60) m^3/s
 BR <- 20/(24 * 60 * 60)

# Q <- 3 * 1/(365 * 24 * 60 * 60)kg/s
 Q <- 3/(365 * 24 * 60 * 60)

 IF <- SummaPixCi * BR/Q

 TalletusKohde <- paste(talletuspolku,"IF_Piste_",Piste,"_Distance.txt",sep="")
 write.table(IF, file = TalletusKohde, row.names = FALSE, col.names = FALSE)
}





#----------------------------------------------------------
#Kolmion pinta-ala, kun tiedetään sivujen pituudet
HeronAla <- function(a,b,c){
p <- (a+b+c)/2
sqrt(p*(p-a)*(p-b)*(p-c))
}

KantaKorkeus <- function(k, h){
1/2 * k * h
}


#--------------------------------------------------------------------------

# Voi tulla ongelmia, jos piste on kulmapisteiden välisellä janalla !

# Jos palautetun vektorin 1. koordinaatti on 1, niin tarkasteltavan
# kolmion tilavuus on nolla. Muutoin ensimmäinen koordinaatti on 1, 2. koord.
# on A pisteen x-koord, 3. koord. on A pisteen y-koord., ..., 8. koord. on 1 jos
# piste on kolmion sisällä, muutoin 8. koord on nolla. 
KolmionSisalla <- function(x1,x2,D1,D2,E1,E2,F1,F2){
xvek <- c(D1,E1,F1)
yvek <- c(D2,E2,F2)
jarjvek <- NULL

paha <- 0
#Huonot tapaukset
#pisteiden D ja E kautta kulkeva suora -> testataan ovatko pisteet samalla suoralla
if(xvek[1] != xvek[2]){
 k <- (yvek[2] - yvek[1])/(xvek[2] - xvek[1])
 y <- k*(xvek[3] - xvek[1]) + yvek[1]
 if(yvek[3] == y){ paha <- 1}
}
# Kaksi pistettä samat tai samalla vaakasuoralla tai samalla pystysuoralla
if( (xvek[1] == xvek[2] & yvek[1] == yvek[2]) |
    (xvek[1] == xvek[3] & yvek[1] == yvek[3]) |
    (xvek[2] == xvek[3] & yvek[2] == yvek[3]) |
    (xvek[1] == xvek[2] & xvek[1] == xvek[3]) |
    (yvek[1] == yvek[2] & yvek[1] == yvek[3])  
   )
{
paha <- 1
}

jarjvek[1] <- paha

if(paha == 0)
{
 apu1 <- min(xvek)

 if( length(xvek[xvek==apu1])==1)
  {
   jarjvek[2] <- xvek[xvek == apu1] 
   jarjvek[3] <- yvek[xvek == apu1]
  
   apu2 <- max(yvek[xvek != apu1])

   if( length(yvek[(xvek!= apu1 & yvek == apu2)]) ==1)
    {
   jarjvek[4] <- xvek[(xvek!= apu1 & yvek == apu2)] 
   jarjvek[5] <- yvek[(xvek!= apu1 & yvek == apu2)]

   jarjvek[6] <- xvek[(xvek!= apu1 & yvek != apu2)] 
   jarjvek[7] <- yvek[(xvek!= apu1 & yvek != apu2)]
   }

   if( length(yvek[(xvek!= apu1 & yvek == apu2)]) ==2)
    {
    apu3 <- min(xvek[xvek!= apu1])
   jarjvek[4] <- xvek[(xvek!= apu1 & xvek == apu3)] 
   jarjvek[5] <- yvek[(xvek!= apu1 & xvek == apu3)]

   jarjvek[6] <- xvek[(xvek!= apu1 & xvek != apu3)] 
   jarjvek[7] <- yvek[(xvek!= apu1 & xvek != apu3)]
   }

   #Jos piste C on janan AB yläpuolella, niin vaihdetaan B ja C
   janaT <- NULL
   janaT[1] <- (jarjvek[6] - jarjvek[2])/(jarjvek[4] - jarjvek[2]) 

   janaY <- NULL 
   janaY[1] <- jarjvek[5] * janaT[1] + (1 - janaT[1])* jarjvek[3]
   if(jarjvek[7] > janaY[1]){
    apu4 <- jarjvek
    jarjvek[4:5] <- apu4[6:7]
    jarjvek[6:7] <- apu4[4:5]
   }
  }

 if( length(xvek[xvek==apu1])==2)
  {
   apu2 <- min(yvek[xvek == apu1])

   jarjvek[2] <- xvek[(xvek == apu1 & yvek == apu2)] 
   jarjvek[3] <- yvek[(xvek == apu1 & yvek == apu2)]
  
   jarjvek[4] <- xvek[(xvek == apu1 & yvek != apu2)] 
   jarjvek[5] <- yvek[(xvek == apu1 & yvek != apu2)]

   jarjvek[6] <- xvek[xvek != apu1 ] 
   jarjvek[7] <- yvek[xvek != apu1 ]
  }

}

if(jarjvek[1]==0)
{
 #Janan koordinaatit aina järjestyksessä AB, BC, AC
 janaT <- NULL
 janaT[1] <- (x1 - jarjvek[2])/(jarjvek[4] - jarjvek[2])
 janaT[2] <- (x1 - jarjvek[4])/(jarjvek[6] - jarjvek[4])
 janaT[3] <- (x1 - jarjvek[2])/(jarjvek[6] - jarjvek[2])

 janaY <- NULL
 janaY[1] <- jarjvek[5] * janaT[1] + (1 - janaT[1])* jarjvek[3]
 janaY[2] <- jarjvek[7] * janaT[2] + (1 - janaT[2])* jarjvek[5]
 janaY[3] <- jarjvek[7] * janaT[3] + (1 - janaT[3])* jarjvek[3]


 if(jarjvek[2]==jarjvek[4])
 {
  jarjvek[8] <- 0
  if((jarjvek[2]<= x1) & (x2 <= janaY[2]) & (janaY[3]<= x2))
  {
   jarjvek[8] <- 1
  }
 }

 if(jarjvek[2]!=jarjvek[4])
 {
  if(jarjvek[4]==jarjvek[6])
  {
   jarjvek[8] <- 0
   if((x1 <= jarjvek[4]) & (x2 <= janaY[1]) & (janaY[3]<= x2))
   {
    jarjvek[8] <- 1
   }
  }
  if(jarjvek[4]!=jarjvek[6])
  {
   jarjvek[8] <- 0
   if( ((jarjvek[6] > jarjvek[4]) & (x2 <= janaY[1]) &
        (x2 <= janaY[2]) & (janaY[3] <= x2)) |
       ((jarjvek[6] < jarjvek[4]) & (x2 <= janaY[1]) &
        (janaY[2] <= x2) & (janaY[3] <= x2))
     )
   {
    jarjvek[8] <- 1
   }
  }

 }
}

jarjvek
}





#------------------------------------------------------------------------
#Pistevektoriin syötetään KolmionSisalla tuottama vektori!!!!!
# x1, x2 on ympyrän keskipiste
PisteenAlue <- function(x1,x2, pistevek){

alue <- c(0,0,0)

# C Kärkenä piste C

A <- pistevek[1:2]
B <- pistevek[3:4]
C <- pistevek[5:6]


t1 <-
((x2 - C[2])*(A[1] - C[1]) - (x1 - C[1])*(A[2] - C[2]))/
((x1 - C[1])*(B[2] - A[2]) - (x2 - C[2])*(B[1] - A[1]))

if(t1==0){
suunta <- -1
if((A[1] - C[1])!=0 & ((x1 - C[1])/(A[1] - C[1]))>=0)
suunta <- 1
if((A[1] - C[1])==0 & ((x2 - C[2])/(A[2] - C[2]))>=0)
suunta <- 1


t2 <- Etaisyys(6371, x1,x2, C[1],C[2],2)/
Etaisyys(6371, A[1],A[2], C[1],C[2],2)* suunta
}

if(t1==1){
suunta <- -1
if((B[1] - C[1])!=0 & ((x1 - C[1])/(B[1] - C[1]))>=0)
suunta <- 1
if((B[1] - C[1])==0 & ((x2 - C[2])/(B[2] - C[2]))>=0)
suunta <- 1


t2 <- Etaisyys(6371, x1,x2, C[1],C[2],2)/
Etaisyys(6371, B[1],B[2], C[1],C[2],2)*suunta
}

if(t1!=0 & t1!=1)
ifelse(x1==C[1], 
{
suunta <- -1
if(sum(C - t1*B + (1 - t1)*A)/sum(C - c(x1,x2))>=0)
suunta <- 1
t2 <- Etaisyys( 6371, x1,x2,C[1],C[2],2)/
Etaisyys(6371,C[1],C[2],t1*B[1] + (1 - t1)*A[1],t1*B[2] + (1 - t1)*A[2],2) *suunta
},
t2 <- 
(x1 - C[1])/(t1*B[1] + (1 - t1)*A[1] - C[1])
)

if(0<= t1 & t1 <= 1 & 1 < t2) 
{ 
 alue[1] <- 1
 alue[2] <- t1
 alue[3] <- t2
}
if(0< t1 & t1 < 1 & t2 < 0)
{
 alue[1] <- 4
 alue[2] <- t1
 alue[3] <- t2
}
if(alue[1] == 0){

# C Kärkenä piste A
#---------------------------

A <- pistevek[3:4]
B <- pistevek[5:6]
C <- pistevek[1:2]


t1 <-
((x2 - C[2])*(A[1] - C[1]) - (x1 - C[1])*(A[2] - C[2]))/
((x1 - C[1])*(B[2] - A[2]) - (x2 - C[2])*(B[1] - A[1]))

if(t1==0){
suunta <- -1
if((A[1] - C[1])!=0 & ((x1 - C[1])/(A[1] - C[1]))>=0)
suunta <- 1
if((A[1] - C[1])==0 & ((x2 - C[2])/(A[2] - C[2]))>=0)
suunta <- 1

t2 <- Etaisyys(6371, x1,x2, C[1],C[2],2)/
Etaisyys(6371, A[1],A[2], C[1],C[2],2)* suunta
}

if(t1==1){
suunta <- -1
if((B[1] - C[1])!=0 & ((x1 - C[1])/(B[1] - C[1]))>=0)
suunta <- 1
if((B[1] - C[1])==0 & ((x2 - C[2])/(B[2] - C[2]))>=0)
suunta <- 1


t2 <- Etaisyys(6371, x1,x2, C[1],C[2],2)/
Etaisyys(6371, B[1],B[2], C[1],C[2],2)*suunta
}

if(t1!=0 & t1!=1)
ifelse(x1==C[1], 
{
suunta <- -1
if(sum(C - t1*B + (1 - t1)*A)/sum(C - c(x1,x2))>=0)
suunta <- 1

t2 <- Etaisyys( 6371, x1,x2,C[1],C[2],2)/
Etaisyys(6371,C[1],C[2],t1*B[1] + (1 - t1)*A[1],t1*B[2] + (1 - t1)*A[2],2) *suunta
},
t2 <- 
(x1 - C[1])/(t1*B[1] + (1 - t1)*A[1] - C[1])
)

#--------------------------------------
if(0<= t1 & t1 <= 1 & 1 < t2)
{
 alue[1] <- 3
 alue[2] <- t1
 alue[3] <- t2
}
if(0< t1 & t1 < 1 & t2 < 0) 
{
 alue[1] <- 6
 alue[2] <- t1
 alue[3] <- t2
}

if(alue[1] == 0){

# C Kärkenä piste B
#-----------------------------------
A <- pistevek[5:6]
B <- pistevek[1:2]
C <- pistevek[3:4]


t1 <-
((x2 - C[2])*(A[1] - C[1]) - (x1 - C[1])*(A[2] - C[2]))/
((x1 - C[1])*(B[2] - A[2]) - (x2 - C[2])*(B[1] - A[1]))

if(t1==0){
suunta <- -1
if((A[1] - C[1])!=0 & ((x1 - C[1])/(A[1] - C[1]))>=0)
suunta <- 1
if((A[1] - C[1])==0 & ((x2 - C[2])/(A[2] - C[2]))>=0)
suunta <- 1

t2 <- Etaisyys(6371, x1,x2, C[1],C[2],2)/
Etaisyys(6371, A[1],A[2], C[1],C[2],2)* suunta
}

if(t1==1){
suunta <- -1
if((B[1] - C[1])!=0 & ((x1 - C[1])/(B[1] - C[1]))>=0)
suunta <- 1
if((B[1] - C[1])==0 & ((x2 - C[2])/(B[2] - C[2]))>=0)
suunta <- 1



t2 <- Etaisyys(6371, x1,x2, C[1],C[2],2)/
Etaisyys(6371, B[1],B[2], C[1],C[2],2)*suunta
}

if(t1!=0 & t1!=1)
ifelse(x1==C[1], 
{
suunta <- -1
if(sum(C - t1*B + (1 - t1)*A)/sum(C - c(x1,x2))>=0)
suunta <- 1

t2 <- Etaisyys( 6371, x1,x2,C[1],C[2],2)/
Etaisyys(6371,C[1],C[2],t1*B[1] + (1 - t1)*A[1],t1*B[2] + (1 - t1)*A[2],2) *suunta
},
t2 <- 
(x1 - C[1])/(t1*B[1] + (1 - t1)*A[1] - C[1])
)

#---------------------------------------------------


if(0<= t1 & t1 <= 1 & 1 < t2) 
{
 alue[1] <- 5
 alue[2] <- t1
 alue[3] <- t2
}
if(0< t1 & t1 < 1 & t2 < 0) 
{
 alue[1] <- 2 
 alue[2] <- t1
 alue[3] <- t2
}
}

}
alue
}


#-----------------------------------------------------------------------------
#Laskee ympyrän ja kolmion leikkauksen alan, kun ympyrän keskipiste on yksi 
#kolmion kulmista
#Palauttaa vektorin jonka 1. komponentti on leikkausen pinta-ala
# 2. komp. kolmion pinta-ala, 3. on näiden suhde ja 4. kolmion korkeus.

AlaYmpyranKpKolmionKulma <- function(kp1, kp2, r, kv1,kv2,kv3,kv4)
{

ala <- 0
vek <- c(kp1, kp2, r, kv1,kv2,kv3,kv4)
tilavek <- KolmionSisalla(vek[1],vek[2],vek[4],vek[5],vek[6], vek[7],vek[1],vek[2])
if(tilavek[1]==1 ){
 tulos <- c(0,0,0,0)
}

if(tilavek[1]==0){


hNollaraja <- 10^(-12)
kp <- c(kp1,kp2)
A <- tilavek[2:3]
B <- tilavek[4:5]
C <- tilavek[6:7]
if(kp[1]==A[1] & kp[2]==A[2]){
A <- B 
B <- C
}
if(kp[1]==B[1] & kp[2]==B[2]){
B <- A
A <- C
}

#oikea kolmioa
if(tilavek[1]==0)
{
 AB <- Etaisyys(6371,A[1],A[2],B[1],B[2],1)
 kpB <- Etaisyys(6371,kp[1],kp[2],B[1],B[2],1)
 kpA <- Etaisyys(6371,kp[1],kp[2],A[1],A[2],1)

 cosAlpha <- (kpA^2 + AB^2 - kpB^2)/(2*AB*kpA)
 cosBeta <- (kpB^2 + AB^2 - kpA^2)/(2*AB*kpB)
 cosGamma <- (kpA^2 + kpB^2 - AB^2)/(2*kpB*kpA)

 if( cosAlpha >= 0)
 {
  h <- sqrt(kpA^2 - kpA^2*cosAlpha^2)

#A ja B eivät kuulu ympyrään ja ympyrä ei leikkaa niiden välistä janaa. 
  if(r <= h){
   ala <- 1/2 * r^2 * acos(cosGamma)
  }
# A ja B eivät kuulu ympyrään ja ympyrä leikkaa niiden välistä janaa. 
  if(h < r & kpA > r & kpB > r){
   if(cosBeta <=0) ala <- 1/2*r^2*acos(cosGamma)
   if(cosBeta > 0)
    ala <- 1/2 * r^2*(acos(h/kpB) - acos(h/r)) +
           KantaKorkeus(2*sqrt(r^2 - h^2),h) +
           1/2 *  r^2* (acos(h/kpA) - acos(h/r))
  }
# B kuuluu ympyrään ja A ei kuulu
  if(h < r & kpB <= r & kpA > r){
   if(cosBeta > 0)
    x <- sqrt(abs(kpB^2 - h^2)) + sqrt(r^2 - h^2)

   if(cosBeta <= 0)
    x <- sqrt(r^2 - h^2) - sqrt(abs(kpB^2 - h^2))

   ala <- KantaKorkeus(x,h) + 1/2 * r^2 * (acos(h/kpA) - acos(h/r))
  }

# B ei kuulu ympyrään ja A kuuluu
  if(h < r & kpB > r & kpA <= r){
   x <- sqrt(abs(kpA^2 - h^2)) + sqrt(r^2 - h^2)
   ala <- KantaKorkeus(x,h) + 1/2 * r^2 * (acos(h/kpB) - acos(h/r))
  }

# A ja B kuuluvat ympyrään
  if(h < r & kpA <= r & kpB <= r){
   ala <- KantaKorkeus(AB, h)
  }
 tulos <- c(ala, (1/2*h*AB), (ala/(1/2*h*AB)),h)
 }

if( cosAlpha < 0)
{
 h <- sqrt(kpB^2 - kpB^2*cosBeta^2)

# A ja B eivät kuulu ympyrään ja ympyrä ei leikkaa niiden välistä janaa. 
 if(r <= h){
  ala <- 1/2 * r^2 * acos(cosGamma)
 }

# A ja B eivät kuulu ympyrään ja ympyrä leikkaa niiden välistä janaa. 
# Tämä ei ole mahdollista, koska Alpha > 90 astetta

# B kuuluu ympyrään ja A ei kuulu
# Ei mahdollista, koska kpA < kpB

# B ei kuulu ympyrään ja A kuuluu
 if(h < r & kpB > r & kpA <= r){
  x <- sqrt(r^2 - h^2) - sqrt(abs(kpA^2 - h^2))
  ala <- KantaKorkeus(x,h) + 1/2 * r^2 * (acos(h/kpB) - acos(h/r))
 }

# A ja B kuuluvat ympyrään
 if(h < r & kpA <= r & kpB <= r){
  ala <- KantaKorkeus(AB, h)
 }

 tulos <- c(ala, (1/2*h*AB), (ala/(1/2*h*AB)),h)
}
}
}

tulos
}


 
#-----------------------------------------------------------------------------
#Laskee ympyrän ja kolmion leikkauksen alan, kun ympyrän keskipiste on kolmion sisällä
AlaYmpyranKpKolmionSisalla <- function(kp1, kp2, r, kv1,kv2,kv3,kv4,kv5,kv6)
{
tilavek <- KolmionSisalla(kp1,kp2,kv1,kv2,kv3,kv4,kv5,kv6)
hNollaraja <- 10^(-6)
kp <- c(kp1,kp2)
A <- tilavek[2:3]
B <- tilavek[4:5]
C <- tilavek[6:7]
KolmionKulmat <- rbind(A,B,C,A)
#oikea kolmioa ja piste kolmion sisalla
if(tilavek[1]==0 & tilavek[8]==1)
{
AB <- Etaisyys(6371,tilavek[2],tilavek[3],tilavek[4],tilavek[5],1)
BC <- Etaisyys(6371,tilavek[4],tilavek[5],tilavek[6],tilavek[7],1)
AC <- Etaisyys(6371,tilavek[2],tilavek[3],tilavek[6],tilavek[7],1)
etA <- Etaisyys(6371, kp1, kp2, tilavek[2], tilavek[3],1)
etB <- Etaisyys(6371, kp1, kp2, tilavek[4], tilavek[5],1)
etC <- Etaisyys(6371, kp1, kp2, tilavek[6], tilavek[7],1)

etv <- c(etA, etB, etC, etA)
valit <- c(AB, BC, AC)

#h <- c(-100,-100,-100)
h <- c(1,1,1)
# Tapaukset, joissa ymp. kp on kolmion kulmapiste
if( (kp == A)[1] & (kp == A)[2]){ h[1] <- 0 
  h[3] <- 0
}
if( (kp == B)[1] & (kp == B)[2]){ h[1] <- 0 
  h[2] <- 0
}
if( (kp == C)[1] & (kp == C)[2]){ h[2] <- 0 
  h[3] <- 0
}

# Tapaukset, joissa ymp. kp on kolmion kulmapisteiden välisellä suoralla
if(A[1]!=B[1])
{
 k <- (B[2] - A[2])/(B[1] - A[1])
 y <- k* (kp[1] - A[1]) + A[2]
 if( abs(kp[2] - y) < hNollaraja) h[1] <- 0
}
if( (A[1] == B[1] & kp[1]==A[1]))
{
h[1] <- 0
}

if(B[1]!=C[1])
{
 k <- (C[2] - B[2])/(C[1] - B[1])
 y <- k* (kp[1] - B[1]) + B[2]
 if( abs(kp[2] - y) < hNollaraja) h[2] <- 0
}
if( (B[1] == C[1] & kp[1]==B[1]))
{
h[2] <- 0
}

if(C[1]!=A[1])
{
 k <- (A[2] - C[2])/(A[1] - C[1])
 y <- k* (kp[1] - C[1]) + C[2]
 if( abs(kp[2] - y) < hNollaraja) h[3] <- 0
}
if( (C[1] == A[1] & kp[1]==C[1]))
{
h[3] <- 0
}

ala <- NULL
for(i in 1:3){
if(h[i] >= hNollaraja){
ala[i] <- AlaYmpyranKpKolmionKulma(kp[1], kp[2], r, KolmionKulmat[i,1],
KolmionKulmat[i,2],KolmionKulmat[(i+1),1],KolmionKulmat[(i+1),2])[1]
}
if(h[i] < hNollaraja){
ala[i] <- 0
}
}

tulos <- NULL
tulos[1] <- sum(ala)
tulos[2] <- KantaKorkeus(valit[1],
     sqrt(valit[3]^2 - ((valit[1]^2 + valit[3]^2 - valit[2]^2)/(2*valit[1]))^2))
tulos[3] <- tulos[1]/tulos[2]
}

if(tilavek[1]==1 & length(tilavek)==1)
{
tulos <- c(0,0,0)
}
if(length(tilavek) > 1 & tilavek[8]== 0)
{
tulos <- c(-1,-1,-1)
}

tulos
}


#-------------------------------------------------------------------------------------------
#Pisteiden koordinaatit syötettävä kiertosuunnassa vastapäivään!!!!
YmpyraNelikulmioKpSisalla <- function(kpLO, kpLA, r, ALO, ALA, BLO, BLA, CLO, CLA, DLO, DLA)
{
tyovek <- c(kpLO, kpLA, r, ALO, ALA, BLO, BLA, CLO, CLA, DLO, DLA)
kp <- tyovek[1:2]
KolmionKulmat <- rbind(c(tyovek[4:5],tyovek[6:7]),
c(tyovek[6:7],tyovek[8:9]),
c(tyovek[8:9],tyovek[10:11]),
c(tyovek[10:11],tyovek[4:5]))

ala <- NULL
for(i in 1:4){
ala[i] <-  AlaYmpyranKpKolmionKulma(kp[1], kp[2], r,
KolmionKulmat[i,1], KolmionKulmat[i,2],
KolmionKulmat[i,3],KolmionKulmat[i,4])[1]
}

AB <- Etaisyys(6371, tyovek[4],tyovek[5], tyovek[6], tyovek[7],1)
BC <- Etaisyys(6371, tyovek[6], tyovek[7], tyovek[8], tyovek[9],1)
CD <- Etaisyys(6371, tyovek[8], tyovek[9], tyovek[10], tyovek[11],1)
DA <- Etaisyys(6371, tyovek[10], tyovek[11], tyovek[4], tyovek[5],1)
AC <- Etaisyys(6371, tyovek[4], tyovek[5], tyovek[8], tyovek[9],1)
tulos <- NULL
tulos[1] <- sum(ala)
tulos[2] <- HeronAla(AB, BC, AC) + HeronAla(AC, CD, DA)
tulos[3] <- tulos[1]/tulos[2]
tulos
}




#----------------------------------------------------------------------------
#kv on vektori, joka sisaltaa kolmion kulmat
KolmioYmpyranAla <- function(kp1, kp2, r, kv1,kv2,kv3,kv4,kv5,kv6)
{


ala <- NULL
tilavek <- KolmionSisalla(kp1,kp2,kv1,kv2,kv3,kv4,kv5,kv6)

#oikea kolmioa ja piste kolmion sisalla
if(tilavek[1]==0 & tilavek[8]==1)
{
ala <- AlaYmpyranKpKolmionSisalla(kp1, kp2, r, kv1,kv2,kv3,kv4,kv5,kv6)
}

#Oikea kolmio ja piste kolmion ulkopuolella
if(tilavek[1]==0 & tilavek[8]==0)
{
paikkavek <- PisteenAlue(kp1,kp2,tilavek[2:7])
A <- tilavek[2:3]
B <- tilavek[4:5]
C <- tilavek[6:7]
AB <- Etaisyys(6371,A[1],A[2],B[1],B[2],1)
BC <- Etaisyys(6371,B[1],B[2],C[1],C[2],1)
AC <- Etaisyys(6371,A[1],A[2],C[1],C[2],1)

if(paikkavek[1] == 1)
{

# kp = (1 - paikkavek[3])*C + paikkavek[3]*(paikkavek[2]*B + (1-paikkavek[2])*A)
uusiA <- (1 - (paikkavek[3] + 1))*C + (paikkavek[3] + 1)*A
uusiB <- (1 - (paikkavek[3] + 1))*C + (paikkavek[3] + 1)*B
uusiC <- C
#Pisteiden koordinaatit syötettävä kiertosuunnassa vastapäivään!!!!!!!!!!

ala[1] <- AlaYmpyranKpKolmionSisalla(kp1, kp2, r,
 uusiB[1], uusiB[2], uusiA[1], uusiA[2], uusiC[1], uusiC[2])[1]- 
YmpyraNelikulmioKpSisalla(kp1, kp2, r,
   uusiB[1], uusiB[2], uusiA[1], uusiA[2], A[1], A[2], B[1], B[2])[1]
}

if(paikkavek[1] == 2)
{
#iso kolmio - pikkukolmiot
ala[1] <- AlaYmpyranKpKolmionSisalla(kp1, kp2, r, A[1], A[2], C[1],C[2], kp1, kp2)[1]-
AlaYmpyranKpKolmionSisalla(kp1, kp2, r, B[1], B[2], C[1],C[2], kp1, kp2)[1]-
AlaYmpyranKpKolmionSisalla(kp1, kp2, r, A[1], A[2], B[1],B[2], kp1, kp2)[1]
}

if(paikkavek[1] == 3)
{
uusiA <- A
uusiB <- (1 - (paikkavek[3] + 1))*A + (paikkavek[3] + 1)*B
uusiC <- (1 - (paikkavek[3] + 1))*A + (paikkavek[3] + 1)*C
#Pisteiden koordinaatit syötettävä kiertosuunnassa vastapäivään!!!!!!!!!!
ala[1] <- AlaYmpyranKpKolmionSisalla(kp1, kp2, r,
 uusiB[1], uusiB[2], uusiA[1], uusiA[2], uusiC[1], uusiC[2])[1]- 
YmpyraNelikulmioKpSisalla(kp1, kp2, r,
   uusiB[1], uusiB[2], B[1], B[2], C[1], C[2], uusiC[1], uusiC[2])[1]
}

if(paikkavek[1] == 4)
{
#iso kolmio - pikkukolmiot
ala[1] <- AlaYmpyranKpKolmionSisalla(kp1, kp2, r, A[1], A[2], B[1],B[2], kp1, kp2)[1]-
AlaYmpyranKpKolmionSisalla(kp1, kp2, r, B[1], B[2], C[1],C[2], kp1, kp2)[1]-
AlaYmpyranKpKolmionSisalla(kp1, kp2, r, A[1], A[2], C[1],C[2], kp1, kp2)[1]
}

if(paikkavek[1] == 5)
{
uusiA <- (1 - (paikkavek[3] + 1))*B + (paikkavek[3] + 1)*A
uusiB <- B
uusiC <- (1 - (paikkavek[3] + 1))*B + (paikkavek[3] + 1)*C
#Pisteiden koordinaatit syötettävä kiertosuunnassa vastapäivään!!!!!!!!!!

ala[1] <- AlaYmpyranKpKolmionSisalla(kp1, kp2, r,
 uusiB[1], uusiB[2], uusiA[1], uusiA[2], uusiC[1], uusiC[2])[1]- 
YmpyraNelikulmioKpSisalla(kp1, kp2, r,
      uusiA[1], uusiA[2], uusiC[1], uusiC[2], C[1], C[2], A[1], A[2])[1]

}

if(paikkavek[1] == 6)
{
#iso kolmio - pikkukolmiot
ala[1] <- AlaYmpyranKpKolmionSisalla(kp1, kp2, r, C[1], C[2], B[1],B[2], kp1, kp2)[1]-
AlaYmpyranKpKolmionSisalla(kp1, kp2, r, B[1], B[2], A[1],A[2], kp1, kp2)[1]-
AlaYmpyranKpKolmionSisalla(kp1, kp2, r, A[1], A[2], C[1],C[2], kp1, kp2)[1]
}

ala[2] <- HeronAla(AB,BC,AC)
ala[3] <- ala[1]/ala[2]
}

ala
}


#----------------------------------------------------------------------------
# Syötä kulmat kiertojärjestyksessä vastapäivään!!!!!!

YmpyraNelikulmioAla <-function(kp,r,kulmat)
{
tulos <- NULL
apu1 <- KolmioYmpyranAla(kp[1], kp[2], r, kulmat[1],kulmat[2],kulmat[3],kulmat[4],kulmat[7],kulmat[8])
apu2 <- KolmioYmpyranAla(kp[1], kp[2], r, kulmat[5],kulmat[6],kulmat[3],kulmat[4],kulmat[7],kulmat[8])
tulos[1:2] <- apu1[1:2] + apu2[1:2]
tulos[3] <- tulos[1]/tulos[2]
tulos
}









pisteLO <- kojetaulu[4,1]
pisteLA <- kojetaulu[4,2]

leik <- NULL
for(sade in 0:30){
leik[sade] <- YmpyraNelikulmioAla(c(25,60),sade, c(pisteLO,pisteLA,pisteLO+LOLisays,pisteLA,
pisteLO+LOLisays,pisteLA + LALisays,pisteLO,pisteLA+ LALisays))[3]
}

KolmioYmpyranAla(25,61,1,24,59,25,61,26,59)
for(km in seq(0,20,by=0.2)){
cat(km)
cat(" ")
cat(KolmioYmpyranAla(1,1,0.5*km,5,0,2,3,-2,0))
cat("\n")
}



kp1 <- -5
kp2 <- -1500
r <- 1503
tilavek <- KolmionSisalla(kp1,kp2,10,1,-4,1,-4,8)
tilavek
PisteenAlue(kp1,kp2,tilavek[2:7])
KolmioYmpyranAla(kp1,kp2,r,10,1,-4,1,-4,8)





kp1 <- 5
kp2 <- -4
#Ongelma säteellä 9!!!!!!!!!!!!!!!!!!!!!
r <- 9

KolmioYmpyranAla(kp1,kp2,9,5,2,12,2,5,6)
tilavek <- KolmionSisalla(kp1,kp2,5,2,12,2,5,6)
paikkavek <- PisteenAlue(kp1,kp2,tilavek[2:7])

x1 <- kp1
x2 <- kp2
pistevek <- tilavek[2:7]

kp <- c(1,500)
kp <- c(-4,4)
PisteenAlue(kp[1],kp[2],KolmionSisalla(kp[1],kp[2],A[1],A[2],B[1],B[2],C[1],C[2])[2:7])
x1 <- kp[1]
x2 <- kp[2]
pistevek <- KolmionSisalla(kp[1],kp[2],A[1],A[2],B[1],B[2],C[1],C[2])[2:7]

#-----------------------------------------------------------------------
#Maan pinnalla sijaitsevien pisteiden etäisyyden laskeva funktio
#-----------------------------------------------------------------------

#maan säde n. 6371 km 
#jos laskutapa on 1, niin lasketaan 2 pisteen välimatka pallon pintaa pitkin
# jos laskutapa =2, niin lasketaan pisteiden (p1LO, p1LA) ja (p2LO, p2LA) euklidinen etäisyys,
# tällöin säteellä ei ole merkitystä 

Etaisyys <- function(sade, p1LO, p1LA, p2LO, p2LA, laskutapa)
{

if(laskutapa == 1){
piste1LA <- p1LA * (2*pi/360)
piste1LO <- p1LO * (2*pi/360)
piste2LA <- p2LA * (2*pi/360)
piste2LO <- p2LO * (2*pi/360)
#pallokoorinaatistoa käyttäen
SuoraValiMatka <- sade * sqrt(  (cos(piste1LO)*cos(piste1LA) - cos(piste2LO)*cos(piste2LA))^2
                              + (sin(piste1LO)*cos(piste1LA) - sin(piste2LO)*cos(piste2LA))^2
                              + (sin(piste1LA) - sin(piste2LA))^2 )
#cosini lause
Kulma <- acos(1 - SuoraValiMatka^2/(2*sade^2))
#kaaren pituus on kulma(radiaaneina)*säde
tulos <- sade * Kulma
}
if(laskutapa == 2)
{
tulos <- sqrt((p1LO - p2LO)^2 + (p1LA - p2LA)^2)
}
tulos
}



polku <- "N:\\Huippuyksikko\\Tutkimus\\R19_Kopra\\Data\\Model_Converted\\"
talletuspolku <- "N:\\Huippuyksikko\\Tutkimus\\R19_Kopra\\Mallit\\IntakeFraction\\Distance\\"

DataRuudut <- 10
# HUOM! tiedostoissa pääte .out on 50km ja .out.1 on 10km

if (DataRuudut == 10) paate <- ".out.1"
if (DataRuudut == 50) paate <- ".out"

PopulaatioData <- NULL
PopulaatioData <- 
read.table(paste(polku,"Population_Europe\\",DataRuudut,"km\\ciesin_center_63lon_",DataRuudut,"km_2",sep=""))



#----------------------------------------------------------------------------
# PopulaatioData siirtävä funktio
#-------------------------------------------------------------------------

LohkoPituus <- 354
Lohkoja <-  dim(PopulaatioData)[1]/LohkoPituus
# sPD <- PDSiirto(0,0,2,0,PopulaatioData, Lohkoja, LohkoPituus)

PDSiirto <- function(poh, ita, ete, lan, PopulaatioData, Lohkoja, LohkoPituus){
UusiPD <- NULL
UusiPD <- PopulaatioData

#Siirto etelään
if(ete > 0){
#Siirretään dataa
for( i in 1:(Lohkoja-ete)){
UusiPD[((i-1)*LohkoPituus + 1):(i*LohkoPituus),3] <- 
PopulaatioData[((i + ete - 1)*LohkoPituus + 1):((i+ete)*LohkoPituus),3]
}
#Täytetään lopusta nollilla
for(i in 1:ete){
UusiPD[((Lohkoja-ete)*LohkoPituus+1):((Lohkoja - ete +1)*LohkoPituus),3] <- t(t(rep(0,LohkoPituus)))
}
}

if(lan > 0){
#Siirto länteen
#Siirretään dataa
for( i in 1:Lohkoja){
apu1 <- NULL
apu1 <- PopulaatioData[((i-1)*LohkoPituus + 1):(i*LohkoPituus),3]
UusiPD[((i-1)*LohkoPituus + 1):(i*LohkoPituus),3] <- 
t(t(c(apu1[(1+lan):LohkoPituus],rep(0,lan))))
}
}

if(poh > 0){
#Siirto pohjoiseen
#Siirretään dataa
for( i in (poh + 1):(Lohkoja)){
UusiPD[((i-1)*LohkoPituus + 1):(i*LohkoPituus),3] <- 
PopulaatioData[((i-poh)*LohkoPituus + 1):((i-poh+1)*LohkoPituus),3]
}
#Täytetään lopusta nollilla
for(i in 1:poh){
UusiPD[((i-1)*LohkoPituus+1):((i)*LohkoPituus),3] <- t(t(rep(0,LohkoPituus)))
}
}

if(ita > 0){
#Siirto itään
#Siirretään dataa
for( i in 1:Lohkoja){
apu1 <- NULL
apu1 <- PopulaatioData[((i-1)*LohkoPituus + 1):(i*LohkoPituus),3]
UusiPD[((i-1)*LohkoPituus + 1):(i*LohkoPituus),3] <-
t(t(c(rep(0,ita),apu1[1:(LohkoPituus - ita)])))
}
}

UusiPD
}


kojetaulu2 <- LahimmatRuudut(1,25,60.08,sPD,60,1, 58, 65)
kojetaulu3 <- LahimmatRuudut(1,25,60,sPD,500,1, 52.2,64.8)

#Pietari lähellä pistettä 30.2, 60
PietariTot <- LahimmatRuudut(2.5,30.1,60,PopulaatioData,PopulaatioData,100,1, 0, 80, 1)
PData <- LahimmatRuudut(25,60.08,PopulaatioData,1, 58, 64.8)
apuri <- PopulaatioData[58 <= PopulaatioData[,2] & PopulaatioData[,2] <= 65, ]
apuri[PData,]
PData
#-----------------------------------------------------------
#Lähimpien Ruutujen etsintä
#-----------------------------------------------------------
# Annetaan pisteen Lo ja La koordinaatit, piste itse, Populaatio data ja Pienhiukkasdata. 
# Viimeisin annetaan, jotta ruutujen määrää voidaan vähentää. 
# Voidaan myös rajata antamalla alaraja ja yläraja Latitudelle
#-----------------------------------------------------------


LahimmatRuudut <-function(Piste,LO,LA, PD, PM, Sade,Sateella, LAalaraja, LAylaraja, totuustaulu){
LOLisays <- 0.1981
LOLisaysArvio <- 0.199 
LALisays <- 0.08995
LALisaysArvio <- 0.09
apuPD <- NULL
apuPD <- PD
#apuPD <- PD[(PD[,3]!=0 & LAalaraja <= PD[,2] & PD[,2] <= LAylaraja),]
#apuPD <- PD[(PD[,3]!=0 & (PM[,]!=0 | PM[,]!=0)),]
#apuPD <- PD[PD[,3]!=0 & (PM[,(1 + 2*Piste)]!=0 | PM[,(2 + 2*Piste)]!=0),]
#apuPD <- PD[PD[,3]!=0 & LAalaraja <= PD[,2] & PD[,2] <= LAylaraja,]
#apuPD <- PD[LAalaraja <= PD[,2] & PD[,2] <= LAylaraja,]
apu2 <- NULL

if(Sateella==1){
for(k in 1:dim(apuPD)[1]){
if( (apuPD[k,1] >= LO   &  apuPD[k,2] >= LA) | (apuPD[k,1] < LO   &  apuPD[k,2] < LA ) ){
apu2[k] <- 
(Etaisyys(6371,LO, LA, (apuPD[k,1] + LOLisays), (apuPD[k,2]),1) <= Sade &
Etaisyys(6371,LO, LA, (apuPD[k,1]), (apuPD[k,2] + LALisays),1) <= Sade)
}
if((apuPD[k,1] < LO   &  apuPD[k,2] >= LA) | (apuPD[k,1] >= LO   &  apuPD[k,2] < LA ) ){
apu2[k] <- 
(Etaisyys(6371,LO, LA, (apuPD[k,1]), (apuPD[k,2]),1) <= Sade &
Etaisyys(6371,LO, LA, (apuPD[k,1] + LOLisays), (apuPD[k,2] + LALisays),1) <= Sade)
}
}
}

if(Sateella==0)
{
apu2 <- rep(FALSE,(dim(apuPD)[1])) 
for(k in 1:dim(apuPD)[1]){
if( (apuPD[k,1] <= LO   &  apuPD[k,2] <= LA) & 
(LO < (apuPD[k,1] + LOLisays) &  LA < (apuPD[k,2] + LALisays)) )
{
apu2[k] <- (1==1)
Ruutu <- k
}
}
for(k in 1:dim(apuPD)[1]){
if
(
 ((apuPD[Ruutu,1]-LOLisaysArvio)<= apuPD[k,1]) & (apuPD[k,1] <= (apuPD[Ruutu,1]+LOLisaysArvio)) &
 ((apuPD[Ruutu,2]-LALisaysArvio) <= apuPD[k,2]) & (apuPD[k,2] <= (apuPD[Ruutu,2]+LALisaysArvio))
)
{
apu2[k] <- (1==1)
if(length(apu2[apu2]) == 9){ k <- dim(apuPD)[1] + 1}
}
}
}

if(totuustaulu == 1)
{
tulos <- apu2
}
if(totuustaulu == 0)
{
tulos <- apuPD[apu2,]
}
tulos
}





#------------------------------------------------------
# IF-datan piirtävä koodi
#------------------------------------------------
#polku <- "c:\\omat\\dataa\\concentration_point\\"
polku <- "N:\\Huippuyksikko\\Tutkimus\\R19_Kopra\\Mallit\\IntakeFraction\\Distance\\"
talletuspolku <- "N:\\Huippuyksikko\\Tutkimus\\R19_Kopra\\Mallit\\IntakeFraction\\Distance\\"
SadeLisays <- 1
PieninSade <- 0
MaxSade <- 100
PiirYlaRaja <- 3.5*10^(-6)


for(i in 4:4){	
Piste <- i 
data <- paste(polku,"IF_Piste_",Piste,"_Distance_Osuus_Kokeilu_SadeVali_1_Suurin_100.txt",sep="")
#data <- paste(polku,"IF_Piste_",Piste,"_Distance.txt",sep="")
#data <- paste(polku,"IF_Piste_",Piste,"_Distance_alle_50.txt",sep="")
#data <- paste(polku,"IF_Piste_",Piste,"_Distance_Siirretty.txt",sep="")

#kohde <- paste(talletuspolku, "Piste_",Piste,".ps",sep="")
#kohde <- paste(talletuspolku, "Piste_",Piste,"_Siirretty.ps",sep="")
#kohde <- paste(talletuspolku, "Piste_",Piste,"_Siirretty_Osuus_Kokeilu_SadeVali_1_Suurin_100.ps",sep="")
#kohde <- paste(talletuspolku, "Piste_",Piste,"_alle_50.ps",sep="")
kohde <- paste(talletuspolku, "Piste_",Piste,"_Siirretty_Osuus_Kokeilu_SadeVali_1_Suurin_100_Helsingin_Skaala.ps",sep="")

postscript(file=kohde, width=10,height=10,horizontal=F)

IFDist<-read.table(data)

xAkseli <- seq(PieninSade, MaxSade, by=SadeLisays)

for(kk in 1:12){
yAkseliHI <- NULL
yAkseliLO <- NULL

for(y in 1:(dim(IFDist)[2]/2)){
yAkseliHI[y] <- IFDist[kk,(2*y - 1)]
yAkseliLO[y] <- IFDist[kk,(2*y)]
}

plot(xAkseli,yAkseliHI, xlim=c(1,(MaxSade+ 2* SadeLisays)), ylim=c(0, PiirYlaRaja),
xlab="Säde", ylab="IF (high = punainen, low = musta)", col=rgb(1,0,0), main = paste("Piste ", Piste, " kuukausi ", kk, sep=""))
par(new=T)
plot(xAkseli,yAkseliLO, xlim=c(1,(MaxSade+ 2* SadeLisays)), ylim=c(0, PiirYlaRaja),
xlab="Säde", ylab="IF (high = punainen, low = musta)", col=rgb(0,0,0), main = paste("Piste ", Piste, " kuukausi ", kk, sep=""))
}

# musta on 0,0,0
dev.off()

}




#-------------------------------------------------------------------------------------------------------------


DataRuudut <- 10
# HUOM! tiedostoissa pääte .out on 50km ja .out.1 on 10km

if (DataRuudut == 10) paate <- ".out.1"
if (DataRuudut == 50) paate <- ".out"

SadeLisays <- 50
PieninSade <- 50
MaxSade <- 1700


LOLisays <- 0.1981
LOLisaysArvio <- 0.199 
LALisays <- 0.08995
LALisaysArvio <- 0.09

PisteKK <- c(25,60,22,60.5,22,63,26,62,30,63,24.5,65,29.5,65,24,67.5,28,69)

polku <- "N:\\Huippuyksikko\\Tutkimus\\R19_Kopra\\Data\\Model_Converted\\"
talletuspolku <- "N:\\Huippuyksikko\\Tutkimus\\R19_Kopra\\Mallit\\IntakeFraction\\Distance\\"


PopulaatioData <- NULL
PopulaatioData <- 
read.table(paste(polku,"Population_Europe\\",DataRuudut,"km\\ciesin_center_63lon_",DataRuudut,"km_2",sep=""))

PopulaatioData <- sPD



for( Piste in 1:1){
 LO <- PisteKK[(2*Piste-1)]
 LA <- PisteKK[(2*Piste)]

 SummaPixCi <- NULL

# Seuraavassa loopissa kk on kuukausi
 for(kk in 1:12){
  tiedosto1<-paste("joinfilelist_",kk,"_PM_0_1.joined", paate,sep="")
  tiedosto2<-paste("joinfilelist_",kk,"_PM_0_1_1.joined", paate,sep="")
  tiedosto3<-paste("joinfilelist_",kk,"_PM_1_2_5.joined", paate,sep="")

  PitoisuusData1 <- NULL
  PitoisuusData1 <-read.table(paste(polku,"Concentration_Point\\",tiedosto1,sep=""))

  PitoisuusData2 <- NULL
  PitoisuusData2 <-read.table(paste(polku,"Concentration_Point\\",tiedosto2,sep=""))

  PitoisuusData3 <- NULL
  PitoisuusData3 <-read.table(paste(polku,"Concentration_Point\\",tiedosto3,sep=""))

  PM2.5PitoisuusData<- NULL
  PM2.5PitoisuusData <- PitoisuusData1 + PitoisuusData2 + PitoisuusData3

  apuPM <- NULL
  apuPM <- PM2.5PitoisuusData[(PM2.5PitoisuusData[,(1 + 2*Piste)]!=0 | PM2.5PitoisuusData[,(2 + 2*Piste)]!=0)
  & PopulaatioData[,3]!=0,]

  apuPD <- NULL
  apuPD <- PopulaatioData[PopulaatioData[,3]!=0 &
  (PM2.5PitoisuusData[,(1 + 2*Piste)]!=0 | PM2.5PitoisuusData[,(2 + 2*Piste)]!=0),]

  apuHI <- NULL
  apuLO <- NULL
  apuSumma <- NULL

#apu1 <- rep(FALSE, dim(apuPD)[1])

  sadeindeksi <- 0

  for(Sade in seq(MaxSade, PieninSade, by= -SadeLisays )){
   apu1 <- NULL
   sadeindeksi <- sadeindeksi + 1

   for(k in 1:dim(apuPD)[1]){
    if( (apuPD[k,1] >= LO   &  apuPD[k,2] >= LA) | (apuPD[k,1] < LO   &  apuPD[k,2] < LA ) ){
     apu1[k] <- 
     (Etaisyys(6371,LO, LA, (apuPD[k,1] + LOLisays), (apuPD[k,2]),1) < Sade &
     Etaisyys(6371,LO, LA, (apuPD[k,1]), (apuPD[k,2] + LALisays),1) < Sade)
    }

    if((apuPD[k,1] < LO   &  apuPD[k,2] >= LA) | (apuPD[k,1] >= LO   &  apuPD[k,2] < LA ) ){
     apu1[k] <- 
     (Etaisyys(6371,LO, LA, (apuPD[k,1]), (apuPD[k,2]),1) < Sade &
     Etaisyys(6371,LO, LA, (apuPD[k,1] + LOLisays), (apuPD[k,2] + LALisays),1) < Sade)
    }
   }
   apuPM <- apuPM[apu1,]
   apuPD <- apuPD[apu1,]

   apuHI[sadeindeksi] <- sum( apuPM[,(1 + 2*Piste)] * apuPD[,3] )
   apuLO[sadeindeksi] <- sum( apuPM[,(2 + 2*Piste)] * apuPD[,3] )

   aputaulu <- c(Sade, dim(apuPD)[1])

   TalletusKohde <- paste(talletuspolku,"Piste_",Piste,"\\IF_Piste_",Piste,"_kk_",kk,"_Sade_",Sade,"_Distance.txt",sep="")
   write.table(t(aputaulu), file = TalletusKohde, row.names = FALSE, col.names = FALSE)
  }
 
  ylaraja <- length(apuHI)
  for(i in 1: ylaraja){
   apuSumma[(2*i - 1)] <- apuHI[(ylaraja - i + 1)]
   apuSumma[2*i] <- apuLO[(ylaraja - i + 1)]
  }

  SummaPixCi <- rbind(SummaPixCi, t(apuSumma))
 }

# BR <- 20 m^3/päivä
# BR <- 20/(24 * 60 * 60) m^3/s
 BR <- 20/(24 * 60 * 60)

# Q <- 3 * 1/(365 * 24 * 60 * 60)kg/s
 Q <- 3/(365 * 24 * 60 * 60)

 IF <- SummaPixCi * BR/Q

 TalletusKohde <- paste(talletuspolku,"IF_Piste_",Piste,"_Distance.txt",sep="")
 write.table(IF, file = TalletusKohde, row.names = FALSE, col.names = FALSE)
}


write.table(t(t(PietariTot)), file = paste(talletuspolku,"Pietari_100km_totuus.txt",sep=""), row.names = FALSE, col.names = FALSE)


leikkausData <- NULL
osuus <- NULL
for(i in 1:dim(kojetaulu)[1]){
pisteLO <- NULL
pisteLA <- NULL
pisteLO <- kojetaulu[i,1]
pisteLA <- kojetaulu[i,2]
osuus[i] <- YmpyraNelikulmioAla(c(25,60),150, c(pisteLO,pisteLA,pisteLO+LOLisays,pisteLA,
pisteLO+LOLisays,pisteLA + LALisays,pisteLO,pisteLA+ LALisays))[3]
}

leikkausData <- osuus*kojetaulu[,3]





leik <- NULL
 for(sade in 0:30){
 leik[sade] <- YmpyraNelikulmioAla(c(25,60),sade, c(24.9272,59.9944,24.9272+LOLisays,59.9944,
 24.9272+LOLisays,59.9944 + LALisays,24.9272,59.9944+ LALisays))[3]
 }


pisteLO <- kojetaulu[4,1]
pisteLA <- kojetaulu[4,2]

leik <- NULL
for(sade in 0:30){
leik[sade] <- YmpyraNelikulmioAla(c(25,60),sade, c(pisteLO,pisteLA,pisteLO+LOLisays,pisteLA,
pisteLO+LOLisays,pisteLA + LALisays,pisteLO,pisteLA+ LALisays))[3]
}

KolmioYmpyranA25,61,1,24,59,25,61,26,59)
for(km in seq(0,20,by=0.2)){
cat(km)
cat(" ")
cat(KolmioYmpyranAla(1,1,0.5*km,5,0,2,3,-2,0))
cat("\n")
}



kp1 <- -5
kp2 <- -1500
r <- 1503
tilavek <- KolmionSisalla(kp1,kp2,10,1,-4,1,-4,8)
tilavek
PisteenAlue(kp1,kp2,tilavek[2:7])
KolmioYmpyranAla(kp1,kp2,r,10,1,-4,1,-4,8)





kp1 <- 5
kp2 <- -4
#Ongelma säteellä 9!!!!!!!!!!!!!!!!!!!!!
r <- 9

KolmioYmpyranAla(kp1,kp2,9,5,2,12,2,5,6)
tilavek <- KolmionSisalla(kp1,kp2,5,2,12,2,5,6)
paikkavek <- PisteenAlue(kp1,kp2,tilavek[2:7])

x1 <- kp1
x2 <- kp2
pistevek <- tilavek[2:7]

kp <- c(1,500)
kp <- c(-4,4)
PisteenAlue(kp[1],kp[2],KolmionSisalla(kp[1],kp[2],A[1],A[2],B[1],B[2],C[1],C[2])[2:7])
x1 <- kp[1]
x2 <- kp[2]
pistevek <- KolmionSisalla(kp[1],kp[2],A[1],A[2],B[1],B[2],C[1],C[2])[2:7]

#-----------------------------------------------------------------------
#Maan pinnalla sijaitsevien pisteiden etäisyyden laskeva funktio
#-----------------------------------------------------------------------

# maan säteen arviona voidaan käyttää 6371 km 
# jos laskutapa on 1, niin lasketaan 2 pisteen välimatka pallon pintaa pitkin
# jos laskutapa on 2, niin lasketaan pisteiden (p1LO, p1LA) ja (p2LO, p2LA) euklidinen etäisyys,
# tällöin säteellä ei ole merkitystä 

Etaisyys <- function(sade, p1LO, p1LA, p2LO, p2LA, laskutapa)
{

if(laskutapa == 1){
piste1LA <- p1LA * (2*pi/360)
piste1LO <- p1LO * (2*pi/360)
piste2LA <- p2LA * (2*pi/360)
piste2LO <- p2LO * (2*pi/360)
#pallokoorinaatistoa käyttäen
SuoraValiMatka <- sade * sqrt(  (cos(piste1LO)*cos(piste1LA) - cos(piste2LO)*cos(piste2LA))^2
                              + (sin(piste1LO)*cos(piste1LA) - sin(piste2LO)*cos(piste2LA))^2
                              + (sin(piste1LA) - sin(piste2LA))^2 )
#cosini lause
Kulma <- acos(1 - SuoraValiMatka^2/(2*sade^2))
#kaaren pituus on kulma(radiaaneina)*säde
tulos <- sade * Kulma
}
if(laskutapa == 2)
{
tulos <- sqrt((p1LO - p2LO)^2 + (p1LA - p2LA)^2)
}
tulos
}
#----------------------------------------------------------------------------
#----------------------------------------------------------
#Kolmion pinta-ala, kun tiedetään sivujen pituudet
HeronAla <- function(a,b,c){
p <- (a+b+c)/2
sqrt(p*(p-a)*(p-b)*(p-c))
}

#Tavallisella tavalla laskettava kolmion pinta-ala
KantaKorkeus <- function(k, h){
1/2 * k * h
}

#--------------------------------------------------------------------------
#----------------------------------------------------------------------------
# Voi tulla ongelmia, jos piste on kulmapisteiden välisellä janalla !

# Jos palautetun vektorin 1. koordinaatti on 1, niin tarkasteltavan
# kolmion tilavuus on nolla. Muutoin ensimmäinen koordinaatti on 0, 2. koord.
# on A pisteen x-koord, 3. koord. on A pisteen y-koord., ..., 8. koord. on 1 jos
# piste on kolmion sisällä, muutoin 8. koord on nolla. 
KolmionSisalla <- function(x1,x2,D1,D2,E1,E2,F1,F2){
xvek <- c(D1,E1,F1)
yvek <- c(D2,E2,F2)
jarjvek <- NULL

paha <- 0
#Huonot tapaukset
#pisteiden D ja E kautta kulkeva suora -> testataan ovatko pisteet samalla suoralla
if(xvek[1] != xvek[2]){
 k <- (yvek[2] - yvek[1])/(xvek[2] - xvek[1])
 y <- k*(xvek[3] - xvek[1]) + yvek[1]
 if(yvek[3] == y){ paha <- 1}
}
# Kaksi pistettä samat tai samalla vaakasuoralla tai samalla pystysuoralla
if( (xvek[1] == xvek[2] & yvek[1] == yvek[2]) |
    (xvek[1] == xvek[3] & yvek[1] == yvek[3]) |
    (xvek[2] == xvek[3] & yvek[2] == yvek[3]) |
    (xvek[1] == xvek[2] & xvek[1] == xvek[3]) |
    (yvek[1] == yvek[2] & yvek[1] == yvek[3])  
   )
{
paha <- 1
}

jarjvek[1] <- paha

if(paha == 0)
{
 apu1 <- min(xvek)

 if( length(xvek[xvek==apu1])==1)
  {
   jarjvek[2] <- xvek[xvek == apu1] 
   jarjvek[3] <- yvek[xvek == apu1]
  
   apu2 <- max(yvek[xvek != apu1])

   if( length(yvek[(xvek!= apu1 & yvek == apu2)]) ==1)
    {
   jarjvek[4] <- xvek[(xvek!= apu1 & yvek == apu2)] 
   jarjvek[5] <- yvek[(xvek!= apu1 & yvek == apu2)]

   jarjvek[6] <- xvek[(xvek!= apu1 & yvek != apu2)] 
   jarjvek[7] <- yvek[(xvek!= apu1 & yvek != apu2)]
   }

   if( length(yvek[(xvek!= apu1 & yvek == apu2)]) ==2)
    {
    apu3 <- min(xvek[xvek!= apu1])
   jarjvek[4] <- xvek[(xvek!= apu1 & xvek == apu3)] 
   jarjvek[5] <- yvek[(xvek!= apu1 & xvek == apu3)]

   jarjvek[6] <- xvek[(xvek!= apu1 & xvek != apu3)] 
   jarjvek[7] <- yvek[(xvek!= apu1 & xvek != apu3)]
   }

   #Jos piste C on janan AB yläpuolella, niin vaihdetaan B ja C
   janaT <- NULL
   janaT[1] <- (jarjvek[6] - jarjvek[2])/(jarjvek[4] - jarjvek[2]) 

   janaY <- NULL 
   janaY[1] <- jarjvek[5] * janaT[1] + (1 - janaT[1])* jarjvek[3]
   if(jarjvek[7] > janaY[1]){
    apu4 <- jarjvek
    jarjvek[4:5] <- apu4[6:7]
    jarjvek[6:7] <- apu4[4:5]
   }
  }

 if( length(xvek[xvek==apu1])==2)
  {
   apu2 <- min(yvek[xvek == apu1])

   jarjvek[2] <- xvek[(xvek == apu1 & yvek == apu2)] 
   jarjvek[3] <- yvek[(xvek == apu1 & yvek == apu2)]
  
   jarjvek[4] <- xvek[(xvek == apu1 & yvek != apu2)] 
   jarjvek[5] <- yvek[(xvek == apu1 & yvek != apu2)]

   jarjvek[6] <- xvek[xvek != apu1 ] 
   jarjvek[7] <- yvek[xvek != apu1 ]
  }

}

if(jarjvek[1]==0)
{
 #Janan koordinaatit aina järjestyksessä AB, BC, AC
 janaT <- NULL
 janaT[1] <- (x1 - jarjvek[2])/(jarjvek[4] - jarjvek[2])
 janaT[2] <- (x1 - jarjvek[4])/(jarjvek[6] - jarjvek[4])
 janaT[3] <- (x1 - jarjvek[2])/(jarjvek[6] - jarjvek[2])

 janaY <- NULL
 janaY[1] <- jarjvek[5] * janaT[1] + (1 - janaT[1])* jarjvek[3]
 janaY[2] <- jarjvek[7] * janaT[2] + (1 - janaT[2])* jarjvek[5]
 janaY[3] <- jarjvek[7] * janaT[3] + (1 - janaT[3])* jarjvek[3]


 if(jarjvek[2]==jarjvek[4])
 {
  jarjvek[8] <- 0
  if((jarjvek[2]<= x1) & (x2 <= janaY[2]) & (janaY[3]<= x2))
  {
   jarjvek[8] <- 1
  }
 }

 if(jarjvek[2]!=jarjvek[4])
 {
  if(jarjvek[4]==jarjvek[6])
  {
   jarjvek[8] <- 0
   if((x1 <= jarjvek[4]) & (x2 <= janaY[1]) & (janaY[3]<= x2))
   {
    jarjvek[8] <- 1
   }
  }
  if(jarjvek[4]!=jarjvek[6])
  {
   jarjvek[8] <- 0
   if( ((jarjvek[6] > jarjvek[4]) & (x2 <= janaY[1]) &
        (x2 <= janaY[2]) & (janaY[3] <= x2)) |
       ((jarjvek[6] < jarjvek[4]) & (x2 <= janaY[1]) &
        (janaY[2] <= x2) & (janaY[3] <= x2))
     )
   {
    jarjvek[8] <- 1
   }
  }

 }
}

jarjvek
}

#----------------------------------------------------------------------------
#------------------------------------------------------------------------
#Pistevektoriin syötetään KolmionSisalla tuottama vektori!!!!!
# x1, x2 on ympyrän keskipiste
PisteenAlue <- function(x1,x2, pistevek){

alue <- c(0,0,0)

# C Kärkenä piste C

A <- pistevek[1:2]
B <- pistevek[3:4]
C <- pistevek[5:6]


t1 <-
((x2 - C[2])*(A[1] - C[1]) - (x1 - C[1])*(A[2] - C[2]))/
((x1 - C[1])*(B[2] - A[2]) - (x2 - C[2])*(B[1] - A[1]))

if(t1==0){
suunta <- -1
if((A[1] - C[1])!=0 & ((x1 - C[1])/(A[1] - C[1]))>=0)
suunta <- 1
if((A[1] - C[1])==0 & ((x2 - C[2])/(A[2] - C[2]))>=0)
suunta <- 1


t2 <- Etaisyys(6371, x1,x2, C[1],C[2],2)/
Etaisyys(6371, A[1],A[2], C[1],C[2],2)* suunta
}

if(t1==1){
suunta <- -1
if((B[1] - C[1])!=0 & ((x1 - C[1])/(B[1] - C[1]))>=0)
suunta <- 1
if((B[1] - C[1])==0 & ((x2 - C[2])/(B[2] - C[2]))>=0)
suunta <- 1


t2 <- Etaisyys(6371, x1,x2, C[1],C[2],2)/
Etaisyys(6371, B[1],B[2], C[1],C[2],2)*suunta
}

if(t1!=0 & t1!=1)
ifelse(x1==C[1], 
{
suunta <- -1
if(sum(C - t1*B + (1 - t1)*A)/sum(C - c(x1,x2))>=0)
suunta <- 1
t2 <- Etaisyys( 6371, x1,x2,C[1],C[2],2)/
Etaisyys(6371,C[1],C[2],t1*B[1] + (1 - t1)*A[1],t1*B[2] + (1 - t1)*A[2],2) *suunta
},
t2 <- 
(x1 - C[1])/(t1*B[1] + (1 - t1)*A[1] - C[1])
)

if(0<= t1 & t1 <= 1 & 1 < t2) 
{ 
 alue[1] <- 1
 alue[2] <- t1
 alue[3] <- t2
}
if(0< t1 & t1 < 1 & t2 < 0)
{
 alue[1] <- 4
 alue[2] <- t1
 alue[3] <- t2
}
if(alue[1] == 0){

# C Kärkenä piste A
#---------------------------

A <- pistevek[3:4]
B <- pistevek[5:6]
C <- pistevek[1:2]


t1 <-
((x2 - C[2])*(A[1] - C[1]) - (x1 - C[1])*(A[2] - C[2]))/
((x1 - C[1])*(B[2] - A[2]) - (x2 - C[2])*(B[1] - A[1]))

if(t1==0){
suunta <- -1
if((A[1] - C[1])!=0 & ((x1 - C[1])/(A[1] - C[1]))>=0)
suunta <- 1
if((A[1] - C[1])==0 & ((x2 - C[2])/(A[2] - C[2]))>=0)
suunta <- 1

t2 <- Etaisyys(6371, x1,x2, C[1],C[2],2)/
Etaisyys(6371, A[1],A[2], C[1],C[2],2)* suunta
}

if(t1==1){
suunta <- -1
if((B[1] - C[1])!=0 & ((x1 - C[1])/(B[1] - C[1]))>=0)
suunta <- 1
if((B[1] - C[1])==0 & ((x2 - C[2])/(B[2] - C[2]))>=0)
suunta <- 1


t2 <- Etaisyys(6371, x1,x2, C[1],C[2],2)/
Etaisyys(6371, B[1],B[2], C[1],C[2],2)*suunta
}

if(t1!=0 & t1!=1)
ifelse(x1==C[1], 
{
suunta <- -1
if(sum(C - t1*B + (1 - t1)*A)/sum(C - c(x1,x2))>=0)
suunta <- 1

t2 <- Etaisyys( 6371, x1,x2,C[1],C[2],2)/
Etaisyys(6371,C[1],C[2],t1*B[1] + (1 - t1)*A[1],t1*B[2] + (1 - t1)*A[2],2) *suunta
},
t2 <- 
(x1 - C[1])/(t1*B[1] + (1 - t1)*A[1] - C[1])
)

#--------------------------------------
if(0<= t1 & t1 <= 1 & 1 < t2)
{
 alue[1] <- 3
 alue[2] <- t1
 alue[3] <- t2
}
if(0< t1 & t1 < 1 & t2 < 0) 
{
 alue[1] <- 6
 alue[2] <- t1
 alue[3] <- t2
}

if(alue[1] == 0){

# C Kärkenä piste B
#-----------------------------------
A <- pistevek[5:6]
B <- pistevek[1:2]
C <- pistevek[3:4]


t1 <-
((x2 - C[2])*(A[1] - C[1]) - (x1 - C[1])*(A[2] - C[2]))/
((x1 - C[1])*(B[2] - A[2]) - (x2 - C[2])*(B[1] - A[1]))

if(t1==0){
suunta <- -1
if((A[1] - C[1])!=0 & ((x1 - C[1])/(A[1] - C[1]))>=0)
suunta <- 1
if((A[1] - C[1])==0 & ((x2 - C[2])/(A[2] - C[2]))>=0)
suunta <- 1

t2 <- Etaisyys(6371, x1,x2, C[1],C[2],2)/
Etaisyys(6371, A[1],A[2], C[1],C[2],2)* suunta
}

if(t1==1){
suunta <- -1
if((B[1] - C[1])!=0 & ((x1 - C[1])/(B[1] - C[1]))>=0)
suunta <- 1
if((B[1] - C[1])==0 & ((x2 - C[2])/(B[2] - C[2]))>=0)
suunta <- 1



t2 <- Etaisyys(6371, x1,x2, C[1],C[2],2)/
Etaisyys(6371, B[1],B[2], C[1],C[2],2)*suunta
}

if(t1!=0 & t1!=1)
ifelse(x1==C[1], 
{
suunta <- -1
if(sum(C - t1*B + (1 - t1)*A)/sum(C - c(x1,x2))>=0)
suunta <- 1

t2 <- Etaisyys( 6371, x1,x2,C[1],C[2],2)/
Etaisyys(6371,C[1],C[2],t1*B[1] + (1 - t1)*A[1],t1*B[2] + (1 - t1)*A[2],2) *suunta
},
t2 <- 
(x1 - C[1])/(t1*B[1] + (1 - t1)*A[1] - C[1])
)

#---------------------------------------------------


if(0<= t1 & t1 <= 1 & 1 < t2) 
{
 alue[1] <- 5
 alue[2] <- t1
 alue[3] <- t2
}
if(0< t1 & t1 < 1 & t2 < 0) 
{
 alue[1] <- 2 
 alue[2] <- t1
 alue[3] <- t2
}
}

}
alue
}
#----------------------------------------------------------------------------


#kp1, kp2, r, A[1],A[2],C[1],C[2]
#kulm <- c(B[1],B[2],C[1],C[2])
#kv1 <- kulm[1]
#kv2 <- kulm[2]
#kv3 <- kulm[3]
#kv4 <- kulm[4]

#----------------------------------------------------------------------------
#-----------------------------------------------------------------------------
#Laskee ympyrän ja kolmion leikkauksen alan, kun ympyrän keskipiste on yksi 
#kolmion kulmista
#Palauttaa vektorin jonka 1. komponentti on leikkausen pinta-ala
# 2. komp. kolmion pinta-ala, 3. on näiden suhde ja 4. kolmion korkeus.

AlaYmpyranKpKolmionKulma <- function(kp1, kp2, r, kv1,kv2,kv3,kv4)
{

ala <- 0
vek <- c(kp1, kp2, r, kv1,kv2,kv3,kv4)
tilavek <- KolmionSisalla(vek[1],vek[2],vek[4],vek[5],vek[6], vek[7],vek[1],vek[2])
if(tilavek[1]==1 ){
 tulos <- c(0,0,0,0)
}

# Oikea kolmio
if(tilavek[1]==0){


hNollaraja <- 10^(-12)
kp <- c(kp1,kp2)
A <- tilavek[2:3]
B <- tilavek[4:5]
C <- tilavek[6:7]
if(kp[1]==A[1] & kp[2]==A[2]){
A <- B 
B <- C
}
if(kp[1]==B[1] & kp[2]==B[2]){
B <- A
A <- C
}

 AB <- Etaisyys(6371,A[1],A[2],B[1],B[2],1)
 kpB <- Etaisyys(6371,kp[1],kp[2],B[1],B[2],1)
 kpA <- Etaisyys(6371,kp[1],kp[2],A[1],A[2],1)

 cosAlpha <- (kpA^2 + AB^2 - kpB^2)/(2*AB*kpA)
 cosBeta <- (kpB^2 + AB^2 - kpA^2)/(2*AB*kpB)
 cosGamma <- (kpA^2 + kpB^2 - AB^2)/(2*kpB*kpA)

 if( cosAlpha >= 0)
 {
  h <- sqrt(kpA^2 - kpA^2*cosAlpha^2)

#A ja B eivät kuulu ympyrään ja ympyrä ei leikkaa niiden välistä janaa. 
  if(r <= h | (h < r & r <= kpB & cosBeta < 0)){
   ala <- 1/2 * r^2 * acos(cosGamma)
  }
# A ja B eivät kuulu ympyrään ja ympyrä leikkaa niiden välistä janaa. 
  if(h < r & kpA > r & kpB > r){
   if(cosBeta <=0) ala <- 1/2*r^2*acos(cosGamma)
   if(cosBeta > 0)
    ala <- 1/2 * r^2*(acos(h/kpB) - acos(h/r)) +
           KantaKorkeus(2*sqrt(r^2 - h^2),h) +
           1/2 *  r^2* (acos(h/kpA) - acos(h/r))
  }
# B kuuluu ympyrään ja A ei kuulu
  if(h < r & kpB <= r & kpA > r){
   if(cosBeta > 0)
    x <- sqrt(abs(kpB^2 - h^2)) + sqrt(r^2 - h^2)

   if(cosBeta <= 0)
    x <- sqrt(r^2 - h^2) - sqrt(abs(kpB^2 - h^2))

   ala <- KantaKorkeus(x,h) + 1/2 * r^2 * (acos(h/kpA) - acos(h/r))
  }

# B ei kuulu ympyrään ja A kuuluu
  if(h < r & kpB > r & kpA <= r){
   x <- sqrt(abs(kpA^2 - h^2)) + sqrt(r^2 - h^2)
   ala <- KantaKorkeus(x,h) + 1/2 * r^2 * (acos(h/kpB) - acos(h/r))
  }

# A ja B kuuluvat ympyrään
  if(h < r & kpA <= r & kpB <= r){
   ala <- KantaKorkeus(AB, h)
  }
 tulos <- c(ala, (1/2*h*AB), (ala/(1/2*h*AB)),h)
 }

if( cosAlpha < 0)
{
 h <- sqrt(kpB^2 - kpB^2*cosBeta^2)

# A ja B eivät kuulu ympyrään ja ympyrä ei leikkaa niiden välistä janaa. 
 if(r <= h | (h < r & r <= kpA)){
  ala <- 1/2 * r^2 * acos(cosGamma)
 }

# A ja B eivät kuulu ympyrään ja ympyrä leikkaa niiden välistä janaa. 
# Tämä ei ole mahdollista, koska Alpha > 90 astetta

# B kuuluu ympyrään ja A ei kuulu
# Ei mahdollista, koska kpA < kpB

# B ei kuulu ympyrään ja A kuuluu
 if(h < r & kpB > r & kpA <= r){
  x <- sqrt(r^2 - h^2) - sqrt(abs(kpA^2 - h^2))
  ala <- KantaKorkeus(x,h) + 1/2 * r^2 * (acos(h/kpB) - acos(h/r))
 }

# A ja B kuuluvat ympyrään
 if(h < r & kpA <= r & kpB <= r){
  ala <- KantaKorkeus(AB, h)
 }
 
 tulos <- c(ala, (1/2*h*AB), (ala/(1/2*h*AB)),h)
}
}

tulos
}
#----------------------------------------------------------------------------



#leik <- NULL
# for(sade in 0:30){
# leik[sade+1] <- YmpyraNelikulmioAla(c(25,60),sade, c(24.9272,59.9944,24.9272+LOLisays,59.9944,
# 24.9272+LOLisays,59.9944 + LALisays,24.9272,59.9944+ LALisays))[3]
# }
#leik

#pisteLO <- kojetaulu[4,1]
#pisteLA <- kojetaulu[4,2]

#leik <- NULL
#for(sade in 0:30){
#leik[sade+1] <- YmpyraNelikulmioAla(c(25,60),sade, c(pisteLO,pisteLA,pisteLO+LOLisays,pisteLA,
#pisteLO+LOLisays,pisteLA + LALisays,pisteLO,pisteLA+ LALisays))[3]
#}
#leik

#kp <- c(25,60)
#kulmat <- c(pisteLO,pisteLA,pisteLO+LOLisays,pisteLA,
#pisteLO+LOLisays,pisteLA + LALisays,pisteLO,pisteLA+ LALisays)
#r <- 1
#KolmioYmpyranAla(kp[1], kp[2], r, kulmat[1],kulmat[2],kulmat[3],kulmat[4],kulmat[7],kulmat[8])
#KolmioYmpyranAla(kp[1], kp[2], r, kulmat[5],kulmat[6],kulmat[3],kulmat[4],kulmat[7],kulmat[8])
#kulmavek <- c(kulmat[1],kulmat[2],kulmat[3],kulmat[4],kulmat[7],kulmat[8])
#r <- 1
#kp1 <- kp[1]
#kp2 <- kp[2]
#kv1 <- kulmavek[1]
#kv2 <- kulmavek[2]
#kv3 <- kulmavek[3]
#kv4 <- kulmavek[4]
#kv5 <- kulmavek[5]
#kv6 <- kulmavek[6]


#----------------------------------------------------------------------------
#----------------------------------------------------------------------------
#kv on vektori, joka sisaltaa kolmion kulmat
KolmioYmpyranAla <- function(kp1, kp2, r, kv1,kv2,kv3,kv4,kv5,kv6)
{

 ala <- NULL
 tilavek <- KolmionSisalla(kp1,kp2,kv1,kv2,kv3,kv4,kv5,kv6)

 A <- tilavek[2:3]
 B <- tilavek[4:5]
 C <- tilavek[6:7]
 AB <- Etaisyys(6371,A[1],A[2],B[1],B[2],1)
 BC <- Etaisyys(6371,B[1],B[2],C[1],C[2],1)
 AC <- Etaisyys(6371,A[1],A[2],C[1],C[2],1)

#Lasketaan tarvittavien kolmioiden leikkaus
 alaKPAB <- AlaYmpyranKpKolmionKulma(kp1, kp2, r, A[1],A[2],B[1],B[2])[1]
 alaKPBC <- AlaYmpyranKpKolmionKulma(kp1, kp2, r, B[1],B[2],C[1],C[2])[1]
 alaKPAC <- AlaYmpyranKpKolmionKulma(kp1, kp2, r, A[1],A[2],C[1],C[2])[1]

#oikea kolmio ja piste kolmion sisalla
 if(tilavek[1]==0 & tilavek[8]==1)
 {
  ala[1] <- alaKPAB + alaKPAC + alaKPBC
  ala[2] <- HeronAla(AB,BC,AC)

  if(ala[1] < 0) ala[1] <- 0
  if(ala[1] > ala[2]) ala[1] <- ala[2]

  ala[3] <- ala[1]/ala[2]
 }

#Oikea kolmio ja piste kolmion ulkopuolella
 if(tilavek[1]==0 & tilavek[8]==0)
 {
  paikkavek <- PisteenAlue(kp1,kp2,tilavek[2:7])

  if(paikkavek[1] == 1)
  {
   ala[1] <- alaKPAC + alaKPBC - alaKPAB
  }

  if(paikkavek[1] == 2)
  {
   ala[1] <- alaKPAC - alaKPAB - alaKPBC
  }

  if(paikkavek[1] == 3)
  {
   ala[1] <- alaKPAB + alaKPAC - alaKPBC
  }

  if(paikkavek[1] == 4)
  {
   ala[1] <- alaKPAB - alaKPAC - alaKPBC
  }

  if(paikkavek[1] == 5)
  {
   ala[1] <- alaKPAB + alaKPBC - alaKPAC
  }

  if(paikkavek[1] == 6)
  {
   ala[1] <- alaKPBC - alaKPAC - alaKPAB
  }
  ala[2] <- HeronAla(AB,BC,AC)

  if(ala[1] < 0) ala[1] <- 0
  if(ala[1] > ala[2]) ala[1] <- ala[2]

  ala[3] <- ala[1]/ala[2]
 }

 ala
}


#----------------------------------------------------------------------------
#----------------------------------------------------------------------------
# Syötä kulmat kiertojärjestyksessä vastapäivään!!!!!!

YmpyraNelikulmioAla <-function(kp,r,kulmat)
{
tulos <- NULL
apu1 <- KolmioYmpyranAla(kp[1], kp[2], r, kulmat[1],kulmat[2],kulmat[3],kulmat[4],kulmat[7],kulmat[8])
apu2 <- KolmioYmpyranAla(kp[1], kp[2], r, kulmat[5],kulmat[6],kulmat[3],kulmat[4],kulmat[7],kulmat[8])
tulos[1:2] <- apu1[1:2] + apu2[1:2]
tulos[3] <- tulos[1]/tulos[2]
tulos
}
#-----------------------------------------------------------------------------













#leikkausData <- NULL
#osuus <- NULL
#for(i in 1:dim(kojetaulu3)[1]){
#for(i in 1:100){
#pisteLO <- NULL
#pisteLA <- NULL
#pisteLO <- kojetaulu[i,1]
#pisteLA <- kojetaulu[i,2]
#osuus[i] <- YmpyraNelikulmioAla(c(25,60),25, c(pisteLO,pisteLA,pisteLO+LOLisays,pisteLA,
#pisteLO+LOLisays,pisteLA + LALisays,pisteLO,pisteLA+ LALisays))[3]
#}
#leikkausData <- osuus*kojetaulu[,3]
#cbind(kojetaulu[,3],t(t(leikkausData)))
#summary(cbind(kojetaulu3[,3],t(t(leikkausData)))[,2])


#-----------------------------------------------------------------------------
#Laskeetaan pienempiä ympyröitä osuuksien avulla.
DataRuudut <- 10
# HUOM! tiedostoissa pääte .out on 50km ja .out.1 on 10km

if (DataRuudut == 10) paate <- ".out.1"
if (DataRuudut == 50) paate <- ".out"

SadeLisays <- 50
PieninSade <- 0
MaxSade <- 1700


LOLisays <- 0.1981
LOLisaysArvio <- 0.199 
LALisays <- 0.08995
LALisaysArvio <- 0.09

PisteKK <- c(25,60,22,60.5,22,63,26,62,30,63,24.5,65,29.5,65,24,67.5,28,69)

polku <- "N:\\Huippuyksikko\\Tutkimus\\R19_Kopra\\Data\\Model_Converted\\"
talletuspolku <- "N:\\Huippuyksikko\\Tutkimus\\R19_Kopra\\Mallit\\IntakeFraction\\Distance\\"


PopulaatioData <- NULL
PopulaatioData <- 
read.table(paste(polku,"Population_Europe\\",DataRuudut,"km\\ciesin_center_63lon_",DataRuudut,"km_2",sep=""))

PopulaatioData <- sPD



for( Piste in 4:4){
 LO <- PisteKK[(2*Piste-1)]
 LA <- PisteKK[(2*Piste)]

 SummaPixCi <- NULL

# Seuraavassa loopissa kk on kuukausi
 for(kk in 1:12){
  tiedosto1<-paste("joinfilelist_",kk,"_PM_0_1.joined", paate,sep="")
  tiedosto2<-paste("joinfilelist_",kk,"_PM_0_1_1.joined", paate,sep="")
  tiedosto3<-paste("joinfilelist_",kk,"_PM_1_2_5.joined", paate,sep="")

  PitoisuusData1 <- NULL
  PitoisuusData1 <-read.table(paste(polku,"Concentration_Point\\",tiedosto1,sep=""))

  PitoisuusData2 <- NULL
  PitoisuusData2 <-read.table(paste(polku,"Concentration_Point\\",tiedosto2,sep=""))

  PitoisuusData3 <- NULL
  PitoisuusData3 <-read.table(paste(polku,"Concentration_Point\\",tiedosto3,sep=""))

  PM2.5PitoisuusData<- NULL
  PM2.5PitoisuusData <- PitoisuusData1 + PitoisuusData2 + PitoisuusData3

  apuPM <- NULL
  apuPM <- PM2.5PitoisuusData[(PM2.5PitoisuusData[,(1 + 2*Piste)]!=0 | PM2.5PitoisuusData[,(2 + 2*Piste)]!=0)
  & PopulaatioData[,3]!=0,]

  apuPD <- NULL
  apuPD <- PopulaatioData[PopulaatioData[,3]!=0 &
  (PM2.5PitoisuusData[,(1 + 2*Piste)]!=0 | PM2.5PitoisuusData[,(2 + 2*Piste)]!=0),]

  apuHI <- NULL
  apuLO <- NULL
  apuSumma <- NULL


  sadeindeksi <- 0

  for(Sade in seq(MaxSade, PieninSade, by= -SadeLisays )){
   apu1 <- NULL


   sadeindeksi <- sadeindeksi + 1
#-------------------------------------------------------------
#Etäisyyden laskentatapa

   for(k in 1:dim(apuPD)[1]){
    if( (apuPD[k,1] >= LO   &  apuPD[k,2] >= LA) | (apuPD[k,1] < LO   &  apuPD[k,2] < LA ) ){
     apu1[k] <- 
     (Etaisyys(6371,LO, LA, (apuPD[k,1] + LOLisays), (apuPD[k,2]),1) < Sade &
     Etaisyys(6371,LO, LA, (apuPD[k,1]), (apuPD[k,2] + LALisays),1) < Sade)
    }

    if((apuPD[k,1] < LO   &  apuPD[k,2] >= LA) | (apuPD[k,1] >= LO   &  apuPD[k,2] < LA ) ){
     apu1[k] <- 
     (Etaisyys(6371,LO, LA, (apuPD[k,1]), (apuPD[k,2]),1) < Sade &
     Etaisyys(6371,LO, LA, (apuPD[k,1] + LOLisays), (apuPD[k,2] + LALisays),1) < Sade )
    }
   }
#  for(k in 1:dim(apuPD)[1]){
#    if( (apuPD[k,1] >= LO   &  apuPD[k,2] >= LA) | (apuPD[k,1] < LO   &  apuPD[k,2] < LA ) ){
#     apu1[k] <- 
#     (Etaisyys(6371,LO, LA, (apuPD[k,1] + LOLisays), (apuPD[k,2]),1) < Sade+25 &
#     Etaisyys(6371,LO, LA, (apuPD[k,1]), (apuPD[k,2] + LALisays),1) < Sade+25)
#    }

#    if((apuPD[k,1] < LO   &  apuPD[k,2] >= LA) | (apuPD[k,1] >= LO   &  apuPD[k,2] < LA ) ){
#     apu1[k] <- 
#     (Etaisyys(6371,LO, LA, (apuPD[k,1]), (apuPD[k,2]),1) < Sade+25 &
#     Etaisyys(6371,LO, LA, (apuPD[k,1] + LOLisays), (apuPD[k,2] + LALisays),1) < Sade+25)
#    }
#   }
   apuPM <- apuPM[apu1,]
   apuPD <- apuPD[apu1,]

#   osuus <- NULL
#   apuPDOsuus <- NULL
#   for(i in 1:dim(apuPD)[1]){
#    pisteLO <- NULL
#    pisteLA <- NULL
#    pisteLO <- apuPD[i,1]
#    pisteLA <- apuPD[i,2]
#    osuus[i] <- YmpyraNelikulmioAla(c(LO,LA),Sade, c(pisteLO,pisteLA,pisteLO+LOLisays,pisteLA,
#    pisteLO+LOLisays,pisteLA + LALisays,pisteLO,pisteLA+ LALisays))[3]
#   }
#   apuPDOsuus <- apuPD[,3]*t(t(osuus))
#
#   apuHI[sadeindeksi] <- sum( apuPM[,(1 + 2*Piste)] * apuPDOsuus )
#   apuLO[sadeindeksi] <- sum( apuPM[,(2 + 2*Piste)] * apuPDOsuus )

   apuHI[sadeindeksi] <- sum( apuPM[,(1 + 2*Piste)] * apuPD[,3] )
   apuLO[sadeindeksi] <- sum( apuPM[,(2 + 2*Piste)] * apuPD[,3] )

#Etäisyyden laskentatapa
#--------------------------------------------------------------------------
   aputaulu <- c(Sade, dim(apuPD)[1])

#   TalletusKohde <- paste(talletuspolku,"Piste_",Piste,"\\IF_Piste_",Piste,"_kk_",kk,"_Sade_",Sade,"_Distance_Osuus_kokeilu.txt",sep="")
   TalletusKohde <- paste(talletuspolku,"Piste_",Piste,"\\IF_Piste_",Piste,"_kk_",kk,"_Sade_",Sade,"_Distance.txt",sep="")
   write.table(t(aputaulu), file = TalletusKohde, row.names = FALSE, col.names = FALSE)
  }
 
  ylaraja <- length(apuHI)
  for(i in 1: ylaraja){
   apuSumma[(2*i - 1)] <- apuHI[(ylaraja - i + 1)]
   apuSumma[2*i] <- apuLO[(ylaraja - i + 1)]
  }

  SummaPixCi <- rbind(SummaPixCi, t(apuSumma))
 }

# BR <- 20 m^3/päivä
# BR <- 20/(24 * 60 * 60) m^3/s
 BR <- 20/(24 * 60 * 60)

# Q <- 3 * 1/(365 * 24 * 60 * 60)kg/s
 Q <- 3/(365 * 24 * 60 * 60)

 IF <- SummaPixCi * BR/Q

 TalletusKohde <- paste(talletuspolku,"IF_Piste_",Piste,"_Distance.txt",sep="")
# TalletusKohde <- paste(talletuspolku,"IF_Piste_",Piste,"_Distance_Osuus_Kokeilu_SadeVali_1_Suurin_100.txt",sep="")
 write.table(IF, file = TalletusKohde, row.names = FALSE, col.names = FALSE)
}



#-----------------------------------------------------------------------------
# valitaan ruudut annetulla säteellä. 
DataRuudut <- 10
MaxSade
# HUOM! tiedostoissa pääte .out on 50km ja .out.1 on 10km
if (DataRuudut == 10) paate <- ".out.1"
if (DataRuudut == 50) paate <- ".out"
SadeLisays <- 5
PieninSade <- 0
MaxSade <- 100
LOLisays <- 0.1981
LOLisaysArvio <- 0.199 
LALisays <- 0.08995
LALisaysArvio <- 0.09
PisteKK <- c(25,60,22,60.5,22,63,26,62,30,63,24.5,65,29.5,65,24,67.5,28,69)
polku <- "N:\\Huippuyksikko\\Tutkimus\\R19_Kopra\\Data\\Model_Converted\\"
talletuspolku <- "N:\\Huippuyksikko\\Tutkimus\\R19_Kopra\\Mallit\\IntakeFraction\\Distance\\"
PopulaatioData <- NULL
PopulaatioData <- 
read.table(paste(polku,"Population_Europe\\",DataRuudut,"km\\ciesin_center_63lon_",DataRuudut,"km_2",sep=""))
  apuPD <- NULL
  apuPD <- PopulaatioData[PopulaatioData[,3]!=0 ,]
  aputaulu <- NULL

  for(Sade in seq(MaxSade, PieninSade, by= -SadeLisays )){
   apu1 <- NULL
   apu1[1] <- Sade
   for(k in 1:dim(apuPD)[1]){
    if( (apuPD[k,1] >= LO   &  apuPD[k,2] >= LA) | (apuPD[k,1] < LO   &  apuPD[k,2] < LA ) ){
     apu1[k] <- 
     (Etaisyys(6371,LO, LA, (apuPD[k,1] + LOLisays), (apuPD[k,2]),1) < Sade+25 &
     Etaisyys(6371,LO, LA, (apuPD[k,1]), (apuPD[k,2] + LALisays),1) < Sade+25)
    }
    if((apuPD[k,1] < LO   &  apuPD[k,2] >= LA) | (apuPD[k,1] >= LO   &  apuPD[k,2] < LA ) ){
     apu1[k] <- 
     (Etaisyys(6371,LO, LA, (apuPD[k,1]), (apuPD[k,2]),1) < Sade+25 &
     Etaisyys(6371,LO, LA, (apuPD[k,1] + LOLisays), (apuPD[k,2] + LALisays),1) < Sade+25)
    }
   }
   aputaulu <- rbind(aputaulu, apu1)
   TalletusKohde <- paste(talletuspolku,"PopulData_Sateittain.txt",sep="")
   write.table(t(aputaulu), file = TalletusKohde, row.names = FALSE, col.names = FALSE)
  }
#---------------------------------------------------------------------------------------------------------




kulmat <-  c(24.9272,59.9944,24.9272+LOLisays,59.9944,
 24.9272+LOLisays,59.9944 + LALisays,24.9272,59.9944+ LALisays)
kp <- c(25,60)
r <- 1

leik <- NULL
 for(sade in 0:50){
 leik[sade] <- YmpyraNelikulmioAla(c(25,60),sade, c(24.9272,59.9944,24.9272+LOLisays,59.9944,
 24.9272+LOLisays,59.9944 + LALisays,24.9272,59.9944+ LALisays))[3]
 }
leik

kulmat <- c(24.9272,59.9944,24.9272+LOLisays,59.9944,
 24.9272+LOLisays,59.9944 + LALisays,24.9272,59.9944+ LALisays)
kp <- c(25,60)

r<- 7



  leikkausData <- NULL
  leikkausData <- t(t(kojetaulu2[,3]))
  for(Sade in seq(1, 40, by= 1 )){
   osuus <- NULL
   for(i in 1:dim(leikkausData)[1]){
    pisteLO <- NULL
    pisteLA <- NULL
    pisteLO <- kojetaulu2[i,1]
    pisteLA <- kojetaulu2[i,2]
    osuus[i] <- YmpyraNelikulmioAla(c(LO,LA),Sade, c(pisteLO,pisteLA,pisteLO+LOLisays,pisteLA,
    pisteLO+LOLisays,pisteLA + LALisays,pisteLO,pisteLA+ LALisays))[3]
   }
   leikkausData <- cbind(leikkausData,kojetaulu2[,3]*t(t(osuus)))
  }

alue <- NULL
for(i in 1:41){
alue[i] <- sum(leikkausData[,i]) 
}



#---------------------------------------------------------
#---------------------------------------------------------
#
# Vanhaa tavaraa!!!!


#-----------------------------------------------------------------------------
#Laskee ympyrän ja kolmion leikkauksen alan, kun ympyrän keskipiste on kolmion sisällä
AlaYmpyranKpKolmionSisalla <- function(kp1, kp2, r, kv1,kv2,kv3,kv4,kv5,kv6)
{
tilavek <- KolmionSisalla(kp1,kp2,kv1,kv2,kv3,kv4,kv5,kv6)
hNollaraja <- 10^(-6)
kp <- c(kp1,kp2)
A <- tilavek[2:3]
B <- tilavek[4:5]
C <- tilavek[6:7]
KolmionKulmat <- rbind(A,B,C,A)
#oikea kolmioa ja piste kolmion sisalla
if(tilavek[1]==0 & tilavek[8]==1)
{
AB <- Etaisyys(6371,tilavek[2],tilavek[3],tilavek[4],tilavek[5],1)
BC <- Etaisyys(6371,tilavek[4],tilavek[5],tilavek[6],tilavek[7],1)
AC <- Etaisyys(6371,tilavek[2],tilavek[3],tilavek[6],tilavek[7],1)
etA <- Etaisyys(6371, kp1, kp2, tilavek[2], tilavek[3],1)
etB <- Etaisyys(6371, kp1, kp2, tilavek[4], tilavek[5],1)
etC <- Etaisyys(6371, kp1, kp2, tilavek[6], tilavek[7],1)

etv <- c(etA, etB, etC, etA)
valit <- c(AB, BC, AC)

#h <- c(-100,-100,-100)
h <- c(1,1,1)
# Tapaukset, joissa ymp. kp on kolmion kulmapiste
if( (kp == A)[1] & (kp == A)[2]){ h[1] <- 0 
  h[3] <- 0
}
if( (kp == B)[1] & (kp == B)[2]){ h[1] <- 0 
  h[2] <- 0
}
if( (kp == C)[1] & (kp == C)[2]){ h[2] <- 0 
  h[3] <- 0
}

# Tapaukset, joissa ymp. kp on kolmion kulmapisteiden välisellä suoralla
if(A[1]!=B[1])
{
 k <- (B[2] - A[2])/(B[1] - A[1])
 y <- k* (kp[1] - A[1]) + A[2]
 if( abs(kp[2] - y) < hNollaraja) h[1] <- 0
}
if( (A[1] == B[1] & kp[1]==A[1]))
{
h[1] <- 0
}

if(B[1]!=C[1])
{
 k <- (C[2] - B[2])/(C[1] - B[1])
 y <- k* (kp[1] - B[1]) + B[2]
 if( abs(kp[2] - y) < hNollaraja) h[2] <- 0
}
if( (B[1] == C[1] & kp[1]==B[1]))
{
h[2] <- 0
}

if(C[1]!=A[1])
{
 k <- (A[2] - C[2])/(A[1] - C[1])
 y <- k* (kp[1] - C[1]) + C[2]
 if( abs(kp[2] - y) < hNollaraja) h[3] <- 0
}
if( (C[1] == A[1] & kp[1]==C[1]))
{
h[3] <- 0
}

ala <- NULL
for(i in 1:3){
if(h[i] >= hNollaraja){
ala[i] <- AlaYmpyranKpKolmionKulma(kp[1], kp[2], r, KolmionKulmat[i,1],
KolmionKulmat[i,2],KolmionKulmat[(i+1),1],KolmionKulmat[(i+1),2])[1]
}
if(h[i] < hNollaraja){
ala[i] <- 0
}
}

tulos <- NULL
tulos[1] <- sum(ala)
tulos[2] <- KantaKorkeus(valit[1],
     sqrt(valit[3]^2 - ((valit[1]^2 + valit[3]^2 - valit[2]^2)/(2*valit[1]))^2))
tulos[3] <- tulos[1]/tulos[2]
}

if(tilavek[1]==1 & length(tilavek)==1)
{
tulos <- c(0,0,0)
}
if(length(tilavek) > 1 & tilavek[8]== 0)
{
tulos <- c(-1,-1,-1)
}

tulos
}


#-------------------------------------------------------------------------------------------
#Pisteiden koordinaatit syötettävä kiertosuunnassa vastapäivään!!!!
YmpyraNelikulmioKpSisalla <- function(kpLO, kpLA, r, ALO, ALA, BLO, BLA, CLO, CLA, DLO, DLA)
{
tyovek <- c(kpLO, kpLA, r, ALO, ALA, BLO, BLA, CLO, CLA, DLO, DLA)
kp <- tyovek[1:2]
KolmionKulmat <- rbind(c(tyovek[4:5],tyovek[6:7]),
c(tyovek[6:7],tyovek[8:9]),
c(tyovek[8:9],tyovek[10:11]),
c(tyovek[10:11],tyovek[4:5]))

ala <- NULL
for(i in 1:4){
ala[i] <-  AlaYmpyranKpKolmionKulma(kp[1], kp[2], r,
KolmionKulmat[i,1], KolmionKulmat[i,2],
KolmionKulmat[i,3],KolmionKulmat[i,4])[1]
}

AB <- Etaisyys(6371, tyovek[4],tyovek[5], tyovek[6], tyovek[7],1)
BC <- Etaisyys(6371, tyovek[6], tyovek[7], tyovek[8], tyovek[9],1)
CD <- Etaisyys(6371, tyovek[8], tyovek[9], tyovek[10], tyovek[11],1)
DA <- Etaisyys(6371, tyovek[10], tyovek[11], tyovek[4], tyovek[5],1)
AC <- Etaisyys(6371, tyovek[4], tyovek[5], tyovek[8], tyovek[9],1)
tulos <- NULL
tulos[1] <- sum(ala)
tulos[2] <- HeronAla(AB, BC, AC) + HeronAla(AC, CD, DA)
tulos[3] <- tulos[1]/tulos[2]
tulos
}


#----------------------------------------------------------------------------
#kv on vektori, joka sisaltaa kolmion kulmat
KolmioYmpyranAla <- function(kp1, kp2, r, kv1,kv2,kv3,kv4,kv5,kv6)
{


ala <- NULL
tilavek <- KolmionSisalla(kp1,kp2,kv1,kv2,kv3,kv4,kv5,kv6)

#oikea kolmioa ja piste kolmion sisalla
if(tilavek[1]==0 & tilavek[8]==1)
{
ala <- AlaYmpyranKpKolmionSisalla(kp1, kp2, r, kv1,kv2,kv3,kv4,kv5,kv6)
}

#Oikea kolmio ja piste kolmion ulkopuolella
if(tilavek[1]==0 & tilavek[8]==0)
{
paikkavek <- PisteenAlue(kp1,kp2,tilavek[2:7])
A <- tilavek[2:3]
B <- tilavek[4:5]
C <- tilavek[6:7]
AB <- Etaisyys(6371,A[1],A[2],B[1],B[2],1)
BC <- Etaisyys(6371,B[1],B[2],C[1],C[2],1)
AC <- Etaisyys(6371,A[1],A[2],C[1],C[2],1)

if(paikkavek[1] == 1)
{

# kp = (1 - paikkavek[3])*C + paikkavek[3]*(paikkavek[2]*B + (1-paikkavek[2])*A)
uusiA <- (1 - (paikkavek[3] + 1))*C + (paikkavek[3] + 1)*A
uusiB <- (1 - (paikkavek[3] + 1))*C + (paikkavek[3] + 1)*B
uusiC <- C
#Pisteiden koordinaatit syötettävä kiertosuunnassa vastapäivään!!!!!!!!!!

ala[1] <- AlaYmpyranKpKolmionSisalla(kp1, kp2, r,
 uusiB[1], uusiB[2], uusiA[1], uusiA[2], uusiC[1], uusiC[2])[1]- 
YmpyraNelikulmioKpSisalla(kp1, kp2, r,
   uusiB[1], uusiB[2], uusiA[1], uusiA[2], A[1], A[2], B[1], B[2])[1]
}

if(paikkavek[1] == 2)
{
#iso kolmio - pikkukolmiot
ala[1] <- AlaYmpyranKpKolmionSisalla(kp1, kp2, r, A[1], A[2], C[1],C[2], kp1, kp2)[1]-
AlaYmpyranKpKolmionSisalla(kp1, kp2, r, B[1], B[2], C[1],C[2], kp1, kp2)[1]-
AlaYmpyranKpKolmionSisalla(kp1, kp2, r, A[1], A[2], B[1],B[2], kp1, kp2)[1]
}

if(paikkavek[1] == 3)
{
uusiA <- A
uusiB <- (1 - (paikkavek[3] + 1))*A + (paikkavek[3] + 1)*B
uusiC <- (1 - (paikkavek[3] + 1))*A + (paikkavek[3] + 1)*C
#Pisteiden koordinaatit syötettävä kiertosuunnassa vastapäivään!!!!!!!!!!
ala[1] <- AlaYmpyranKpKolmionSisalla(kp1, kp2, r,
 uusiB[1], uusiB[2], uusiA[1], uusiA[2], uusiC[1], uusiC[2])[1]- 
YmpyraNelikulmioKpSisalla(kp1, kp2, r,
   uusiB[1], uusiB[2], B[1], B[2], C[1], C[2], uusiC[1], uusiC[2])[1]
}

if(paikkavek[1] == 4)
{
#iso kolmio - pikkukolmiot
ala[1] <- AlaYmpyranKpKolmionSisalla(kp1, kp2, r, A[1], A[2], B[1],B[2], kp1, kp2)[1]-
AlaYmpyranKpKolmionSisalla(kp1, kp2, r, B[1], B[2], C[1],C[2], kp1, kp2)[1]-
AlaYmpyranKpKolmionSisalla(kp1, kp2, r, A[1], A[2], C[1],C[2], kp1, kp2)[1]
}

if(paikkavek[1] == 5)
{
uusiA <- (1 - (paikkavek[3] + 1))*B + (paikkavek[3] + 1)*A
uusiB <- B
uusiC <- (1 - (paikkavek[3] + 1))*B + (paikkavek[3] + 1)*C
#Pisteiden koordinaatit syötettävä kiertosuunnassa vastapäivään!!!!!!!!!!

ala[1] <- AlaYmpyranKpKolmionSisalla(kp1, kp2, r,
 uusiB[1], uusiB[2], uusiA[1], uusiA[2], uusiC[1], uusiC[2])[1]- 
YmpyraNelikulmioKpSisalla(kp1, kp2, r,
      uusiA[1], uusiA[2], uusiC[1], uusiC[2], C[1], C[2], A[1], A[2])[1]

}

if(paikkavek[1] == 6)
{
#iso kolmio - pikkukolmiot
ala[1] <- AlaYmpyranKpKolmionSisalla(kp1, kp2, r, C[1], C[2], B[1],B[2], kp1, kp2)[1]-
AlaYmpyranKpKolmionSisalla(kp1, kp2, r, B[1], B[2], A[1],A[2], kp1, kp2)[1]-
AlaYmpyranKpKolmionSisalla(kp1, kp2, r, A[1], A[2], C[1],C[2], kp1, kp2)[1]
}

ala[2] <- HeronAla(AB,BC,AC)
ala[3] <- ala[1]/ala[2]
}

ala
}


#----------------------------------------------------------------------------
# Syötä kulmat kiertojärjestyksessä vastapäivään!!!!!!

YmpyraNelikulmioAla <-function(kp,r,kulmat)
{
tulos <- NULL
apu1 <- KolmioYmpyranAla(kp[1], kp[2], r, kulmat[1],kulmat[2],kulmat[3],kulmat[4],kulmat[7],kulmat[8])
apu2 <- KolmioYmpyranAla(kp[1], kp[2], r, kulmat[5],kulmat[6],kulmat[3],kulmat[4],kulmat[7],kulmat[8])
tulos[1:2] <- apu1[1:2] + apu2[1:2]
tulos[3] <- tulos[1]/tulos[2]
tulos
}

Polku <- "C:\\omat\\dataa\\Concentration_SOX\\"
PaastoLahde <- "SO2"
 

postscript(file=paste(Polku,"IFConcentration",PaastoLahde,".ps",sep=""), width=10,height=10,horizontal=F)


IFmat10<-read.table(paste(Polku,"IFConcentration",PaastoLahde,"10km.txt",sep=""))
IFmat50<-read.table(paste(Polku,"IFConcentration",PaastoLahde,"50km.txt",sep=""))

#Ensin piirretään pisteen High sitten Low (=rgb(0,0,0) = musta)

for(i in 1:9){
piste <- paste("Piste ", i, sep="")
plot(c(1:12),IFmat10[,(2*i - 1)],type="b", xlim=c(1,12), ylim=c(0,1*10^(-5)),
xlab="kuukausi", ylab="IF 10km", col=rgb(1,0,0), main = paste(piste,"High = punainen, Low = musta"))
par(new=T)
plot(c(1:12),IFmat10[,(2*i)],type="b", xlim=c(1,12), ylim=c(0,1*10^(-5)),
xlab="kuukausi", ylab="IF 10km", col=rgb(0,0,0), main = paste(piste,"High = punainen, Low = musta"))
}


for(i in 1:9){
piste <- paste("Piste ", i, sep="")
plot(c(1:12),IFmat50[,(2*i - 1)],type="b", xlim=c(1,12), ylim=c(0,1*10^(-5)),
xlab="kuukausi", ylab="IF 50km", col=rgb(1,0,0), main = paste(piste,"High = punainen, Low = musta"))
par(new=T)
plot(c(1:12),IFmat50[,(2*i)],type="b", xlim=c(1,12), ylim=c(0,1*10^(-5)),
xlab="kuukausi", ylab="IF 50km", col=rgb(0,0,0), main = paste(piste,"High = punainen, Low = musta"))
}

# musta on 0,0,0

dev.off()



#_-----------------------------------------------------------

Polku <- "C:\\omat\\dataa\\Concentration_SOX\\"
PaastoLahde <- "SO4"
 

postscript(file=paste(Polku,"IFConcentration",PaastoLahde,".ps",sep=""), width=10,height=10,horizontal=F)


IFmat10<-read.table(paste(Polku,"IFConcentration",PaastoLahde,"10km.txt",sep=""))
IFmat50<-read.table(paste(Polku,"IFConcentration",PaastoLahde,"50km.txt",sep=""))

#Ensin piirretään pisteen High sitten Low (=rgb(0,0,0) = musta)

for(i in 1:9){
piste <- paste("Piste ", i, sep="")
plot(c(1:12),IFmat10[,(2*i - 1)],type="b", xlim=c(1,12), ylim=c(0,1.5*10^(-4)),
xlab="kuukausi", ylab="IF 10km", col=rgb(1,0,0), main = paste(piste,"High = punainen, Low = musta"))
par(new=T)
plot(c(1:12),IFmat10[,(2*i)],type="b", xlim=c(1,12), ylim=c(0,1.5*10^(-4)),
xlab="kuukausi", ylab="IF 10km", col=rgb(0,0,0), main = paste(piste,"High = punainen, Low = musta"))
}


for(i in 1:9){
piste <- paste("Piste ", i, sep="")
plot(c(1:12),IFmat50[,(2*i - 1)],type="b", xlim=c(1,12), ylim=c(0,1.5*10^(-4)),
xlab="kuukausi", ylab="IF 50km", col=rgb(1,0,0), main = paste(piste,"High = punainen, Low = musta"))
par(new=T)
plot(c(1:12),IFmat50[,(2*i)],type="b", xlim=c(1,12), ylim=c(0,1.5*10^(-4)),
xlab="kuukausi", ylab="IF 50km", col=rgb(0,0,0), main = paste(piste,"High = punainen, Low = musta"))
}

# musta on 0,0,0

dev.off()

# Datat -> N:\Huippuyksikko\Tutkimus\R19_Kopra\Data\Model_Converted


#--------------


# Concentration_SOX

# 50x50 ruudut

AlkuAika <- Sys.time()

SummaPixCi50kmSO2 <- matrix(nrow=12,ncol=18)
SummaPixCi50kmSO4 <- matrix(nrow=12,ncol=18)

polku <- "N:\\Huippuyksikko\\Tutkimus\\R19_Kopra\\Data\\Model_Converted\\"
talletuspolku <- "c:\\omat\\dataa\\concentration_SOX\\"

PopulaatioData <- NULL
PopulaatioData<-read.table(paste(polku,"Population_Europe\\50km\\ciesin_center_63lon_50km_2",sep=""))


for(i in 1:12){


tiedosto1<-paste("joinfilelist_",i,"_SO2.joined.out",sep="")
tiedosto2<-paste("joinfilelist_",i,"_SO4.joined.out",sep="")



PitoisuusDataSO2 <- NULL
PitoisuusDataSO2 <-read.table(paste(polku,"Concentration_SOX\\",tiedosto1,sep=""))

PitoisuusDataSO4 <- NULL
PitoisuusDataSO4 <-read.table(paste(polku,"Concentration_SOX\\",tiedosto2,sep=""))


for(j in 1:18){
 SummaPixCi50kmSO2[i,j] <- sum( PitoisuusDataSO2[,(j+2)] * PopulaatioData[,3] )
 SummaPixCi50kmSO4[i,j] <- sum( PitoisuusDataSO4[,(j+2)] * PopulaatioData[,3] )
}
}



# BR <- 20 m^3/päivä
# BR <- 20/(24 * 60 * 60) m^3/s

BR <- 20/(24 * 60 * 60)


# Q <- 1 * 1/(365 * 24 * 60 * 60) mol/s

Q <- 1/(365 * 24 * 60 * 60)

IFConcentrationPoint50kmSO2 <- SummaPixCi50kmSO2 * BR/Q
IFConcentrationPoint50kmSO4 <- SummaPixCi50kmSO4 * BR/Q


TalletusKohde <- paste(talletuspolku,"IFConcentrationSO250km.txt",sep="")
write.table(IFConcentrationPoint50kmSO2, file = TalletusKohde, row.names = FALSE, col.names = FALSE)

TalletusKohde <- paste(talletuspolku,"IFConcentrationSO450km.txt",sep="")
write.table(IFConcentrationPoint50kmSO4, file = TalletusKohde, row.names = FALSE, col.names = FALSE)


LoppuAika <- Sys.time()

LoppuAika - AlkuAika


rm(SummaPixCi50kmSO2, SummaPixCi50kmSO4,polku, talletuspolku, PopulaatioData, AlkuAika, LoppuAika, tiedosto1, tiedosto2, PitoisuusDataSO2, PitoisuusDataSO4, Q,
   BR, IFConcentrationPoint50kmSO2, IFConcentrationPoint50kmSO4, TalletusKohde)



#-----------------------------
# 10x10 ruudut

AlkuAika <- Sys.time()
SummaPixCi10kmSO2 <- matrix(nrow=12,ncol=18)
SummaPixCi10kmSO4 <- matrix(nrow=12,ncol=18)
polku <- "N:\\Huippuyksikko\\Tutkimus\\R19_Kopra\\Data\\Model_Converted\\"
talletuspolku <- "c:\\omat\\dataa\\concentration_SOX\\"
PopulaatioData <- NULL
PopulaatioData<-read.table(paste(polku,"Population_Europe\\10km\\ciesin_center_63lon_10km_2",sep=""))

for(i in 1:12){

tiedosto1<-paste("joinfilelist_",i,"_SO2.joined.out.1",sep="")
tiedosto2<-paste("joinfilelist_",i,"_SO4.joined.out.1",sep="")

PitoisuusDataSO2 <- NULL
PitoisuusDataSO2 <-read.table(paste(polku,"Concentration_SOX\\",tiedosto1,sep=""))

PitoisuusDataSO4 <- NULL
PitoisuusDataSO4 <-read.table(paste(polku,"Concentration_SOX\\",tiedosto2,sep=""))

for(j in 1:18){
 SummaPixCi10kmSO2[i,j] <- sum( PitoisuusDataSO2[,(j+2)] * PopulaatioData[,3] )
 SummaPixCi10kmSO4[i,j] <- sum( PitoisuusDataSO4[,(j+2)] * PopulaatioData[,3] )
}
}


# BR <- 20 m^3/päivä
# BR <- 20/(24 * 60 * 60) m^3/s

BR <- 20/(24 * 60 * 60)


# Q <- 1/(365 * 24 * 60 * 60) mol/s

Q <- 1/(365 * 24 * 60 * 60)

IFConcentrationPoint10kmSO2 <- SummaPixCi10kmSO2 * BR/Q
IFConcentrationPoint10kmSO4 <- SummaPixCi10kmSO4 * BR/Q


TalletusKohde <- paste(talletuspolku,"IFConcentrationSO210km.txt",sep="")
write.table(IFConcentrationPoint10kmSO2, file = TalletusKohde, row.names = FALSE, col.names = FALSE)

TalletusKohde <- paste(talletuspolku,"IFConcentrationSO410km.txt",sep="")
write.table(IFConcentrationPoint10kmSO4, file = TalletusKohde, row.names = FALSE, col.names = FALSE)


LoppuAika <- Sys.time()

LoppuAika - AlkuAika


rm(SummaPixCi10kmSO2, SummaPixCi10kmSO4,polku, talletuspolku, PopulaatioData, AlkuAika, LoppuAika, tiedosto1, tiedosto2, PitoisuusDataSO2, PitoisuusDataSO4, Q,
   BR, IFConcentrationPoint10kmSO2, IFConcentrationPoint10kmSO4, TalletusKohde)

#7.3.2008, Marko Tainio

# This model calculates intake fractions (iF) for European
# primary PM2.5 emissions. Description of model:
# - Emission source: European primary PM2.5 emissions from different countries
# - Dispersion: Dispersion over the Europe
# - Population: European population divided to different countries (38)
# (Note: the population data does not include Finnish population!!)

#Work folders - the location of data and where the output is recorded
wrk_input = "N:/Huippuyksikko/Tutkimus/R19_Kopra/Data 2/"
wrk_output = "N:/Huippuyksikko/Tutkimus/R19_Kopra/Mallit/IntakeFraction2/Result/"

#Breathing rate, 20 m^3/day -> unit g/s
BR = 20/(24 * 60 * 60)

# Calculation for different countries. The iF's are calculated separately 
# for all emission sources (countires). There are 48 sources in data.
for (k in 1:48) {

# This reads from file the titles of countries to be used in calculations
maa = read.table("N:/Huippuyksikko/Tutkimus/R19_Kopra/Mallit/IntakeFraction2/EUR_country.txt", sep="")
maa = maa[,k]

#Emission volume - this part estimates emission volume in unit g/s
Q = read.csv("N:/Huippuyksikko/Tutkimus/R19_Kopra/Mallit/IntakeFraction2/Eur_emission_05.csv", header=T)

#Concentration data
#The order in list debends on the prepared emission data
rm(conc)
conc = "N:/Huippuyksikko/Tutkimus/R19_Kopra/Data 2/Concentration_primary_Europe_2000/EMEP_monthly_PM2_5_2000/"
conc = paste(conc,"EMEP_2000_PM_2_5_hirlam_v3_8_1_",maa,"_200012.txt", sep="")
conc = read.table(conc, header=T, sep="")

#Population data:
rm(pop)
pop = paste(wrk_input, "Population_Combined/Pop_eur00_all3.txt", sep="")
pop = read.table(pop, header=T, sep="")

#Co-ordination check-in - this check that co-ordinate system is same
#in population and concentration data.
check = (conc["lon"] - pop["Lo"]) + (conc["lat"] - pop["La"])
sum(check)

#This part slice the concentration data from the data file
#(that is, LO and LA data is sliced oof)
month = c("jan", "feb", "mar", "apr", "may", "jun", "jul", "aug", "sep", "oct", "nov", "dec", "annual")
conc = data.frame(conc[month])

#Unit change in concentration data: kg/m^3 -> g/m^3
conc = conc * 1000

#This part slice the population data from the data file
#(that is, LO and LA data is sliced off)
rm(country_eur)
country_eur = c("Sum_Austri", "Sum_Belgiu", "Sum_Cyprus", "Sum_Czech_", "Sum_Denmar", "Sum_Estoni", "Sum_France", "Sum_German", "Sum_Greece", "Sum_Hungar", "Sum_Irish_", "Sum_Italy", "Sum_Latvia", "Sum_Lithua", "Sum_Luxemb", "Sum_Malta", "Sum_Nether", "Sum_Poland", "Sum_Portug", "Sum_Romani", "Sum_Slovak", "Sum_Sloven", "Sum_Spain", "Sum_UK", "Sum_Sweden", "Sum_Bulgar", "Sum_Norway", "Sum_Belaru", "Sum_Russia", "Sum_Switze", "Sum_Ukrain", "Sum_Moldov", "Sum_Croati", "Sum_Serbia", "Sum_Albani", "Sum_Macedo", "Sum_Turkey", "Sum_bosni", "Sum_Finl", "Sum_all")
pop = data.frame(pop[country_eur])

#X must be defined before for -loop
x = NULL

#Month/concentration - month 13 in annual mean
for (i in 1:13) {

#Removal of huge values
#There are some weird very high values in some concentration datasets.
#This code replace those values with value 0.
y = conc[,i]
y = ifelse(y > 100, 0, y)

#Country/population - there are 40 different population data files
for (j in 1:40) {

x[j] = sum(y * pop[,j])

}

#Calculation of iF
x = (x * BR) / Q[k,i]

#Names for the data
names(x) = c("Sum_Austri", "Sum_Belgiu", "Sum_Cyprus", "Sum_Czech_", "Sum_Denmar", "Sum_Estoni", "Sum_France", "Sum_German", "Sum_Greece", "Sum_Hungar", "Sum_Irish_", "Sum_Italy", "Sum_Latvia", "Sum_Lithua", "Sum_Luxemb", "Sum_Malta", "Sum_Nether", "Sum_Poland", "Sum_Portug", "Sum_Romani", "Sum_Slovak", "Sum_Sloven", "Sum_Spain", "Sum_UK", "Sum_Sweden", "Sum_Bulgar", "Sum_Norway", "Sum_Belaru", "Sum_Russia", "Sum_Switze", "Sum_Ukrain", "Sum_Moldov", "Sum_Croati", "Sum_Serbia", "Sum_Albani", "Sum_Macedo", "Sum_Turkey", "Sum_bosni", "Sum_Finl", "Sum_all")

#Recording of results
write.table(x, paste(wrk_output, file = "European_emission/2000/if_emis_eur_year_2000_country_",maa,"_month_",i,"_pop_all.txt", sep=""), row.names = TRUE, col.names = FALSE)

}

}

#Remove of all data from memory
rm(wrk_input, wrk_output, BR, k, maa, Q, pop, check, conc, month, country_eur, x, i, y, j)

# 18.4.2008 Marko Tainio

# This model change the calculated 1-dimensional data to 2-dimensional so
# that data can be downloaded e.g. to Analytica more easily.
# Description of data:
# - Intake fraction data: iF for European population due to European emissions

#The location of data files
folder="N:/Huippuyksikko/Tutkimus/R19_Kopra/Mallit/IntakeFraction2/Result/European_emission/"

#Month/concentration - month 13 in annual mean
for (i in 1:13) {

#Y must be defined before for -loop
y = NULL

# Calculation for different countries. There are 48 sources in data.
for (j in 1:48) {

# This reads from file the titles of countries to be used in calculations
maa = read.table("N:/Huippuyksikko/Tutkimus/R19_Kopra/Mallit/IntakeFraction2/EUR_country.txt", sep="")
maa = maa[,j]

# This read the intake fraction data files
a = paste(folder,"2000/if_emis_eur_year_2000_country_",maa,"_month_",i,"_pop_all.txt", sep="")
a = read.table(a, sep="")

#This combines all the data to 1-dimensional data
y = c(y, a[,2])

}

#This change the 1-dimensional data to 2-dimensional
# 40 is the number of populations used, 48 is number of emission sources (countries)
dim(y) = c(40,48)

#Recording of results
write.csv(y, paste(folder, file="comb_if_year_2000_country_all_month_",i,".csv", sep=""))

}

#Remove of all data from memory
rm(folder, i, y, maa, a)

#21.4.2008, Marko Tainio

# This model calculates intake fractions (iF) for Finnish
# primary PM2.5 emissions. Description of model:
# - Emission source: Finnish primary PM2.5 emissions from different sectors
# - Dispersion: Dispersion over the Europe - small scale (11 countries)
# - Population: European population (11 countries)

#Work folders - the location of data and where the output is recorded
wrk_input = "N:/Huippuyksikko/Tutkimus/R19_Kopra/Data 2/"
wrk_output = "N:/Huippuyksikko/Tutkimus/R19_Kopra/Mallit/IntakeFraction2/Result/"
wrk_emis = "N:/Huippuyksikko/Tutkimus/R19_Kopra/Mallit/IntakeFraction2/"

#Breathing rate, 20 m^3/day -> unit g/s
BR = 20/(24 * 60 * 60)

#Calculation of different years (2000, 2001, 2002)
for (l in 0:2) {

#The calcultion of sectors - seven sectors are [DOM OTH TRA LPC LPP AGR ALL]
for (k in 1:7) {

#Emission volume - this part copy emission data in unit g/s
Q = paste(wrk_emis,"fin_emission_200",l,"_01.csv", sep="")
Q = read.csv(Q, header=T)

#Concentration data
#The data is read separately for different years because
#year 2000 and 2001,2002 data has different dimensions
rm(conc)
if (l == 0) (conc = paste(wrk_input,"Concentration_primary_sector_200",l,"/monthly_PM_2_5/Fin_200",l,"_PM_hirlam_v3_8_1_", sep=""))
if (l == 1) (conc = paste(wrk_input,"Concentration_primary_sector_200",l,"/monthly_200",l,"_PM_2_5/Fin_200",l,"_PM_EC_v3_8_1_", sep=""))
if (l == 2) (conc = paste(wrk_input,"Concentration_primary_sector_200",l,"/monthly_200",l,"_PM_2_5/Fin_200",l,"_PM_EC_v3_8_1_", sep=""))

#The reading order for the data files
if (k == 1) (conc = paste(conc,"_DOM_200",l,"12.txt", sep=""))
if (k == 2) (conc = paste(conc,"_OTH_200",l,"12.txt", sep=""))
if (k == 3) (conc = paste(conc,"_TRA_200",l,"12.txt", sep=""))
if (k == 4) (conc = paste(conc,"_LPC_200",l,"12.txt", sep=""))
if (k == 5) (conc = paste(conc,"_LPP_200",l,"12.txt", sep=""))
if (k == 6) (conc = paste(conc,"_AGR_200",l,"12.txt", sep=""))
if (k == 7) (conc = paste(conc,"all_srcs.txt", sep=""))
conc = read.table(conc, header=T, sep="")

#Population data
#The data is read separately for different years because
#year 2000 and 2001,2002 data has different dimensions
rm(pop)
if (l == 0) (pop = paste(wrk_input, "Population_Combined/Pop_fin00_all.txt", sep=""))
if (l == 1) (pop = paste(wrk_input, "Population_Combined/Pop_fin01_all.txt", sep=""))
if (l == 2) (pop = paste(wrk_input, "Population_Combined/Pop_fin01_all.txt", sep=""))
pop = read.table(pop, header=T, sep="")

#Co-ordination check-in - this check that co-ordinate system is same
#in population and concentration data.
check = (conc["lon"] - pop["Lo"]) + (conc["lat"] - pop["La"])
sum(check)

#This part slice the concentration data from the data file
#(that is, LO and LA data is sliced off)
month = c("jan", "feb", "mar", "apr", "may", "jun", "jul", "aug", "sep", "oct", "nov", "dec", "annual")
conc = data.frame(conc[month])

#Unit change in concentration data: kg/m^3 -> g/m^3
conc = conc * 1000

#This part slice the population data from the data file
#(that is, LO and LA data is sliced off)
#The data is read separately for different years because
#year 2000 and 2001,2002 data has different dimensions
rm(country_fin)
if (l == 0) (country_fin = c("Sum_ASUKKA", "Sum_Estoni", "Sum_Latvia", "Sum_Lithua", "Sum_Sweden", "Sum_Poland", "Sum_Denmar", "Sum_German", "Sum_Norway", "Sum_Belaru", "Sum_Russia", "Pop_all"))
if (l == 1) (country_fin = c("Finland", "Sum_Estoni", "Sum_Latvia", "Sum_Lithua", "Sum_Poland", "Sum_Sweden", "Sum_Denmar", "Sum_German", "Sum_Norway", "Sum_Belaru", "Sum_Russia", "Pop_all"))
if (l == 2) (country_fin = c("Finland", "Sum_Estoni", "Sum_Latvia", "Sum_Lithua", "Sum_Poland", "Sum_Sweden", "Sum_Denmar", "Sum_German", "Sum_Norway", "Sum_Belaru", "Sum_Russia", "Pop_all"))
pop = data.frame(pop[country_fin])

#X must be defined before for -loop
x = NULL

#Month/concentration - month 13 in annual mean
for (i in 1:13) {

#Removal of huge values
#There are some weird very high values in some concentration datasets.
#This code replace those values with value 0.
y = conc[,i]
y = ifelse(y > 100, 0, y)

#Country/population - iF's are calculated separately for 12 different population
for (j in 1:12) {

x[j] = sum(y * pop[,j])

}

#Calculation of iF
x = (x * BR) / Q[k,i]

#Names for the data
names(x) = c("FIN", "EST", "LAT", "LIT", "SWE", "POL", "DEN", "GER", "NOR", "BEL", "RUS", "All")

#Recording of results
write.table(x, paste(wrk_output, file = "Finnish_emission/200",l,"/if_emis_fin_year_200",l,"_sector_",k,"_month_",i,".txt", sep=""), row.names = TRUE, col.names = FALSE)

}

}

}

rm(wrk_input, wrk_output, wrk_emis, BR, l, k, Q, conc, pop, check, month, country_fin, x, i, y, j)

# 21.4.2008 Marko Tainio

# This model change the calculated 1-dimensional data to 2-dimensional so
# that data can be downloaded e.g. to Analytica more easily.
# Description of data:
# - Intake fraction data: iF for European population due to Finnish emissions

#The location of data files
folder="N:/Huippuyksikko/Tutkimus/R19_Kopra/Mallit/IntakeFraction2/Result/Finnish_emission/"

#Year (2000, 2001, 2002)
for (k in 0:2) {

#The calcultion of sectors - seven sectors are [DOM OTH TRA LPC LPP AGR ALL]
for (j in 1:7) {

#Y must be defined before for -loop
y = NULL

#Month/concentration - month 13 in annual mean
for (i in 1:13) {

# This read the intake fraction data files
a = paste(folder,"200",k,"/if_emis_fin_year_200",k,"_sector_",j,"_month_",i,".txt", sep="")
a = read.table(a, sep="")

#This combines all the data to 1-dimensional data
y = c(y, a[,2])

}
#This change the 1-dimensional data to 2-dimensional
# 12 is the number of populations used, 13 is month (13 in annual mean)
dim(y) = c(12,13)

#Recording of results
write.csv(y, paste(folder, file="comb_if_year_200",k,"_sector_",j,"_month_all.csv", sep=""))

}

}

rm(folder, k, j, y, i, a)

#7.3.2008 Marko Tainio

#This model summarise the population from different population datasets used
#in intake fraction calculation


#######
#Data location for all datasets
data_input = ("C:/Varsova/Kopra/Data 2/")
data_output = ("C:/Varsova/Muutokset/Result/")


#######
#Finnish emissions, year 2000 population data
pop_fin00_all = paste(data_input,"Population_Combined/Pop_fin00_all.txt", sep="")
a_fin00 = read.table(pop_fin00_all, sep="", header=T)
title_fin00 = c("Sum_ASUKKA", "Sum_Estoni", "Sum_Latvia", "Sum_Lithua", "Sum_Sweden", "Sum_Poland", "Sum_Denmar", "Sum_German", "Sum_Norway", "Sum_Belaru", "Sum_Russia", "Pop_all")
a_fin00 = data.frame(a_fin00[title_fin00])

a = NULL

for (i in 1:12) {

a[i] = sum(a_fin00[,i])

}

names(a) = c("FIN", "EST", "LAT", "LIT", "SWE", "POL", "DEN", "GER", "NOR", "BEL", "RUS", "All")
write.csv(a, paste(data_output, file="pop_fin00_sum.csv", sep=""))

rm(pop_fin00_all, a_fin00, title_fin00, a)


#######
#Finnish emissions, year 2001 population data
pop_fin01_all = paste(data_input,"Population_Combined/Pop_fin01_all.txt", sep="")
a_fin01 = read.table(pop_fin01_all, sep="", header=T)
title_fin01 = c("Finland", "Sum_Estoni", "Sum_Latvia", "Sum_Lithua", "Sum_Sweden", "Sum_Poland", "Sum_Denmar", "Sum_German", "Sum_Norway", "Sum_Belaru", "Sum_Russia", "Pop_all")
a_fin01 = data.frame(a_fin01[title_fin01])

a = NULL

for (i in 1:12) {

a[i] = sum(a_fin01[,i])

}

names(a) = c("FIN", "EST", "LAT", "LIT", "SWE", "POL", "DEN", "GER", "NOR", "BEL", "RUS", "All")
write.csv(a, paste(data_output, file="pop_fin01_sum.csv", sep=""))

rm(pop_fin01_all, a_fin01, title_fin01, a)


#######
#European emissions, year 2000 population data
pop_eur00_all = paste(data_input,"Population_Combined/Pop_Eur00_all.txt", sep="")
a_eur00 = read.table(pop_eur00_all, sep="", header=T)
title_eur00 = c("Sum_Austri", "Sum_Belgiu", "Sum_Cyprus", "Sum_Czech_", "Sum_Denmar", "Sum_Estoni", "Sum_France", "Sum_German", "Sum_Greece", "Sum_Hungar", "Sum_Irish_", "Sum_Italy", "Sum_Latvia", "Sum_Lithua", "Sum_Luxemb", "Sum_Malta", "Sum_Nether", "Sum_Poland", "Sum_Portug", "Sum_Romani", "Sum_Slovak", "Sum_Sloven", "Sum_Spain", "Sum_UK", "Sum_Sweden", "Sum_Bulgar", "Sum_Norway", "Sum_Belaru", "Sum_Russia", "Sum_Switze", "Sum_Ukrain", "Sum_Moldov", "Sum_Croati", "Sum_Serbia", "Sum_Albani", "Sum_Macedo", "Sum_Turkey", "pop_all")
a_eur00 = data.frame(a_eur00[title_eur00])

a = NULL

for (i in 1:38) {

a[i] = sum(a_eur00[,i])

}

names(a) = c("Sum_Austri", "Sum_Belgiu", "Sum_Cyprus", "Sum_Czech_", "Sum_Denmar", "Sum_Estoni", "Sum_France", "Sum_German", "Sum_Greece", "Sum_Hungar", "Sum_Irish_", "Sum_Italy", "Sum_Latvia", "Sum_Lithua", "Sum_Luxemb", "Sum_Malta", "Sum_Nether", "Sum_Poland", "Sum_Portug", "Sum_Romani", "Sum_Slovak", "Sum_Sloven", "Sum_Spain", "Sum_UK", "Sum_Sweden", "Sum_Bulgar", "Sum_Norway", "Sum_Belaru", "Sum_Russia", "Sum_Switze", "Sum_Ukrain", "Sum_Moldov", "Sum_Croati", "Sum_Serbia", "Sum_Albani", "Sum_Macedo", "Sum_Turkey", "pop_all")
write.csv(a, paste(data_output, file="pop_eur00_sum.csv", sep=""))

rm(pop_eur00_all, a_eur00, title_eur00, a)


#European emissions, year 2000 population data, Finnish
pop_eur00_fin = paste(data_input,"Population_Combined/Pop_Eur00_fin.txt", sep="")
a_eur00 = read.table(pop_eur00_fin, sep="", header=T)
a_eur00 = sum(a_eur00[,3])

write.csv(a_eur00, paste(data_output, file="pop_eur00_sum_fin.csv", sep=""))

rm(pop_eur00_fin, a_eur00)

#######

rm(data_input, data_output)