Boa tarde colegas programadores! Estou tendo dificuldades em acertar o formato que devo organizar os outputs de um looping que fiz. Tenho um banco de dados com datas e algumas variáveis climáticas. Dessas datas, preciso simular alguns cálculos para dias específicos. São 36 simulações por ano, durante um período de 34 anos, totalizando 1224 simulações. Para alcançar esse objetivo estou trabalhando com 3 comandos for, um dentro do outro. O primeiro seleciona as datas/posições que devem ser feitos os cálculos. O segundo executa os devidos cálculos para as datas selecionadas e o ultimo integra os valores calculados, gerando a informação de interesse. Acredito que a logica de cálculos dentro de cada for() está ok, mas a saída não está puxando os valores que o looping calcula. Segue o script no corpo do email e, em anexo, um exemplo de um ano do banco de dados que está sendo utilizado. Desde já, agradeço a colaboração dos colegas. Abs #--------------------Remover objetos salvos na memória do R--------------------# rm(list = ls()) #--------------------Seleção de diretorios e arquivos de interesse--------------------# a = read.table("x.txt", sep = "\t", header = T) #--------------------Calculo das Variaves Astronomicas--------------------# head(a) a$DATE = as.Date(levels(a$DATE))[a$DATE] ano = format(a$DATE, trim = T, '%Y'); a$ano=ano mes = format(a$DATE, trim = T, '%m') ; a$mes=mes dia = format(a$DATE, trim = T, '%d') ; a$dia=dia dece = ifelse(dia<11, 1,ifelse(dia<21,2,3)); a$decendio=dece lat = -7.53 corr = pi/180 decl = 23.45*sin(corr*((a$DiaJuliano-80)*360/365)); decl hn = 1/corr*acos(-tan(corr*lat)*tan(corr*decl)) a$N = 2*hn/15 s = (1+0.033*cos(a$DiaJuliano*360/365)) t = corr*hn*sin(corr*lat)*sin(corr*decl) u = cos(corr*lat)*cos(corr*decl)*sin(corr*hn) a$Qo = 37.6*s*(t+u) head(a) x = 0.25 b = 0.50 n.est = (a$SRAD/a$Qo-x)*a$N/b; n.est a$n<- ifelse((n.est)<0,1, n.est); head(a) a$n.N = a$n/a$N; head(a) #--------------------Transformar MJ em caloria--------------------# a$Srad.cal = round((a$SRAD*1000000)/41800,1); head(a) #round é para delimitar o numero de casas apos a virgula a$sem = paste0(dia, '-', mes) a$TMED = (((a$TMAX+a$TMIN))/2) head(a) #-------------------- Determinacao das datas de simulacao--------------------# ndc = 140 datas = 1:366 # referentes aos dia juliano dim(a) anos = table(a$ano); anos potencial = matrix(NA, nrow = 366, ncol = 34 ) rownames(potencial) = 1:366 colnames(potencial) = c(1980:2013) head(potencial); dim(potencial) head(a) for (i in (a$DiaJuliano)){ epocas = a[which( a$sem == '01-01'|a$sem == '11-01' |a$sem =='21-01'| a$sem == '01-02'|a$sem == '11-02' |a$sem =='21-02'| a$sem == '01-03'|a$sem == '11-03' |a$sem =='21-03'| a$sem == '01-04'|a$sem == '11-04' |a$sem =='21-04'| a$sem == '01-05'|a$sem == '11-05' |a$sem =='21-05'| a$sem == '01-06'|a$sem == '11-06' |a$sem =='21-06'| a$sem == '01-07'|a$sem == '11-07' |a$sem =='21-07'| a$sem == '01-08'|a$sem == '11-08' |a$sem =='21-08'| a$sem == '01-09'|a$sem == '11-09' |a$sem =='21-09'| a$sem == '01-10'|a$sem == '11-10' |a$sem =='21-10'| a$sem == '01-11'|a$sem == '11-11' |a$sem =='21-11'| a$sem == '01-12'|a$sem == '11-12' |a$sem =='21-12'& a$ano>=1980),] semear = length(epocas[,1]) for (j in length(semear)){ ini = as.numeric(rownames(epocas[j])) fim = (ini+ndc)-1 pac = a[(ini:fim),] #--------------------Inserindo as Produtividades Parciais--------------------# #----------Inserindo influencia da nebulosidade----------# pac$ctc[j] = -4.16+(0.4325*pac$TMED[j])-(0.00725*((a$TMED[j])^2)) pac$ctn[j] = -1.064+(0.173*a$TMED[j])-(0.0029*((a$TMED[j])^2)) pac$PPBc[j] = (107.2+(0.36*pac$Qo[j]*1000000)/41800)*pac$ctc[j]*pac$n.N[j] pac$PPBn[j] = (31.7+(0.219*pac$Qo[j]*1000000)/41800)*pac$ctn[j]*(1-(pac$n.N[j])) pac$PPB[j] = round((pac$PPBc[j])+(pac$PPBn[j]),1); head(pac) #--------------------Calculo da Produtividade - base seca --------------------# pac$c.resp[j] = ifelse(a$TMED[j] < 20, 0.6, 0.5) pac$c.colheita[j] = 0.5 umidade = 13 c.um = (1/(1-0.01*umidade)) pac$iaf[j] = 5 pac$c.iaf[j] = (0.0093+0.185*pac$iaf[j])-(0.0175*pac$iaf[j]^2) pac$PP.Liquida_seca = NA pac$PP.Liquida_seca[1] = round(pac$PPB[1]*pac$c.resp[1]*pac$c.colheita[1]*pac$c.iaf[1],1); head(pac) # em kg/ha for (w in 2:ndc){ pac$PP.Liquida_seca[w] = round(pac$PPB[w]*pac$c.resp[w]*pac$c.colheita[w]*pac$c.iaf[w]+ pac$PP.Liquida_seca[w-1],1) } #--------------------Soma da PP--------------------# pac$PP_final = round(pac$PP.Liquida_seca[w]*c.um,1) dim(pac) PP = pac[ndc,32] potencial[i,j] = PP } print(i) } View(pac) View(potencial) Yury Duarte Engenheiro Agrônomo - ESALQ/USP