diff --git a/modules/35_transport/edge_esm/output.gms b/modules/35_transport/edge_esm/output.gms
new file mode 100644
index 0000000000000000000000000000000000000000..49c34856e554c3fdaabef32622f91e31189bb5c4
--- /dev/null
+++ b/modules/35_transport/edge_esm/output.gms
@@ -0,0 +1,10 @@
+*** |  (C) 2006-2019 Potsdam Institute for Climate Impact Research (PIK)
+*** |  authors, and contributors see CITATION.cff file. This file is part
+*** |  of REMIND and licensed under AGPL-3.0-or-later. Under Section 7 of
+*** |  AGPL-3.0, you are granted additional permissions described in the
+*** |  REMIND License Exception, version 1.0 (see LICENSE file).
+*** |  Contact: remind@pik-potsdam.de
+*** SOF ./modules/35_transport/edge_esm/output.gms
+
+Execute "Rscript EDGE_transport.R --reporting";
+*** EOF ./modules/35_transport/edge_esm/output.gms
diff --git a/modules/35_transport/edge_esm/realization.gms b/modules/35_transport/edge_esm/realization.gms
index 90ede08e76b4840a88170c17195d245e7210f4ad..6ada3a91d37df1332df5cb49b29919f86613e43a 100644
--- a/modules/35_transport/edge_esm/realization.gms
+++ b/modules/35_transport/edge_esm/realization.gms
@@ -25,5 +25,6 @@ $Ifi "%phase%" == "preloop" $include "./modules/35_transport/edge_esm/preloop.gm
 $Ifi "%phase%" == "bounds" $include "./modules/35_transport/edge_esm/bounds.gms"
 $Ifi "%phase%" == "presolve" $include "./modules/35_transport/edge_esm/presolve.gms"
 $Ifi "%phase%" == "postsolve" $include "./modules/35_transport/edge_esm/postsolve.gms"
+$Ifi "%phase%" == "output" $include "./modules/35_transport/edge_esm/output.gms"
 *######################## R SECTION END (PHASES) ###############################
 *** EOF ./modules/35_transport/edge_esm.gms
diff --git a/scripts/iterative/EDGE_transport.R b/scripts/iterative/EDGE_transport.R
index 887e99369871eb4fd539c8b466b6035f29c7428f..1a7f4d340a6ee1f0b9db5062ef46a1ea0f8fd334 100644
--- a/scripts/iterative/EDGE_transport.R
+++ b/scripts/iterative/EDGE_transport.R
@@ -1,19 +1,31 @@
-require(data.table)
-require(gdx)
-require(gdxdt)
-require(edgeTrpLib)
-require(rmndt)
+library(optparse)
+
+opt_parser = OptionParser(
+  description = "Coupled version of EDGE-T, to be run within a REMIND output folder.",
+  option_list = list(
+    make_option(
+      "--reporting", action="store_true", default=FALSE,
+      help="Store output files in subfolder EDGE-T")));
+opt = parse_args(opt_parser);
+
+library(data.table)
+library(gdx)
+library(gdxdt)
+library(edgeTrpLib)
+library(rmndt)
+library(moinput)
+
 ## use cached input data for speed purpose
-require(moinput)
 setConfig(forcecache=T)
 
+data_folder <- "EDGE-T"
 
 mapspath <- function(fname){
     file.path("../../modules/35_transport/edge_esm/input", fname)
 }
 
 datapath <- function(fname){
-    file.path("input_EDGE", fname)
+    file.path(data_folder, fname)
 }
 
 REMINDpath <- function(fname){
@@ -35,12 +47,9 @@ EDGE_scenario <- cfg$gms$cm_EDGEtr_scen
 
 EDGEscenarios <- fread("../../modules/35_transport/edge_esm/input/EDGEscenario_description.csv")[scenario_name == EDGE_scenario]
 
-merge_traccs <<- EDGEscenarios[options == "merge_traccs", switch]
-addvintages <<- EDGEscenarios[options == "addvintages", switch]
-inconvenience <<- EDGEscenarios[options == "inconvenience", switch]
-selfmarket_taxes <<- EDGEscenarios[options == "selfmarket_taxes", switch]
-selfmarket_policypush <<- EDGEscenarios[options == "selfmarket_policypush", switch]
-selfmarket_acceptancy <<- EDGEscenarios[options == "selfmarket_acceptancy", switch]
+inconvenience <- EDGEscenarios[options == "inconvenience", switch]
+selfmarket_policypush <- EDGEscenarios[options == "selfmarket_policypush", switch]
+selfmarket_acceptancy <- EDGEscenarios[options == "selfmarket_acceptancy", switch]
 
 if (EDGE_scenario == "Conservative_liquids") {
   techswitch <<- "Liquids"
@@ -53,22 +62,22 @@ if (EDGE_scenario == "Conservative_liquids") {
   exit()
 }
 
-endogeff <<- EDGEscenarios[options== "endogeff", switch]
-enhancedtech <<- EDGEscenarios[options== "enhancedtech", switch]
-rebates_febates <<- EDGEscenarios[options== "rebates_febates", switch] ##NB THEY ARE ONLY IN PSI! ONLY WORKING IN EUROPE
-savetmpinput <<- FALSE
-smartlifestyle <<- EDGEscenarios[options== "smartlifestyle", switch]
-
-
 
 REMIND2ISO_MAPPING <- fread(REMINDpath(cfg$regionmapping))[, .(iso = CountryCode, region = RegionCode)]
 EDGE2teESmap <- fread(mapspath("mapping_EDGE_REMIND_transport_categories.csv"))
 
 
 ## input data loading
-input_path = paste0("../../modules/35_transport/edge_esm/input/")
+input_folder = paste0("../../modules/35_transport/edge_esm/input/")
+
+if (length(list.files(path = data_folder, pattern = "RDS")) < 7) {
+  createRDS(input_folder, data_folder,
+            SSP_scenario = scenario,
+            EDGE_scenario = EDGE_scenario)
+}
+inputdata <- loadInputData(data_folder)
+
 
-inputdata = createRDS(input_path, SSP_scenario = scenario, EDGE_scenario = EDGE_scenario)
 vot_data = inputdata$vot_data
 sw_data = inputdata$sw_data
 inco_data = inputdata$inco_data
@@ -77,27 +86,21 @@ int_dat = inputdata$int_dat
 nonfuel_costs = inputdata$nonfuel_costs
 price_nonmot = inputdata$price_nonmot
 
-## add learning optional
-setlearning = TRUE
-## add optional vintages
-addvintages = TRUE
 ## optional average of prices
 average_prices = FALSE
-## inconvenience costs instead of preference factors
-inconvenience = TRUE
 
-if (setlearning | addvintages){
-  ES_demand = readREMINDdemand(gdx, REMIND2ISO_MAPPING, EDGE2teESmap, REMINDyears)
-  ## select from total demand only the passenger sm
-  ES_demand = ES_demand[sector == "trn_pass",]
-}
+
+ES_demand_all = readREMINDdemand(gdx, REMIND2ISO_MAPPING, EDGE2teESmap, REMINDyears)
+## select from total demand only the passenger sm
+ES_demand = ES_demand_all[sector == "trn_pass",]
+
 
 
-if (setlearning & file.exists("demand_previousiter.RDS")) {
+if (file.exists(datapath("demand_previousiter.RDS"))) {
   ## load previous iteration number of cars
-  demand_BEVtmp = readRDS("demand_BEV.RDS")
+  demand_BEVtmp = readRDS(datapath("demand_BEV.RDS"))
   ## load previous iteration demand
-  ES_demandpr = readRDS("demand_previousiter.RDS")
+  ES_demandpr = readRDS(datapath("demand_previousiter.RDS"))
   ## calculate non fuel costs and
   nonfuel_costs = applylearning(gdx,REMINDmapping,EDGE2teESmap, demand_BEVtmp, ES_demandpr)
   saveRDS(nonfuel_costs, "nonfuel_costs_learning.RDS")
@@ -111,6 +114,7 @@ REMIND_prices <- merge_prices(
     intensity_data = int_dat,
     nonfuel_costs = nonfuel_costs)
 
+
 ## save prices
 ## read last iteration count
 keys <- c("iso", "year", "technology", "vehicle_type")
@@ -122,13 +126,13 @@ iter <- as.vector(gdxrrw::rgdx(gdx, list(name="o_iterationNumber"))$val)
 REMIND_prices[, iternum := iter]
 
 ## save REMIND prices (before dampening)
-saveRDS(REMIND_prices, paste0("REMINDprices", iter, ".RDS"))
+saveRDS(REMIND_prices, datapath(paste0("REMINDprices", iter, ".RDS")))
 
 
 if(average_prices){
 
   if(max(unique(REMIND_prices$iternum)) >= 20 & max(unique(REMIND_prices$iternum)) <= 30){
-    old_prices <- readRDS(pfile)
+    old_prices <- readRDS(datapath(pfile))
     all_prices <- rbind(old_prices, REMIND_prices)
     setkeyv(all_prices, keys)
     ## apply moving avg
@@ -138,10 +142,10 @@ if(average_prices){
   }else{
     all_prices <- REMIND_prices
   }
-  saveRDS(all_prices, pfile)
+  saveRDS(all_prices, datapath(pfile))
 
   ## save REMIND prices (after dampening)
-  saveRDS(REMIND_prices,paste0("REMINDpricesDampened", iter, ".RDS"))
+  saveRDS(REMIND_prices, datapath(paste0("REMINDpricesDampened", iter, ".RDS")))
 
 }
 
@@ -157,7 +161,9 @@ if (inconvenience) {
     inco_data = inco_data,
     logit_params = logit_params,
     intensity_data = int_dat,
-    price_nonmot = price_nonmot)
+    price_nonmot = price_nonmot,
+    selfmarket_policypush = selfmarket_policypush,
+    selfmarket_acceptancy = selfmarket_acceptancy)
 
 } else{
 
@@ -177,47 +183,65 @@ shares <- logit_data[["share_list"]] ## shares of alternatives for each level of
 mj_km_data <- logit_data[["mj_km_data"]] ## energy intensity at a technology level
 prices <- logit_data[["prices_list"]] ## prices at each level of the logit function, 1990USD/pkm
 
-if(addvintages){
-  ## calculate vintages (new shares, prices, intensity)
-  vintages = calcVint(shares = shares,
-                      totdem_regr = ES_demand,
-                      prices = prices,
-                      mj_km_data = mj_km_data,
-                      years = REMINDyears)
-
-  shares$FV_shares = vintages[["shares"]]$FV_shares
-  prices = vintages[["prices"]]
-  mj_km_data = vintages[["mj_km_data"]]
-}
+
+## calculate vintages (new shares, prices, intensity)
+vintages = calcVint(shares = shares,
+                    totdem_regr = ES_demand,
+                    prices = prices,
+                    mj_km_data = mj_km_data,
+                    years = REMINDyears)
+
+shares$FV_shares = vintages[["shares"]]$FV_shares
+prices = vintages[["prices"]]
+mj_km_data = vintages[["mj_km_data"]]
+
 
 ## use logit to calculate shares and intensities (on tech level)
 EDGE2CESmap <- fread(mapspath("mapping_CESnodes_EDGE.csv"))
 
-shares_intensity_demand <- shares_intensity_and_demand(
-    logit_shares=shares,
-    MJ_km_base=mj_km_data,
-    EDGE2CESmap=EDGE2CESmap,
-    REMINDyears=REMINDyears,
-    scenario=scenario,
-    REMIND2ISO_MAPPING=REMIND2ISO_MAPPING)
 
-demByTech <- shares_intensity_demand[["demand"]] ##in [-]
-intensity <- shares_intensity_demand[["demandI"]] ##in million pkm/EJ
-norm_demand <- shares_intensity_demand$demandF_plot_pkm ## total demand is 1, required for costs
-
-if (setlearning) {
-  demand_BEV=calc_num_vehicles( norm_dem_BEV = norm_demand[technology == "BEV" & ## battery vehicles
-                                                           subsector_L1 == "trn_pass_road_LDV_4W", ## only 4wheelers
-                                                           c("iso", "year", "sector", "vehicle_type", "demand_F") ],
-                                ES_demand = ES_demand)
-
-  ## save number of vehicles for next iteration
-  saveRDS(demand_BEV, "demand_BEV.RDS")
-  ## save the demand for next iteration renaming the column
-  setnames(ES_demand, old ="demand", new = "demandpr")
-  saveRDS(ES_demand, "demand_previousiter.RDS")
+shares_int_dem <- shares_intensity_and_demand(
+  logit_shares=shares,
+  MJ_km_base=mj_km_data,
+  EDGE2CESmap=EDGE2CESmap,
+  REMINDyears=REMINDyears,
+  scenario=scenario,
+  REMIND2ISO_MAPPING=REMIND2ISO_MAPPING,
+  demand_input = if (opt$reporting) ES_demand_all)
+
+demByTech <- shares_int_dem[["demand"]] ##in [-]
+intensity <- shares_int_dem[["demandI"]] ##in million pkm/EJ
+norm_demand <- shares_int_dem[["demandF_plot_pkm"]] ## total demand is 1, required for costs
+
+
+if (opt$reporting) {
+  saveRDS(vintages[["vintcomp"]], file = datapath("vintcomp.RDS"))
+  saveRDS(vintages[["newcomp"]], file = datapath("newcomp.RDS"))
+  saveRDS(shares, file = datapath("shares.RDS"))
+  saveRDS(logit_data$EF_shares, file = datapath("EF_shares.RDS"))
+  saveRDS(logit_data$mj_km_data, file = datapath("mj_km_data.RDS"))
+  saveRDS(logit_data$inconv_cost, file=datapath("inco_costs.RDS"))
+  saveRDS(shares_int_dem$demandF_plot_EJ,
+          file=datapath("demandF_plot_EJ.RDS"))
+  saveRDS(shares_int_dem$demandF_plot_pkm,
+          datapath("demandF_plot_pkm.RDS"))
+  saveRDS(logit_data$annual_sales, file = datapath("annual_sales.RDS"))
+  quit()
 }
 
+demand_BEV=calc_num_vehicles(
+  norm_dem_BEV = norm_demand[
+    technology == "BEV" & ## battery vehicles
+    subsector_L1 == "trn_pass_road_LDV_4W", ## only 4wheelers
+    c("iso", "year", "sector", "vehicle_type", "demand_F") ],
+  ES_demand = ES_demand)
+
+## save number of vehicles for next iteration
+saveRDS(demand_BEV, datapath("demand_BEV.RDS"))
+## save the demand for next iteration renaming the column
+setnames(ES_demand, old ="demand", new = "demandpr")
+saveRDS(ES_demand, datapath("demand_previousiter.RDS"))
+
 
 ## use logit to calculate costs
 budget <- calculate_capCosts(
diff --git a/scripts/output/single/EDGETransportValidation.R b/scripts/output/single/EDGETransportValidation.R
index 927aadbfc68fab3a0d7738efa721e70521bf3df4..2fdfab5fda67f8d52cd162448da9c8530501d319 100644
--- a/scripts/output/single/EDGETransportValidation.R
+++ b/scripts/output/single/EDGETransportValidation.R
@@ -1,4 +1,19 @@
+# |  (C) 2006-2019 Potsdam Institute for Climate Impact Research (PIK)
+# |  authors, and contributors see CITATION.cff file. This file is part
+# |  of REMIND and licensed under AGPL-3.0-or-later. Under Section 7 of
+# |  AGPL-3.0, you are granted additional permissions described in the
+# |  REMIND License Exception, version 1.0 (see LICENSE file).
+# |  Contact: remind@pik-potsdam.de
+
 require(rmarkdown)
+require(lucode)
+
+if(!exists("source_include")) {
+  ## Define arguments that can be read from command line
+  readArgs("outputdir","gdx_name","gdx_ref_name")
+}
+
+load(file.path(outputdir, "config.Rdata"))
 
 ## run EDGE transport validation output if required
 if(cfg$gms$transport == "edge_esm"){
diff --git a/scripts/output/single/notebook_templates/EDGETransportReport.Rmd b/scripts/output/single/notebook_templates/EDGETransportReport.Rmd
index 20e7b74c78da0e3abc926c4181f50af9a0a0feee..4286f64d6baf73942208786c49d5af4eabf8df2c 100644
--- a/scripts/output/single/notebook_templates/EDGETransportReport.Rmd
+++ b/scripts/output/single/notebook_templates/EDGETransportReport.Rmd
@@ -14,645 +14,426 @@ require(remind)
 require(gdxdt)
 require(gdx)
 require(rmndt)
-require(data.table)
-
+require(magclass)
+require(quitte)
+require(ggpubr)
+require(gridExtra)
 require(edgeTrpLib)
 ```
 
 
 ```{r, echo=FALSE, warning=FALSE}
+iso_plot = "DEU"
+output_folder = "EDGE-T/"
 
-dem_shares <- list()
-intensity <- list()
-demand_km <- list()
-demand_ej <- list()
-sw_tech <- list()
-prices_FV <- list()
+cols <- c("NG" = "#d11141",
+          "Liquids" = "#8c8c8c",
+          "Hybrid Liquids" = "#ffc425",
+          "Hybrid Electric" = "#f37735",
+          "BEV" = "#00b159",
+          "FCEV" = "#00aedb")
 
 datapath <- function(fname){
-  file.path("./input_EDGE/", fname)
+  file.path("./EDGE-T/", fname)
 }
 
 mapspath <- function(fname, scenariopath=""){
     file.path("../../modules/35_transport/edge_esm/input", fname)
 }
 
+
 ## Load mappings
-EDGE2CESmap <- fread(mapspath("mapping_CESnodes_EDGE.csv"))
+REMIND2ISO_MAPPING <- fread("../../config/regionmappingH12.csv")[, .(iso = CountryCode, region = RegionCode)]
+EDGE2teESmap <- fread(mapspath("mapping_EDGE_REMIND_transport_categories.csv"))
 
-REMIND2ISO_MAPPING <- fread("../../config/regionmappingH12.csv")[, .(iso = CountryCode,
-                                                                         region = RegionCode)]
+## load input data from last EDGE run
+demand_km <- readRDS(datapath(fname = "demandF_plot_pkm.RDS")) ## detailed energy services demand, million km
+demand_ej <- readRDS(datapath(fname = "demandF_plot_EJ.RDS")) ## detailed final energy demand, EJ
+vintcomp <- readRDS(datapath(fname = "vintcomp.RDS"))
+newcomp <- readRDS(datapath(fname = "newcomp.RDS"))
+shares <- readRDS(datapath(fname = "shares.RDS"))
+inco_tech <- readRDS(datapath(fname = "inco_costs.RDS"))
+EF_shares <- readRDS(datapath(fname = "EF_shares.RDS"))
+annual_sales <- readRDS(datapath(fname = "annual_sales.RDS"))
+mj_km_data <- readRDS(datapath(fname = "mj_km_data.RDS"))
+
+name_mif = list.files(pattern = "REMIND_generic", full.names = F)
+name_mif = name_mif[!grepl("withoutPlu", name_mif)]
+miffile <- as.data.table(read.quitte(name_mif))
+```
 
-EDGE2teESmap <- fread(mapspath("mapping_EDGE_REMIND_transport_categories.csv"))
+# LDVs vintages
 
-years <- c(1990,
-           seq(2005, 2060, by = 5),
-           seq(2070, 2110, by = 10),
-           2130, 2150)
-
-REMINDyears <- c(1990,
-           seq(2005, 2060, by = 5),
-           seq(2070, 2110, by = 10),
-           2130, 2150)
-
-## include the paths to the scenarios you want to compare
-output_folders <- "./"
-
-for(output_folder in output_folders){
-  ## load gdx for fuel prices and demand
-  gdx = file.path(output_folder, "fulldata.gdx")
-  ## load policy scenario
-  load(file.path(output_folder, "config.Rdata"))
-  REMIND_scenario <- cfg$gms$cm_GDPscen
-  EDGE_scenario <- cfg$gms$cm_EDGEtr_scen
-  policy_scenario <- cfg$gms$c_expname
-  
-  scen <- paste0(REMIND_scenario, "-", EDGE_scenario, "-", policy_scenario)
-  
-  ## load demand
-  ES_demand = readREMINDdemand(gdx, REMIND2ISO_MAPPING, EDGE2teESmap, years)
-
-  ## load input data
-  int_dat <- readRDS(datapath("harmonized_intensities.RDS"))
-  nonfuel_costs <- readRDS(datapath("UCD_NEC_iso.RDS"))
-  sw_data <- readRDS(datapath("SW.RDS"))
-  vot_data <- readRDS(datapath("VOT_iso.RDS"))
-  logit_params <- readRDS(datapath("logit_exp.RDS"))
-  price_nonmot <- readRDS(datapath("price_nonmot.RDS"))
-  
-    ## FIXME: hotfix to make the (empty) vot_data$value_time_VS1 with the right column types. Probably there is another way to do that, did not look for it.
-  vot_data$value_time_VS1$iso = as.character(vot_data$value_time_VS1$iso)
-  vot_data$value_time_VS1$subsector_L1 = as.character(vot_data$value_time_VS1$subsector_L1)
-  vot_data$value_time_VS1$vehicle_type = as.character(vot_data$value_time_VS1$vehicle_type)
-  vot_data$value_time_VS1$year = as.numeric(vot_data$value_time_VS1$year)
-  vot_data$value_time_VS1$time_price = as.numeric(vot_data$value_time_VS1$time_price)
-  
-  ## calculate prices
-  REMIND_prices <- merge_prices(
-    gdx = gdx,
-    REMINDmapping = REMIND2ISO_MAPPING,
-    REMINDyears = REMINDyears,
-    intensity_data = int_dat,
-    nonfuel_costs = nonfuel_costs)
-
-  ## calculates logit
-  logit_data <- calculate_logit(
-    REMIND_prices[tot_price > 0],
-    REMIND2ISO_MAPPING,
-    vot_data = vot_data,
-    sw_data = sw_data,
-    logit_params = logit_params,
-    intensity_data = int_dat,
-    price_nonmot = price_nonmot)
-
-
-  shares <- logit_data[["share_list"]] ## shares of alternatives for each level of the logit function
-  mj_km_data <- logit_data[["mj_km_data"]] ## energy intensity at a technology level
-  prices_FV[[scen]] <- REMIND_prices[, EDGE_scenario := scen] ## prices at each level of the logit function, 1990USD/pkm
-
-  ## calculate energy intensity and FE demand at a REMIND-region level for the desired level of aggregation
-  res <- shares_intensity_and_demand(
-    logit_shares=shares,
-    MJ_km_base=mj_km_data,
-    REMIND2ISO_MAPPING=REMIND2ISO_MAPPING,
-    EDGE2CESmap=EDGE2CESmap,
-    REMINDyears=REMINDyears,
-    demand_input=ES_demand)
-
-  dem_shares[[scen]] <- res$demand[, EDGE_scenario := scen]
-  intensity[[scen]] <- res$demandI[, EDGE_scenario := scen]
-  demand_km[[scen]] <- res$demandF_plot_pkm[, EDGE_scenario := scen]
-  demand_ej[[scen]] <- res$demandF_plot_EJ[, EDGE_scenario := scen]
-  sw_tech[[scen]] <- sw_data$FV_final_SW[, EDGE_scenario := scen]
-}
+```{r, echo=FALSE, warning=FALSE}
 
-dem_shares <- rbindlist(dem_shares)
-intensity <- rbindlist(intensity)
-demand_km <- rbindlist(demand_km)
-demand_ej <- rbindlist(demand_ej)
-sw_tech <- rbindlist(sw_tech)
-prices_FV <- rbindlist(prices_FV)
-```
+plotVint = function(vintcomp, newcomp, sharesVS1){
+  vintcomp = vintcomp[,.(totdem, iso, subsector_L1, year, technology,vehicle_type, sector, sharetech_vint)]
+  newcomp = newcomp[,.(iso, subsector_L1, year, technology,vehicle_type, sector, sharetech_new)]
+
+  allfleet = merge(newcomp, vintcomp, all =TRUE, by = c("iso", "sector", "subsector_L1", "vehicle_type", "technology",  "year"))
+  allfleet = merge(allfleet, sharesVS1[,.(shareVS1 = share, iso, year, vehicle_type, subsector_L1)], all.x=TRUE, by = c("iso", "year", "vehicle_type", "subsector_L1"))
+  allfleet[,vintdem:=totdem*sharetech_vint*shareVS1]
+  allfleet[,newdem:=totdem*sharetech_new*shareVS1]
+  allfleet=melt(allfleet, id.vars = c("iso", "sector", "subsector_L1", "vehicle_type", "technology",
+                                      "year"), measure.vars = c("vintdem", "newdem"))
+  allfleet[,alpha:=ifelse(variable == "vintdem", 0, 1)]
+
+  load_factor = 2
+  annual_mileage = 15000
+  allfleet = allfleet[,.(value = sum(value/load_factor/annual_mileage)), by = c("iso", "technology", "variable", "year")]
+
+  allfleet = merge(allfleet, REMIND2ISO_MAPPING, by = "iso")
+  allfleet = allfleet[,.(value = sum(value)), by = c("region", "technology", "variable", "year")]
+  allfleet[,alphaval := ifelse(variable =="vintdem", 1,0)]
+
+
+  p = ggplot()+
+  geom_bar(data = allfleet[year %in% c(2015,2030,2050)],
+  aes(x=as.character(year),y=value, group=interaction(variable, technology),
+  fill = technology), alpha = 0.5, position="stack", stat = "identity", width = 0.5)+
+  geom_bar(data = allfleet[year %in% c(2015,2030,2050)],
+  aes(x=as.character(year),y=value, group=interaction(variable, technology),
+  fill = technology, alpha = factor(alphaval)), position="stack", stat = "identity", width = 0.5, color = "black")+
+  guides(fill = guide_legend(reverse=TRUE))+
+  theme_minimal()+
+    facet_wrap(~region, nrow = 4)+
+  theme(axis.text.x = element_text(angle = 90, hjust = 1),
+        axis.text = element_text(size=7),
+        title = element_text(size=8),
+        legend.text = element_text(size=8))+
+      scale_x_discrete(breaks = c(2015,2030,2050))+
+  scale_alpha_discrete(breaks = c(1,0), name = "Status", labels = c("Vintages","New additions")) +
+  guides(linetype=FALSE,
+         fill=guide_legend(reverse=FALSE, title="Transport mode"))+
+  scale_fill_manual(values = cols)+
+  labs(y = "LDV fleet [million Veh]", x="")
 
+  return(p)
+  }
 
 
-```{r, echo=FALSE}
-## plot settings
-years_plot = c(2010,2015,2020,2025,2030,2040,2050) ## in bar charts, these are the time steps that are represented
-year_single = 2050
-region_plot = "NEU" ## in case is a region specific plot, this region is represented
-sector_plot ="trn_pass" ## in case is a sector specific plot, this sector is represented
+p = plotVint(vintcomp, newcomp, shares$VS1_shares)
 
+p
 
-##conversion rate 2005->1990 USD
-CONV_2005USD_1990USD=0.67
+```
 
-# print(paste0("Scenario: ", REMIND_scenario))
-print(paste0("Regional plots are about ",region_plot))
-print(paste0("Sectoral plots are about ",sector_plot))
+# Inconvenience cost trend
 
-## maps
-cesmap <- data.table(CES_parent=c("_p_sm", "_p_lo", "_f_lo", "_f_sm"),
-                     CES_label=c("Passenger, Short-to-Medium Distances",
-                                 "Passenger, Long Distances",
-                                 "Freight, Long Distances",
-                                 "Freight, Short-to-Medium Distances"))
+```{r, echo=FALSE, warning=FALSE}
 
-EDGE_sectormap <- data.table(sector=c("trn_pass", "trn_freight", "trn_aviation_intl", "trn_shipping_intl"),
-                     CES_label=c("Passenger, Short-to-Medium Distances",
-                                 "Passenger, Long Distances",
-                                 "Freight, Long Distances",
-                                 "Freight, Short-to-Medium Distances"))
+p=ggplot()+
+  geom_bar(data = inco_tech[iso == iso_plot & subsector_L1 == "trn_pass_road_LDV_4W" & vehicle_type == "Large Car and SUV" & year<=2100 & year>=2010], aes(x = as.character(year), y = pinco, group = technology, fill = technology), position = position_stack(), stat = "identity")+
+  facet_grid(~technology)+
+    theme_minimal()+
+  scale_fill_manual(values = cols)+
+  expand_limits(y = c(0,0.8))+
+      scale_x_discrete(breaks = c(2015,2050,2100))+
+  theme(axis.text.x = element_text(angle = 90, vjust = +0.1),
+        legend.position = "none",
+        strip.background = element_rect(color = "grey"))+
+  labs(x = "", y = "Inconvenience cost [$/pkm]", title = paste0("Example of ", iso_plot))
+
+p
 
 ```
 
+# Endogenous intensity for Liquids
 
-```{r, echo=FALSE}
-## aggregate demands to REMIND regions
-demandF_plot_EJ <- demand_ej[,c("EDGE_scenario", "sector","subsector_L3","subsector_L2",
-                                "subsector_L1","vehicle_type","technology", "iso","year","demand_EJ")]
-demandF_plot_pkm <- demand_km[,c("EDGE_scenario", "sector","subsector_L3","subsector_L2",
-                                 "subsector_L1","vehicle_type","technology","iso","year","demand_F")]
-
-demandF_plot_EJ=aggregate_dt(demandF_plot_EJ,REMIND2ISO_MAPPING,
-                             datacols = c("EDGE_scenario", "sector", "subsector_L3", "subsector_L2", "subsector_L1",
-                                          "vehicle_type", "technology"),
-                             valuecol = "demand_EJ")
-
-demandF_plot_pkm=aggregate_dt(demandF_plot_pkm,REMIND2ISO_MAPPING,
-                              datacol = c("EDGE_scenario", "sector","subsector_L3","subsector_L2",
-                                          "subsector_L1","vehicle_type","technology"),
-                              valuecol = "demand_F")
-## add GLO
-glo <- demandF_plot_pkm[,.(region="GLO", demand_F=sum(demand_F)),
-                        by=eval(names(demandF_plot_pkm)[2:9])]
-demandF_plot_pkm <- rbind(demandF_plot_pkm, glo)
-glo <- demandF_plot_EJ[,.(region="GLO", demand_EJ=sum(demand_EJ)),
-                        by=eval(names(demandF_plot_EJ)[2:9])]
-demandF_plot_EJ <- rbind(demandF_plot_EJ, glo)
+```{r, echo=FALSE, message=FALSE, warning=FALSE}
+## Choice of the energy intensity (of the new sales)
+intcompplotf = function(EF_shares, FV_shares, VS1_shares){
+  EF_shares = EF_shares[,c("iso", "year", "technology", "vehicle_type", "subsector_L1", "subsector_L2", "subsector_L3", "sector", "share","type")]
+  setnames(EF_shares, old="share", new = "shareINT")
+  FV_shares = FV_shares[iso == iso_plot & subsector_L1 == "trn_pass_road_LDV_4W" & technology == "Liquids"]
+  setnames(FV_shares, old="share", new = "shareF")
+  VS1_shares = VS1_shares[iso == iso_plot & subsector_L1 == "trn_pass_road_LDV_4W"]
+  shares_LDV = merge(FV_shares, EF_shares, all = FALSE, by = c("iso", "year", "technology", "vehicle_type", "subsector_L1"))
+
+  shares_LDV[, shareIF := shareF*shareINT]
+  shares_LDV <- shares_LDV[,.(shareIF=sum(shareIF)),by=c("iso","technology","type","vehicle_type","subsector_L1", "year")]
+
+  shares_LDV = merge(shares_LDV, VS1_shares, all = TRUE, by = c("iso", "year", "vehicle_type", "subsector_L1"))
+  shares_LDV[, shareIS1 := shareIF*share]
+  shares_LDV <- shares_LDV[,.(shareIS1=sum(shareIS1)),by=c("iso","type", "technology","subsector_L1","year")]
+  p = ggplot()+
+    geom_bar(data = shares_LDV[year<=2100  & year>=2025], aes(x=year,y=shareIS1, group = technology, fill = technology), alpha = 0.5, position = position_fill(), stat = "identity")+
+    geom_bar(data = shares_LDV[year<=2100  & year>=2025], aes(x=year,y=shareIS1, group = technology, fill = technology, alpha = type), position = position_fill(), stat = "identity")+
+    facet_wrap("technology")+
+    theme_minimal()+
+    expand_limits(y = c(0,1))+
+    scale_fill_manual("technology", values = cols)+
+    scale_alpha_discrete("Type")+
+    labs(y = "Share [%]", title = paste0("Energy intensity new sales of Liquids, example for ", iso_plot))
 
-```
+  return(p)
 
-## ES
-
-
-```{r, echo=FALSE}
-##chunk of code that plots the ES
-ES_modes_bar1=function(demandpkm){
-  #group by subsector_L3 and summarise the demand
-  df=demandpkm[, .(demand_F=sum(demand_F)
-                   ), by = c("EDGE_scenario", "region", "year","sector","subsector_L1")]
-  df[,demand_F:=demand_F   ## in millionkm
-     *1e-6                 ## in trillion km
-     ]
-  df=df[order(year)]
-  # #filter only 2020, 2050 and 2100
-  df=df[year %in% years_plot,]
-  #separate into passenger and freight categories
-  pass=c("trn_pass","trn_aviation_intl")
-  freight=c("trn_freight","trn_shipping_intl")
-  ## give proper names to the categories  
-  df=merge(df, edgeTrpLib::L1mapping, all.x=TRUE, by="subsector_L1")
-  #plot
-  plot_p=ggplot()+
-    geom_bar(data=df%>%filter(sector %in% pass, year %in% years_plot, region == region_plot),
-             aes(x=year,y=demand_F,group=mode,fill=mode),position=position_stack(),stat="identity")+
-    facet_wrap(~EDGE_scenario)+
-    ggtitle("Energy Services Demand - Passenger Transport Modes")+
-    theme_light()+
-    theme(axis.text.x = element_text(angle = 90, hjust = 1))+
-    scale_x_continuous(breaks=years_plot)+
-    xlab("Year")+
-    ylab("Energy Services Demand (trillion pkm)")+ 
-    guides(fill=guide_legend(title="Transport mode"))+
-    theme(axis.text.x = element_text(angle = 90, hjust = 1),
-          axis.text = element_text(size=13),
-          title = element_text(size=13),
-          legend.text = element_text(size=13))
-  
-  plot_f=ggplot()+
-    geom_bar(data=df%>%filter(sector %in% freight,year %in% years_plot, region == region_plot),
-             aes(x=year,y=demand_F,group=mode,fill=mode),position=position_stack(),stat="identity")+
-    facet_wrap(~EDGE_scenario)+
-    ggtitle("Energy Services Demand - Freight Transport Modes")+
-    theme_light()+
-    theme(axis.text.x = element_text(angle = 90, hjust = 1))+
-    scale_x_continuous(breaks=years_plot)+
-    xlab("Year")+
-    ylab("Energy Services Demand (trillion tkm)")+ 
-    guides(fill=guide_legend(title="Transport mode"))+
-    theme(axis.text.x = element_text(angle = 90, hjust = 1),
-          axis.text = element_text(size=13),
-          title = element_text(size=13),
-          legend.text = element_text(size=13))
-  
-  plot=list(plot_p,plot_f)
-  return(plot)
 }
 
-p=ES_modes_bar1(demandpkm=demandF_plot_pkm)
-p[[1]]
-p[[2]]
-
+intcompplotf(EF_shares, shares$FV_shares, shares$VS1_shares)
 ```
 
-```{r, echo=FALSE}
-##chunk of code that plots the ES
-ES_modes_bar=function(demandpkm){
-  demandpkm[technology == "LA-BEV", technology := "BEV"]
-  ## use proper non-fuel mode names
-  demandpkm[technology %in% c("Cycle_tmp_technology", "Walk_tmp_technology"), technology := "Human Powered"]
-  #group by subsector_L3 and summarise the demand
-  df=demandpkm[, .(demand_F=sum(demand_F)),
-               by = c("EDGE_scenario", "region", "year","sector","subsector_L1", "technology")]
-  df[,demand_F:=demand_F   ## in millionkm
-     *1e-6                 ## in trillion km
-     ]
-  df=df[order(year)]
-  #separate into passenger and freight categories
-  pass=c("trn_pass","trn_aviation_intl")
-  freight=c("trn_freight","trn_shipping_intl")
-  ## give proper names to the categories
-  df=merge(df, edgeTrpLib::L1mapping, all.x=TRUE, by="subsector_L1")
-  #plot
-  plot_p=ggplot()+
-    geom_bar(data=df%>%filter(sector %in% pass,year %in% years_plot,region==region_plot),
-             aes(x=year,y=demand_F,group=technology,fill=technology),
-             position=position_stack(),stat="identity")+
-    ggtitle("Energy Services Demand by Technology, Passenger Transport, EUR")+
-    theme_light()+
-    facet_wrap(~EDGE_scenario)+
-    scale_x_continuous(breaks=years_plot)+
-    xlab("Year")+
-    ylab("Energy Services Demand (trillion pkm)")+ 
-    guides(fill=guide_legend(title="Technology"))+
-    theme(axis.text.x = element_text(angle = 90, hjust = 1),
-          axis.text = element_text(size=13),
-          title = element_text(size=13),
-          legend.text = element_text(size=13))
-  
-  plot_f=ggplot()+
-    geom_bar(data=df%>%filter(sector %in% freight,year %in% years_plot, region==region_plot),
-             aes(x=year,y=demand_F, group=technology, fill=technology),
-             position=position_stack(),stat="identity")+
-    ggtitle(paste0("Energy Services Demand by Technology, Freight Transport, ", region_plot))+
-    facet_wrap(~EDGE_scenario) +
-    theme_light()+
-    scale_x_continuous(breaks=years_plot)+
-    xlab("Year")+
-    ylab("Energy Services Demand (trllion tkm)")+ 
-    guides(fill=guide_legend(title="Technology"))+
-    theme(axis.text.x = element_text(angle = 90, hjust = 1),
-          axis.text = element_text(size=13),
-          title = element_text(size=13),
-          legend.text = element_text(size=13))
-  
-  plot=list(plot_p,plot_f)
-  return(plot)
+# Sales of LDVs
+
+```{r, echo=FALSE, warning=FALSE}
+salesplot = function(annual_sales){
+  annual_sales = unique(annual_sales[,c("iso","year", "technology", "shareFS1")])
+  annual_sales <- annual_sales[,.(shareFS1=sum(shareFS1)),by=c("iso","technology","year")]
+
+  p = ggplot()+
+    geom_bar(data = annual_sales[year<=2050  & year>=2015 & iso == iso_plot], aes(x=as.numeric(as.character(year)),y=shareFS1, group = technology, fill = technology), position = position_stack(), stat = "identity")+
+    theme_minimal()+
+    scale_fill_manual("Technology", values = cols)+
+    expand_limits(y = c(0,1))+
+    scale_x_continuous(breaks = c(2015,2030,2050))+
+    theme(axis.text.x = element_text(angle = 90, vjust = +0.1),
+          strip.background = element_rect(color = "grey"),
+          legend.position = "none")+
+    labs(x = "", y = "Market share on LDVs [%]", title = paste0("Sales composition, example of ", iso_plot))
+
+  return(p)
 }
 
-p=ES_modes_bar(demandpkm=demandF_plot_pkm)
-p[[1]]
-p[[2]]
 
+salesplot(annual_sales)
 ```
 
-```{r, echo=FALSE}
-## plot ES for LDVs only divided by fuel
-ES_modes_LDV_bar=function(demandpkm){
-  demandpkm[technology == "LA-BEV", technology := "BEV"]
-  ## use proper non-fuel mode names
-  demandpkm[technology %in% c("Cycle_tmp_technology", "Walk_tmp_technology"), technology := "Human Powered"]
-  #group by subsector_L3 and summarise the demand
-  df=demandpkm[, .(demand_F=sum(demand_F)),
-               by = c("EDGE_scenario", "region", "year","sector","subsector_L1", "technology")]
-  df[,demand_F:=demand_F   ## in millionkm
-     *1e-6                 ## in trillion km
-     ]
-  df=df[order(year)]
-  ## give proper names to the categories
-  df=merge(df, edgeTrpLib::L1mapping, all.x=TRUE, by="subsector_L1")
-  
-  
-  ## select order of facets
-  df$technology = factor(df$technology, levels=c("Liquids","Hybrid Liquids","NG","BEV","FCEV"))
-  
-  #plot
-  plot_LDV=ggplot()+
-    geom_bar(data=df%>%filter(year %in% years_plot,region==region_plot, mode %in% c("4W","2W")),
-             aes(x=year,y=demand_F,group=technology,fill=technology), alpha = 0.9,
-             position=position_stack(),stat="identity")+
-    ggtitle(paste0("Energy Services Demand by Technology, LDVs", region_plot))+
-    theme_light()+
-    facet_wrap(~EDGE_scenario)+
-    theme(axis.text.x = element_text(angle = 90, hjust = 1))+
-    xlab("Year")+
-    ylab("Energy Services Demand (trillion pkm)")+ 
-    guides(fill=guide_legend(title="Technology"))+
-    theme(axis.text.x = element_text(angle = 90, hjust = 1),
-          axis.text = element_text(size=13),
-          title = element_text(size=13),
-          legend.text = element_text(size=13),
-          strip.text.x = element_text(size = 13, color = "black"),
-          strip.background=element_rect(fill="white"))+
-    scale_x_continuous(breaks=years_plot)+
-    scale_fill_brewer(palette = "Set1")
-    
-  
-  return(plot_LDV)
-}
+# Final energy demand
 
-p=ES_modes_LDV_bar(demandpkm=demandF_plot_pkm)
+```{r, echo=FALSE, warning=FALSE}
+demandEJplotf = function(demandEJ){
+  ## EDGE results
+  demandEJ <- demandEJ[, c("sector", "subsector_L3", "subsector_L2", "subsector_L1", "vehicle_type", "technology", "iso", "year", "demand_EJ")]
+
+  ## attribute aggregated mode and vehicle names for plotting purposes, and aggregate
+  demandEJ[, aggr_mode := ifelse(subsector_L2 == "trn_pass_road_LDV", "LDV", NA)]
+  demandEJ[, aggr_mode := ifelse(subsector_L3 %in% c("Passenger Rail", "HSR", "International Aviation", "Domestic Aviation"), "Pass non LDV", aggr_mode)]
+  demandEJ[, aggr_mode := ifelse(subsector_L2 %in% c("trn_pass_road_bus", "Bus"), "Pass non LDV", aggr_mode)]
+  demandEJ[, aggr_mode := ifelse(is.na(aggr_mode), "Freight", aggr_mode)]
+  demandEJ[, veh := ifelse(vehicle_type %in% c("Truck (0-1t)", "Truck (0-3.5t)", "Truck (0-2.7t)", "Truck (0-2t)"), "Trucks (<3.5t)", NA)]
+  demandEJ[, veh := ifelse(vehicle_type %in% c("Truck (16-32t)", "Truck (3.5-16t)", "Truck (6-15t)", "Truck (4.5-12t)", "Truck (2.7-4.5t)", "Truck (4.5-15t)"), "Trucks (3.5t-16)", veh)]
+  demandEJ[, veh := ifelse(vehicle_type %in% c("Truck (>15t)", "Truck (16-32t)", "Truck (>32t)" ), "Trucks (>16)", veh)]
+  demandEJ[, veh := ifelse(grepl("Large|SUV|Midsize|Multipurpose Vehicle|Van|3W Rural", vehicle_type), "Large Cars", veh)]
+  demandEJ[, veh := ifelse(grepl("Subcompact|Compact|Mini|Three-Wheeler", vehicle_type), "Small Cars", veh)]
+  demandEJ[, veh := ifelse(grepl("Motorcycle|Moped|Scooter", vehicle_type), "Motorbikes", veh)]
+  demandEJ[, veh := ifelse(grepl("bus|Bus", vehicle_type), "Bus", veh)]
+  demandEJ[, veh := ifelse(grepl("Freight Rail_tmp_vehicletype", vehicle_type), "Freight Rail", veh)]
+  demandEJ[, veh := ifelse(grepl("Passenger Rail|HSR", vehicle_type), "Passenger Rail", veh)]
+  demandEJ[, veh := ifelse(subsector_L3 == "Domestic Ship", "Domestic Shipping", veh)]
+  demandEJ[, veh := ifelse(subsector_L3 == "International Ship", "International Shipping", veh)]
+  demandEJ[, veh := ifelse(subsector_L3 == "Domestic Aviation", subsector_L3, veh)]
+  demandEJ[, veh := ifelse(subsector_L3 == "International Aviation", subsector_L3, veh)]
+  demandEJ[, veh := ifelse(is.na(veh), vehicle_type, veh)]
+  demandEJ = demandEJ[,.(demand_EJ = sum(demand_EJ)), by = c("iso", "year", "aggr_mode", "veh")]
+
+  demandEJ[, vehicle_type_plot := factor(veh, levels = c("LDV","Freight Rail", "Trucks (<3.5t)", "Trucks (3.5t-16)", "Truck (>12t)", "Trucks (>16)", "Trucks","Domestic Shipping", "International Shipping",
+                                                                "Motorbikes", "Small Cars", "Large Cars", "Van",
+                                                                "Domestic Aviation", "International Aviation", "Bus", "Passenger Rail",
+                                                                "Freight", "Freight (Inland)", "Pass non LDV", "Pass non LDV (Domestic)"))]
+
+  legend_ord <- c("Freight Rail", "Trucks (<3.5t)", "Trucks (3.5t-16)", "Truck (>12t)", "Trucks (>16)", "International Shipping","Domestic Shipping",  "Trucks",
+                  "Motorbikes", "Small Cars", "Large Cars", "Van",
+                  "International Aviation", "Domestic Aviation","Bus", "Passenger Rail",
+                  "Freight", "LDV", "Pass non LDV", "Freight (Inland)", "Pass non LDV (Domestic)")
+
+  demandEJ = merge(demandEJ, REMIND2ISO_MAPPING, by = "iso")
+  demandEJ = demandEJ[,.(demand_EJ= sum(demand_EJ)), by = c("region", "year", "vehicle_type_plot", "aggr_mode")]
 
-p
+  p=ggplot()+
+    geom_area(data = demandEJ[year > 2010], aes(x=year, y=demand_EJ, group = interaction(vehicle_type_plot,aggr_mode), fill = vehicle_type_plot), color = "black", position= position_stack())+
+    facet_wrap(~region, nrow = 4)
+    labs(x = "", y = "Final Energy demand [EJ]")+
+    theme_minimal()+
+    # scale_fill_manual("Vehicle Type",values = cols, breaks=legend_ord)+
+    theme(axis.text.x = element_text(size = 8),
+          axis.text.y = element_text(size=8),
+          axis.title = element_text(size = 9),
+          title = element_text(size = 9),
+          legend.text = element_text(size = 9),
+          legend.title = element_text(size =9),
+          strip.text = element_text(size=9))
 
-```
 
-## FE
-
-```{r, echo=FALSE}
-FE_modes_bar=function(demandEJ){
-  #group by subsector_L1 and summarise the demand
-  df=demandEJ[, .(demand_EJ=sum(demand_EJ)),
-              by = c("EDGE_scenario", "region", "year","subsector_L1","subsector_L3","sector")]
-  df=df[order(year)]
-  df=df[year %in% years_plot,]
-  #separate into passenger and freight categories
-  pass=c("trn_pass","trn_aviation_intl")
-  freight=c("trn_freight","trn_shipping_intl")
-  ## give proper names to the categories
-  df <- merge(df, edgeTrpLib::L1mapping, all.x=TRUE, by="subsector_L1")
-  #plot
-  plot_p=ggplot()+
-    geom_bar(data=df%>%filter(sector %in% pass, region==region_plot),
-             aes(x=year,y=demand_EJ,group=mode,fill=mode),position=position_stack(),stat="identity",color="black")+
-    facet_wrap(~EDGE_scenario)+
-    ylab("Energy (EJ)") +
-    ggtitle(paste0("Final Energy Demand by Mode, Passenger, ", region_plot))+
-    theme(axis.text.x = element_text(angle = 90))+
-    scale_x_continuous(breaks=years_plot)
-  
-  plot_f=ggplot()+
-    geom_bar(data=df%>%filter(sector %in% freight,region==region_plot),
-             aes(x=year,y=demand_EJ,group=mode,fill=mode),position=position_stack(),stat="identity",color="black")+
-    facet_wrap(~EDGE_scenario) +
-    ggtitle(paste0("Final Energy Demand by Mode, Freight, ", region_plot))+
-    ylab("Energy (EJ)") +
-    theme(axis.text.x = element_text(angle = 90))+
-    scale_x_continuous(breaks=years_plot)
-  
-  plot=list(plot_p,plot_f)
-  return(plot)
+  return(p)
 }
 
-p=FE_modes_bar(demandEJ = demandF_plot_EJ)
-p[[1]]
-p[[2]]
+## Final Energy demand
+demandEJplotf(demand_ej)
 ```
 
-```{r, echo=FALSE}
-## function that calculates FE split and splits out the liquids by source
-FE_modes_bar_oilcomponent=function(demandEJ, msect="trn_pass", region_plot){
-  if(msect == "trn_pass")
-    sector_display = "Passenger"
-  if(msect == "trn_freight")
-    sector_display = "Freight"
-    
-  #group by subsector_L1 and summarise the demand
-  df=demandEJ[, .(demand_EJ=sum(demand_EJ)),
-              by = c("EDGE_scenario", "region", "year","subsector_L1","subsector_L3","sector","technology")]
-  df=df[order(year)]
-  df=df[year %in% years_plot,]
-  ## give proper names to the categories
-  df=merge(df, edgeTrpLib::L1mapping, all.x=TRUE, by="subsector_L1")
-  
-  df[,technology := ifelse(technology == "LA-BEV", "BEV", technology)]
-  df[,technology := ifelse(technology == "Electric", "El. Trains", technology)]
-  ## select order of facerts
-  df$technology = factor(df$technology, levels=c("BEV","FCEV","Hybrid Liquids", "El. Trains", "NG","Liquids", "Coal"))
-
-  #plot
-  plot_psm = ggplot()+
-    geom_bar(data=df%>%filter(sector == msect, region == region_plot),
-             aes(x=year,y=demand_EJ,group=technology,fill=technology),
-             position=position_stack(),stat="identity", alpha = 0.9)+
-    theme_light()+
-    ggtitle(paste0("Final Energy Demand by Tech, ", sector_display, ", ", region_plot))+
-    theme(axis.text.x = element_text(angle = 90),
-        strip.text.x = element_text(size = 13, color = "black"))+
-    scale_x_continuous(breaks=years_plot)+
-    scale_fill_brewer(palette = "Set1")+
-    xlab("Year")+
-    ylab("Final energy demand [EJ]")+ 
-    facet_wrap(~EDGE_scenario) +
-    guides(fill=guide_legend(title="Technology"))
-  
-    plot_psm_LDV = ggplot()+
-    geom_bar(data=df%>%filter(sector == msect, region == region_plot, mode == "4W"),
-             aes(x=year,y=demand_EJ,group=technology,fill=technology),
-             position=position_stack(),stat="identity", alpha = 0.9)+
-    theme_light()+
-    ggtitle(paste0("Final Energy Demand by Tech, LDVs, ", region_plot))+
-    theme(axis.text.x = element_text(angle = 90),
-          strip.text.x = element_text(size = 13, color = "black"),
-          strip.background=element_rect(fill="white"),
-          axis.text = element_text(size=13),
-          title = element_text(size=13),
-          legend.text = element_text(size=13))+
-    scale_x_continuous(breaks=years_plot)+
-    scale_fill_brewer(palette = "Set1")+
-    xlab("Year")+
-    ylab("Final energy demand [EJ]")+ 
-    facet_wrap(~EDGE_scenario) +
-    guides(fill=guide_legend(title="Technology"))  
-  
-    plot_list = list(plot_psm, plot_psm_LDV)
-  return(plot_list)
-}
+# LDVs final energy demand
 
+```{r, echo=FALSE, warning=FALSE}
+## demand EJ for LDV, divided by fuel type
+
+demandEJLDVplotf <- function(demandEJ){
+  demandEJ = demandEJ[subsector_L1 == "trn_pass_road_LDV_4W",]
+  demandEJ <- demandEJ[, c("sector", "subsector_L3", "subsector_L2", "subsector_L1", "vehicle_type", "technology", "iso", "year", "demand_EJ")]
+
+  demandEJ = merge(demandEJ, REMIND2ISO_MAPPING, by = "iso")
+  demandEJ[technology == "Hybrid Liquids", technology := "Liquids"]
+  demandEJ[technology == "FCEV", technology := "Hydrogen"]
+  demandEJ[technology == "BEV", technology := "Electricity"]
+  demandEJ = demandEJ[, .(demand_EJ = sum(demand_EJ)), by = c("region", "year", "technology")]
+
+    p = ggplot()+
+    geom_area(data = demandEJ[year > 2010], aes(x=year, y=demand_EJ, group = technology, fill = technology), color="black",position= position_stack())+
+    labs(x = "", y = "Final energy demand for LDVs [EJ]")+
+    facet_wrap(~region, nrow = 4)
+    theme_minimal()+
+    # scale_fill_manual("Vehicle Type",values = cols, breaks=legend_ord)+
+    theme(axis.text.x = element_text(size = 7),
+          axis.text.y = element_text(size=7),
+          axis.title = element_text(size = 8),
+          title = element_text(size = 8),
+          legend.text = element_text(size = 8),
+          legend.title = element_text(size = 8),
+          strip.text = element_text(size=8))
+
+    return(p)
 
+}
 
-FE_modes_bar_oilcomponent(demandEJ = demandF_plot_EJ, msect="trn_pass", region=region_plot)
-FE_modes_bar_oilcomponent(demandEJ = demandF_plot_EJ, msect="trn_pass", region="GLO")
 
+demandEJLDVplotf(demand_ej)
 ```
 
+# Energy services demand
+
+```{r, echo=FALSE, warning=FALSE}
+demandpkmplotf = function(demandpkm){
+  ## REMIND-EDGE results
+  demandpkm <- demandpkm[,c("sector","subsector_L3","subsector_L2",
+                            "subsector_L1","vehicle_type","technology", "iso","year","demand_F")]
+  demandpkm[,demand_F:=demand_F   ## in millionkm
+            *1e-6                      ## in trillion km
+            ]
+  ## attribute aggregated mode and vehicle names for plotting purposes, and aggregate
+  demandpkm[, aggr_mode := ifelse(subsector_L1 %in% c("Three-Wheeler", "trn_pass_road_LDV_4W"), "LDV", NA)]
+  demandpkm[, aggr_mode := ifelse(sector %in% c("trn_freight", "trn_shipping_intl"), "Freight", aggr_mode)]
+  demandpkm[, aggr_mode := ifelse(sector %in% c("trn_aviation_intl"), "Pass. non LDV", aggr_mode)]
+  demandpkm[, aggr_mode := ifelse(subsector_L2 %in% c("trn_pass_road_bus", "HSR_tmp_subsector_L2", "Passenger Rail_tmp_subsector_L2", "Cycle_tmp_subsector_L2", "Walk_tmp_subsector_L2", "Domestic Aviation_tmp_subsector_L2", "Bus") | subsector_L1 %in% c("trn_pass_road_LDV_2W"), "Pass. non LDV", aggr_mode)]
+
+  demandpkm[, veh := ifelse(vehicle_type %in% c("Truck (0-1t)", "Truck (0-3.5t)"), "Trucks (<3.5t)", "Trucks (<3.5t)")]
+  demandpkm[, veh := ifelse(vehicle_type %in% c("Truck (16-32t)", "Truck (3.5-16t)", "Truck (6-15t)"), "Trucks (3.5t-16)", veh)]
+  demandpkm[, veh := ifelse(vehicle_type %in% c("Truck (>15t)", "Truck (16-32t)", "Truck (>32t)" ), "Trucks (>16)", veh)]
+  demandpkm[, veh := ifelse(grepl("Large|SUV|Midsize|Multipurpose Vehicle|Van|3W Rural", vehicle_type), "Large Cars", veh)]
+  demandpkm[, veh := ifelse(grepl("Subcompact|Compact|Mini|Three-Wheeler_tmp_vehicletype", vehicle_type), "Small Cars", veh)]
+  demandpkm[, veh := ifelse(grepl("Motorcycle|Moped|Scooter", vehicle_type), "Motorbikes", veh)]
+  demandpkm[, veh := ifelse(grepl("bus|Bus", vehicle_type), "Bus", veh)]
+  demandpkm[, veh := ifelse(subsector_L3 == "Domestic Aviation", "Domestic Aviation", veh)]
+  demandpkm[, veh := ifelse(subsector_L3 == "International Aviation", "International Aviation", veh)]
+  demandpkm[, veh := ifelse(grepl("Freight Rail", vehicle_type), "Freight Rail", veh)]
+  demandpkm[, veh := ifelse(grepl("Passenger Rail|HSR", vehicle_type), "Passenger Rail", veh)]
+  demandpkm[, veh := ifelse(grepl("Ship", vehicle_type), "Shipping", veh)]
+  demandpkm[, veh := ifelse(grepl("Cycle|Walk", subsector_L3), "Non motorized", veh)]
+  demandpkm = demandpkm[,.(demand_F = sum(demand_F)), by = c("iso", "year", "aggr_mode", "veh")]
+  setnames(demandpkm, old = "veh", new = "vehicle_type")
+
+
+  demandpkm[, vehicle_type_plot := factor(vehicle_type, levels = c("LDV","Freight Rail", "Trucks (<3.5t)", "Trucks (3.5t-16)",  "Trucks (>16)", "Trucks",
+                                                                "Motorbikes", "Small Cars", "Large Cars", "Van",
+                                                                "Domestic Aviation", "International Aviation","Bus", "Passenger Rail",
+                                                                "Freight", "Non motorized", "Shipping"))]
+
+
+  demandpkm[, mode := ifelse(vehicle_type %in% c("Freight", "Freight Rail", "Trucks", "Trucks (3.5t-16)",  "Trucks (>16)", "Shipping"),"freight", "pass")]
+
+  demandpkm = merge(demandpkm, REMIND2ISO_MAPPING, by = "iso")
+  demandpkm = demandpkm[, .(demand_F = sum(demand_F)), by = c("region", "year", "vehicle_type_plot", "aggr_mode", "mode")]
+
+    demandpkm = demandpkm[order(aggr_mode)]
+
+  p = ggplot()+
+    geom_area(data = demandpkm[mode =="pass"& year > 2010], aes(x=year, y=demand_F, group = interaction(vehicle_type_plot,aggr_mode), fill = vehicle_type_plot), color="black",position= position_stack())+
+    labs(x = "", y = "Energy Services demand [trillion pkm]")+
+    facet_wrap(~region, nrow = 4)
+    theme_minimal()+
+    # scale_fill_manual("Vehicle Type",values = cols, breaks=legend_ord)+
+    theme(axis.text.x = element_text(size = 7),
+          axis.text.y = element_text(size=7),
+          axis.title = element_text(size = 8),
+          title = element_text(size = 8),
+          legend.text = element_text(size = 8),
+          legend.title = element_text(size = 8),
+          strip.text = element_text(size=8))
+
+
+  return(p)
 
-```{r, echo=FALSE}
-FE_modes_bar1=function(demandEJ, msect="trn_pass", region_plot){
-  if(msect == "trn_pass")
-    sector_display = "Passenger"
-  if(msect == "trn_freight")
-    sector_display = "Freight"
-
-  ## group by subsector_L1 and summarise the demand
-  df=demandEJ[, .(demand_EJ=sum(demand_EJ)),
-              by = c("EDGE_scenario", "region", "year","subsector_L1","sector")]
-  df=df[order(year) & year %in% years_plot]
-  ## give proper names to the categories
-  df=merge(df, edgeTrpLib::L1mapping, all.x=TRUE, by="subsector_L1")
-  #plot
-  plot_p=ggplot()+
-    geom_bar(data=df%>%filter(sector == msect, region == region_plot),
-             aes(x=year,y=demand_EJ,group=mode,fill=mode),
-             position=position_stack(),stat="identity",color="black")+
-    facet_wrap(~EDGE_scenario) +
-    ylab("Energy (EJ)") +
-    ggtitle(paste0("Final Energy Demand by Mode, ", sector_display, ", ", region_plot))+
-    theme(axis.text.x = element_text(angle = 90))+
-    scale_x_continuous(breaks=years_plot)
-  
-  return(plot_p)
 }
 
-FE_modes_bar1(demandEJ = demandF_plot_EJ, region_plot = region_plot)
+## energy services demand
+demandpkmplotf(demand_km)
 ```
 
+# CO2 intensity of new sales
 
+```{r, echo=FALSE, warning=FALSE}
 
-```{r, echo=FALSE}
-FE_modes_bar_oilVSelec=function(demandEJ, region_plot){
-  sector_display = "Total transport"
-  ## group by subsector_L1 and summarise the demand
-  df=demandEJ[, .(demand_EJ=sum(demand_EJ)),
-              by = c("EDGE_scenario", "region", "year","subsector_L1","sector","technology")]
-  df[, tech_plot := ifelse(technology %in% c("BEV","Electric"), "Electriciy", NA)]
-  df[, tech_plot := ifelse(technology %in% c("Liquids", "Hybrid Liquids"), "Liquids", tech_plot)]
-  df=df[!is.na(tech_plot),] ## only liquids and electric driven entries interesting
-  df=df[order(year) & year %in% years_plot]
-  ## give proper names to the categories
-  df=merge(df, edgeTrpLib::L1mapping, all.x=TRUE, by="subsector_L1")
-  df[,short_names:=ifelse(mode %in% c("Buses","Rail Passenger","High Speed Rail"),"Other Passenger",NA)]
-  df[,short_names:=ifelse(mode %in% c("2W","4W","Three Wheelers"),"LDV",short_names)]
-  df[,short_names:=ifelse(mode %in% c("Domestic Aviation","International Aviation"),"Aviation",short_names)]
-  df[,short_names:=ifelse(mode %in% c("International Shipping","Domestic Shipping"),"Shipping",short_names)]
-  df[,short_names:=ifelse(mode %in% c("Road Freight","Rail Freight"),"Road and Rail Freight",short_names)]
-  
-  #plot
-  plot_p=ggplot()+
-    geom_bar(data=df%>%filter(region == region_plot, year ==2050),
-             aes(x=tech_plot,y=demand_EJ,group=short_names,fill=short_names),
-             position=position_stack(),stat="identity",alpha=0.95)+
-    facet_wrap(~EDGE_scenario) +
-    ylab("Energy (EJ)") +
-    ggtitle(paste0("Final Energy Demand by Mode in 2050, total transport ", region_plot))+
-    theme_light()+
-    theme(axis.text.x = element_text(angle = 90),
-          axis.title.x = element_blank(),
-          strip.text.x = element_text(size = 13, color = "black"),
-          strip.background=element_rect(fill="white"),
-          axis.text = element_text(size=13),
-          title = element_text(size=13),
-          legend.text = element_text(size=13))+
-    scale_fill_brewer(palette = "Set2")+
-    guides(fill=guide_legend(title="Transport mode"))
-  
-  return(plot_p)
-}
+CO2km_intensity_newsalesplotf = function(annual_sales, mj_km_data, sharesVS1, shares_source_liquids){
+  shares_source_liquids[, technology := ifelse(variable %in% c("FE|Transport|Liquids|Oil", "FE|Transport|Liquids|Coal"), "Oil", "Biodiesel")]
+  shares_source_liquids = shares_source_liquids[,.(value = sum(value)), by = c("model","scenario","region", "period", "unit","technology")]
+  shares_source_liquids = shares_source_liquids[region != "World"]
+  shares_source_liquids[, region:=as.character(region)]
+  shares_source_liquids[, year := period]
+  shares_source_liquids[, period:=NULL]
+  gdp <- getRMNDGDP(scenario = "SSP2", usecache = T)
+  shares_source_liquids <- disaggregate_dt(shares_source_liquids, REMIND2ISO_MAPPING,
+                                           valuecol="value",
+                                           datacols=c("model","scenario", "unit","technology"),
+                                           weights=gdp)
 
-FE_modes_bar_oilVSelec(demandEJ = demandF_plot_EJ, region_plot = region_plot)
-```
+  shares_source_liquids[, shareliq := value/sum(value),by=c("iso", "year")]
 
+  # ## CO2 content
+  # CO2_petrol = 3.1 ## gCO2/gFUEL
+  # CO2_biodiesel = 2.7 ## TODO this number is made up!
+  # CO2_cng = 2.7 ## gCO2/gFUEL
 
-## FE composition
-
-```{r, echo=FALSE}
-
-FE_modeshares_area=function(demandEJ){
-  #group by subsector_L3 and summarise the demand
-  df=demandEJ[, .(demand_EJ=sum(demand_EJ)),
-              by = c("EDGE_scenario", "region", "year","subsector_L3")]
-  #order by year
-  df=df[order(year)]
-  df=df[year>=2005,]
-  #plot
-  plot=ggplot()+
-    geom_area(data=df %>% filter(year <= max(years_plot), region == region_plot),
-              aes(x=year,y=demand_EJ,group=subsector_L3,fill=subsector_L3),
-              color="black")+
-    facet_wrap(~EDGE_scenario)+
-    ylab("Energy (EJ)") +
-    ggtitle("Final Energy Demand, Mode composition")+
-    theme(axis.text.x = element_text(angle = 90, hjust = 1))
-  
-  return(plot)
-}
+  ## TODO of CO2 content of biodiesel is made up! gCO2/gFUEL
+  emi_fuel = data.table(technology = c("Oil", "Biodiesel", "NG"), ei_gF_MJ = c(20, 20, 20), emi_cGO2_gF = c(3.1, 3.1, 2.7))
 
-p=FE_modeshares_area(demandEJ = demandF_plot_EJ)
-p
-## ggsave("FE_modeshares.png")
-```
+  emi_liquids = merge(shares_source_liquids, emi_fuel, all.x = TRUE, by = "technology")
+  emi_liquids = emi_liquids[, .(ei_gF_MJ = sum(shareliq*ei_gF_MJ), emi_cGO2_gF = sum(shareliq*emi_cGO2_gF)), by = c("iso", "year")][, technology := "Liquids"]
+  emi_NG = cbind(emi_fuel[technology == "NG"], unique(shares_source_liquids[,c("year", "iso")]))
 
+  emi_fuel = rbind(emi_NG, emi_liquids)
+  emi_fuel[, gCO2_MJ := ei_gF_MJ*emi_cGO2_gF]
 
+  emi_fuel = merge(mj_km_data[subsector_L1 == "trn_pass_road_LDV_4W"], emi_fuel, all.x = TRUE, by = c("iso", "year", "technology"))
+  emi_fuel[is.na(gCO2_MJ) & !technology %in% c("Liquids", "NG"), gCO2_MJ := 0]
+  emi_fuel[, gCO2_km := MJ_km * gCO2_MJ]
 
-```{r, echo=FALSE}
-## fuel use by sector
-fuel_shares_area=function(demandEJ, msect="trn_pass", region_plot){
-  if(msect == "trn_pass")
-    sector_display = "Passenger"
-  if(msect == "trn_freight")
-    sector_display = "Freight"
-  ##group by sector and technology and summarise demand
-  df=demandEJ[, .(demand_EJ=sum(demand_EJ)),
-              by = c("EDGE_scenario", "region", "year","technology","sector")]
-
-  df=df[order(year) & year>=2005,]
-  #plot
-  plot1=ggplot()+
-    geom_area(data=df%>%filter(sector == msect, year <= max(years_plot), region == region_plot),
-              aes(x=year,y=demand_EJ,group=technology,fill=technology),position="fill")+
-    facet_wrap(~EDGE_scenario)+
-    ggtitle(paste0("Final Energy Demand, ", sector_display, ", Fuel Composition"))+
-    ylab("Share")
-    theme(axis.text.x = element_text(angle = 90, hjust = 1))
-    
-  return(plot1)
-}
+  totalemi = merge(emi_fuel, annual_sales, all.y = TRUE, by = c("iso", "year", "technology", "vehicle_type", "subsector_L1"), all.x = TRUE)
+  totalemi = totalemi[!is.na(share) & !is.na(gCO2_km)]
+  totalemi[, gCO2_km_ave := gCO2_km*share]
 
-p=fuel_shares_area(demandEJ = demandF_plot_EJ, region_plot = region_plot)
-p
 
-```
+  ##totalemi = merge(totalemi, demand_ej_plot)
 
 
-```{r, echo=FALSE}
+  totalemi = totalemi[,.(gCO2_km_ave = sum(gCO2_km_ave)), by = c("year", "iso", "vehicle_type")]
 
-SW_trend_plot = function(FV_SW,sector_plot){
-  if (sector_plot == "trn_pass") {
-  FV_SW=FV_SW[iso=="DEU" & vehicle_type =="Large Car and SUV",]
-  } else if (sector_plot =="trn_freight"){
-  FV_SW=FV_SW[iso=="DEU" & vehicle_type =="Truck (16-32t)",]
-  } else if (sector_plot =="trn_aviation_intl"){
-  FV_SW=FV_SW[iso=="DEU" & subsector_L3 =="International Aviation",]
-  } else if (sector_plot =="trn_shipping_intl"){
-  FV_SW=FV_SW[iso=="DEU" & subsector_L3 =="International Ship",]
-  }
-  FV_SW[,type:=ifelse(technology=="Liquids", "Conventional ICE (Liquid fuels)",NA)]
-  FV_SW[,type:=ifelse(technology=="NG", "Natural Gas ICE",type)]
-  FV_SW[,type:=ifelse(technology=="BEV", "Alternative fuels: BEV",type)]
-  FV_SW[,type:=ifelse(technology=="FCEV", "Alternative fuels: FCEV",type)]
-  FV_SW[,type:=ifelse(technology=="Hybrid Liquids", "Unconventional ICE (Hybrid)",type)]
-  
-  p=ggplot()+
-    geom_line(data=FV_SW%>%filter(year>= min(years_plot), year<=max(years_plot)),aes(x=year,y=sw,group=type,color=type),alpha = 0.8,size=1.5)+
-    theme_light()+
-    facet_wrap(~EDGE_scenario)+
-    theme(axis.text.x = element_text(angle = 90, hjust = 1),
-          axis.text = element_text(size=13),
-          title = element_text(size=14),
-          legend.text = element_text(size=13))+
-    scale_x_continuous(breaks=years_plot)+
-    xlab("Year")+
-    ylab ("Preference factors tech. types [-]")+
-    ggtitle(paste0("Preference factors trend for tech. types for ", sector_plot, " [-]"))+
-    theme(strip.text.x = element_text(size=13,color="black"),
-          strip.background = element_rect(fill="white",color = "black"))+
-    scale_color_discrete(name="Technology type")
-  return(p)
-  }
+  totalemi = merge(totalemi, sharesVS1, all.x = TRUE, by = c("iso", "year", "vehicle_type"))
 
-p=SW_trend_plot(FV_SW=sw_tech,sector_plot)
-p
-```
+  totalemi = totalemi[,.(gCO2_km_ave = sum(gCO2_km_ave*share)), by = c("iso", "year", "subsector_L1")]
+
+
+  totalemi = merge(totalemi, REMIND2ISO_MAPPING, by="iso")
+  totalemi = merge(totalemi, gdp, all.x=TRUE, by = c("iso", "year"))
+  totalemi[, share := weight/sum(weight), by = c("year", "region")]
+  totalemi = totalemi[,.(gCO2_km_ave = sum(gCO2_km_ave*share)), by = c("year", "region")]
 
+  p = ggplot()+
+    geom_line(data = totalemi, aes(x = year, y = gCO2_km_ave))+
+    labs(title = "gCO2/km average", y = "Average gCO2/km LDVs new additions")+
+    facet_wrap(~region, nrow = 4)+
+    theme_minimal()
+
+  return(p)
+}
 
+shares_source_liquids = miffile[variable %in% c("FE|Transport|Liquids|Biomass", "FE|Transport|Liquids|Coal", "FE|Transport|Liquids|Oil"),]
+CO2km_intensity_newsalesplotf(annual_sales, mj_km_data, sharesVS1 = shares$VS1_shares, shares_source_liquids)
+```
\ No newline at end of file