# Campo magnético em corte sagital library(ggplot2) dados <- data.frame( d_cm = c(0, 2, 4, 6, 13, 21, 26, 30, 39, 43, 51, 54, 56, 58, 60), B_uT = c(26.3, 45.7, 84.5, 110.5, 50.9, 29.7, 28.6, 28.8, 40.9, 55.3, 112.9, 41.9, 25.2, 17.5, 3.8) ) d1 <- 6 d2 <- 51 B_wire <- function(x, y, y0, I = 1, k = 1, eps = 0.8) { dx <- x dy <- y - y0 r2 <- dx^2 + dy^2 + eps^2 Bx <- -k * I * dy / r2 By <- k * I * dx / r2 list(Bx = Bx, By = By) } B_pair <- function(x, y, y1, y2, I = 1, k = 1, eps = 0.8) { b1 <- B_wire(x, y, y0 = y1, I = I, k = k, eps = eps) b2 <- B_wire(x, y, y0 = y2, I = -I, k = k, eps = eps) list(Bx = b1$Bx + b2$Bx, By = b1$By + b2$By) } x_seq <- seq(-20, 20, by = 2) y_seq <- seq(0, 60, by = 2) grid <- expand.grid(x = x_seq, y = y_seq) b <- B_pair(grid$x, grid$y, y1 = d1, y2 = d2, I = 1, k = 1, eps = 1.0) grid$Bx <- b$Bx grid$By <- b$By grid$Bmag <- sqrt(grid$Bx^2 + grid$By^2) x_meas <- 0 linha <- data.frame(x = x_meas, y = dados$d_cm) bl <- B_pair(linha$x, linha$y, y1 = d1, y2 = d2, I = 1, k = 1, eps = 1.0) linha$Bmag_model <- sqrt(bl$Bx^2 + bl$By^2) scale_fac <- sum(linha$Bmag_model * dados$B_uT) / sum(linha$Bmag_model^2) grid$Bx_uT <- scale_fac * grid$Bx grid$By_uT <- scale_fac * grid$By grid$Bmag_uT <- scale_fac * grid$Bmag arrow_fac <- 0.20 grid$xend <- grid$x + arrow_fac * grid$Bx_uT grid$yend <- grid$y + arrow_fac * grid$By_uT fios <- data.frame(x = 0, y = c(d1, d2), fio = c("fio 1", "fio 2")) pts <- data.frame(x = 0, y = dados$d_cm, B_uT = dados$B_uT) ggplot() + geom_segment( data = grid, aes(x = x, y = y, xend = xend, yend = yend), arrow = arrow(length = unit(1.5, "mm")), linewidth = 0.25, alpha = 0.7 ) + geom_point(data = fios, aes(x = x, y = y), size = 3) + geom_line(data = data.frame(x = x_meas, y = y_seq), aes(x = x, y = y), linewidth = 0.5, linetype = "dashed") + geom_point(data = pts, aes(x = x, y = y, size = B_uT), alpha = 0.7) + scale_size_continuous(name = "B (µT)", range = c(1, 6)) + coord_equal() + labs( x = "x (cm) — lateral (só para visualização)", y = "d (cm) — posição ao longo da régua", title = "Representação vectorial (modelo) do campo de dois condutores antiparalelos\n(escalado para µT com os pontos do gráfico)" ) + theme_minimal() # Campo magnético linearizado library(ggplot2) dados <- data.frame( d_cm = c(0, 2, 4, 6, 13, 21, 26, 30, 39, 43, 51, 54, 56, 58, 60), B_uT = c(26.3, 45.7, 84.5, 110.5, 50.9, 29.7, 28.6, 28.8, 40.9, 55.3, 112.9, 41.9, 25.2, 17.5, 3.8) ) y1 <- 6 y2 <- 51 Bz_pair <- function(y, y1, y2, K = 1, eps = 0.8) { r1 <- sqrt((y - y1)^2 + eps^2) r2 <- sqrt((y - y2)^2 + eps^2) K * ( (y - y1)/r1^2 - (y - y2)/r2^2 ) } fitK <- function(K) { pred <- abs(Bz_pair(dados$d_cm, y1, y2, K = K, eps = 1.0)) sum((pred - dados$B_uT)^2) } Kopt <- optimize(fitK, interval = c(0, 1e6))$minimum yy <- seq(0, 60, by = 1) Bz <- Bz_pair(yy, y1, y2, K = Kopt, eps = 1.0) arrow_fac <- 0.01 dfv <- data.frame(x = 0, y = yy, z = Bz) dfv$xend <- dfv$x dfv$yend <- dfv$y + arrow_fac * dfv$z fios <- data.frame(x = 0, y = c(y1, y2)) ggplot() + geom_segment( data = dfv, aes(x = x, y = y, xend = xend, yend = yend), arrow = arrow(length = unit(1.5, "mm")), linewidth = 0.35, alpha = 0.7 ) + geom_point(data = fios, aes(x = x, y = y), size = 3) + geom_point(data = dados, aes(x = 0.5, y = d_cm, size = B_uT), alpha = 0.7) + scale_size_continuous(name = "B (µT)", range = c(1, 6)) + coord_equal(xlim = c(-2, 4)) + labs( x = "(só para separar camadas)", y = "d (cm) ao longo da régua", title = "Campo vectorial correcto: B aponta ±z (fora/dentro do papel)\nSetas verticais = sinal de Bz; tamanho dos pontos = |B| medido" ) + theme_minimal() rad_juntos<-read.table("~/Downloads/RADPi0.dat", header=TRUE, sep="\t", dec=".") X <- rad_juntos[c("potmin1","potmin2","potmin3")] rad_juntos$Bmean <- rowMeans(X, na.rm=TRUE) dat <- subset(rad_juntos, Bmean > 0, select=c(dist, Bmean)) dat$x <- dat$dist / 100 x0_start <- dat$x[which.max(dat$Bmean)] B0_start <- min(dat$Bmean) H_start <- max(dat$Bmean) - B0_start w_start <- 0.01 p_start <- 0.8 fit <- nls( Bmean ~ B0 + H / (1 + ((x - x0)/w)^2)^p, data = dat, start = list(B0=B0_start, H=H_start, x0=x0_start, w=w_start, p=p_start), algorithm = "port", lower = c(B0=-Inf, H=0, x0=min(dat$x), w=1e-6, p=0.1), upper = c(B0= Inf, H= Inf, x0=max(dat$x), w= Inf, p=5) ) coef(fit) plot(dat$dist, dat$Bmean, pch=21, bg="red", xlab="distância (cm)", ylab = expression("|B|"~"(" * n * "T)")) ord <- order(dat$dist) lines(dat$dist[ord], predict(fit)[ord]) box(lwd=3) win <- subset(rad_juntos, dist >= 25 & dist <= 35) Bmin <- rowMeans(win[c("potmin1","potmin2","potmin3")], na.rm=TRUE) Bmed <- rowMeans(win[c("potmed1","potmed2","potmed3")], na.rm=TRUE) Bmax <- rowMeans(win[c("potmax1","potmax2","potmax3")], na.rm=TRUE) auc_trapz <- function(x, y) { o <- order(x); x <- x[o]; y <- y[o] sum(diff(x) * (head(y, -1) + tail(y, -1)) / 2) } metrics <- function(dist, y) { i <- which.max(y) data.frame( Bpeak_nT = y[i], dist_pico_cm = dist[i], AUC_nT_cm = auc_trapz(dist, y), media_25_35_nT = mean(y), row.names = NULL ) } resumo <- rbind( cbind(potencia="min (calor=1)", metrics(win$dist, Bmin)), cbind(potencia="med (calor=2)", metrics(win$dist, Bmed)), cbind(potencia="max (calor=3)", metrics(win$dist, Bmax)) ) num_cols <- c("Bpeak_nT","dist_pico_cm","AUC_nT_cm","media_25_35_nT") resumo[num_cols] <- lapply(resumo[num_cols], as.numeric) resumo$rel_Bpeak_vs_min <- resumo$Bpeak_nT / resumo$Bpeak_nT[resumo$potencia=="min (calor=1)"] resumo$rel_AUC_vs_min <- resumo$AUC_nT_cm / resumo$AUC_nT_cm[resumo$potencia=="min (calor=1)"] resumo mu0 <- 4*pi*1e-7 d <- 0.003 win <- subset(rad_juntos, dist >= 25 & dist <= 35) x_m <- win$dist/100 Bmin <- rowMeans(win[c("potmin1","potmin2","potmin3")], na.rm=TRUE) * 1e-9 Bmed <- rowMeans(win[c("potmed1","potmed2","potmed3")], na.rm=TRUE) * 1e-9 Bmax <- rowMeans(win[c("potmax1","potmax2","potmax3")], na.rm=TRUE) * 1e-9 dat <- data.frame( x = rep(x_m, 3), B = c(Bmin, Bmed, Bmax), cond = factor(rep(c("min","med","max"), each=length(x_m))) ) x0_start <- x_m[which.max(Bmax)] B0_start <- min(dat$B) b2_start <- (0.005)^2 Bpk_min <- max(Bmin); Bpk_med <- max(Bmed); Bpk_max <- max(Bmax) Kmin_start <- max((Bpk_min - B0_start) * b2_start, 0) Kmed_start <- max((Bpk_med - B0_start) * b2_start, 0) Kmax_start <- max((Bpk_max - B0_start) * b2_start, 0) fit_phys <- nls( B ~ B0 + (Kmin*(cond=="min") + Kmed*(cond=="med") + Kmax*(cond=="max")) / ((x - x0)^2 + b2), data = dat, start = list(B0=B0_start, x0=x0_start, b2=b2_start, Kmin=Kmin_start, Kmed=Kmed_start, Kmax=Kmax_start), algorithm="port", lower = c(B0=-Inf, x0=min(x_m), b2=1e-10, Kmin=0, Kmed=0, Kmax=0) ) p <- coef(fit_phys) xgrid <- seq(min(x_m), max(x_m), length.out=400) Bfit_min <- (p["B0"] + p["Kmin"]/((xgrid - p["x0"])^2 + p["b2"])) * 1e9 Bfit_med <- (p["B0"] + p["Kmed"]/((xgrid - p["x0"])^2 + p["b2"])) * 1e9 Bfit_max <- (p["B0"] + p["Kmax"]/((xgrid - p["x0"])^2 + p["b2"])) * 1e9 Bmin_nT <- Bmin * 1e9 Bmed_nT <- Bmed * 1e9 Bmax_nT <- Bmax * 1e9 o <- order(win$dist) ylim <- range(c(Bmin_nT, Bmed_nT, Bmax_nT, Bfit_min, Bfit_med, Bfit_max), na.rm=TRUE) plot(win$dist[o], Bmin_nT[o], type="n", xlab="distância (cm)", ylab="|B| (nT)", ylim=ylim, main="25–35 cm: dados (pontos) e ajuste de modelo (linhas)") #lines(xgrid*100, Bfit_min, lty=1) lines(xgrid*100, Bfit_med, lty=1, lwd=2, col="grey") lines(xgrid*100, Bfit_max, lty=1, lwd=2, col="grey") #points(win$dist[o], Bmin_nT[o], pch=16) points(win$dist[o], Bmed_nT[o], pch=21, bg="pink") points(win$dist[o], Bmax_nT[o], pch=21, bg="red") legend("topright", legend=c("méd dados","máx dados", "méd ajuste","máx ajuste"), pch=c(21,21,NA,NA), pt.bg=c("pink","red",NA,NA), col=c("black","black","grey", "grey"), lty=c(NA,NA,1,1), lwd=c(NA,NA,2,2), bty="n") box(lwd=3) # Claude: Modelo teórico de |B| (RMS) para dois fios paralelos (fase e neutro) # Transecto ortogonal: 15cm antes do 1.º cabo + 15cm entre cabos + 15cm depois do 2.º cabo # # Convenção geométrica: # - eixo x ao longo do transecto # - 1.º cabo (fase) em x = 0 # - 2.º cabo (neutro) em x = d (d = 0.15 m) # - sensor a uma altura efectiva h acima do centro do condutor (m) # # Saída: |B|_RMS em microtesla (uT), por omissão. B_transect_uT <- function(x, d = 0.15, # separação entre cabos (m) P = 3000, # potência (W) V = 230, # tensão (V RMS) pf = 1.0, # factor de potência (1 = carga resistiva) h = 0.005 + 0.0015, # altura efectiva (m): 5mm + isolamento (exemplo) B0_uT = 0, # fundo (uT) return_T = FALSE) { # se TRUE, devolve em tesla (T) # corrente RMS I_rms <- P / (V * pf) # distâncias ao eixo de cada fio r1 <- sqrt((x - 0)^2 + h^2) r2 <- sqrt((x - d)^2 + h^2) # |B| RMS (tesla) para correntes opostas: | mu0 I /(2pi) * (1/r1 - 1/r2) | mu0_over_2pi <- 2e-7 # mu0/(2*pi) = 2e-7 (SI) B_T <- mu0_over_2pi * I_rms * abs(1/r1 - 1/r2) # adiciona fundo B_uT <- (B_T * 1e6) + B0_uT if (return_T) return(B_uT / 1e6) return(B_uT) } campo_magnetico_transecto <- function() { mu0 <- 4 * pi * 1e-7 V <- 230 P <- 3000 I <- P / V # ≈ 13.04 A RMS x1 <- 0.15 # fio 1 (fase) [m] x2 <- 0.30 # fio 2 (neutro) [m] h <- 0.00669 # distância vertical sensor–centro do fio [m] x <- seq(0, 0.45, length.out = 2000) r1 <- sqrt((x - x1)^2 + h^2) r2 <- sqrt((x - x2)^2 + h^2) # Prefactor de cada fio B1 <- (mu0 * I) / (2 * pi * r1) # módulo do campo do fio 1 B2 <- (mu0 * I) / (2 * pi * r2) # módulo do campo do fio 2 # Versor perpendicular ao raio, para cada fio (componentes x e z) # Fio 1: corrente +I → regra da mão direita B1x <- B1 * ( h / r1) B1z <- B1 * (-(x - x1) / r1) # Fio 2: corrente -I (sentido oposto ao fio 1) B2x <- -B2 * ( h / r2) B2z <- -B2 * (-(x - x2) / r2) # Superposição vectorial Bx_tot <- B1x + B2x Bz_tot <- B1z + B2z # Módulo total B_uT <- sqrt(Bx_tot^2 + Bz_tot^2) * 1e6 # em µT # --- Plot --- B_max <- max(B_uT) x_max <- x[which.max(B_uT)] ICNIRP_pub <- 200 par(mar = c(5, 5, 4, 2), bg = "grey97") plot(x * 100, B_uT, type = "l", lwd = 2.5, col = "#1a6faf", xlab = "Posição no transecto (cm)", ylab = expression("|B| (µT)"), main = expression("|B| ao longo do transecto — fase + neutro, 230 V, 3000 W"), xlim = c(0, 45), ylim = c(0, B_max * 1.25), las = 1, bty = "l") abline(h = ICNIRP_pub, lty = 2, col = "firebrick", lwd = 1.5) text(44, ICNIRP_pub * 1.07, paste0("ICNIRP ", ICNIRP_pub, " µT"), col = "firebrick", adj = 1, cex = 0.8) abline(v = c(x1, x2) * 100, lty = 3, col = c("darkorange", "darkgreen"), lwd = 1.5) text(x1 * 100 + 0.4, B_max * 1.18, "Fio 1\n(fase)", col = "darkorange", adj = 0, cex = 0.78) text(x2 * 100 + 0.4, B_max * 1.18, "Fio 2\n(neutro)", col = "darkgreen", adj = 0, cex = 0.78) points(x_max * 100, B_max, pch = 19, col = "#1a6faf", cex = 1.3) text(x_max * 100, B_max * 1.08, sprintf("máx = %.1f µT\n@ %.1f cm", B_max, x_max * 100), col = "#1a6faf", cex = 0.78, adj = 0.5) legend("topright", legend = c( sprintf("I ≈ %.2f A (RMS)", I), sprintf("h ≈ %.1f mm", h * 1000), "Distância entre fios = 15 cm"), bty = "n", cex = 0.78, text.col = "grey30") invisible(data.frame(x_cm = x * 100, B_uT = B_uT)) } dados <- campo_magnetico_transecto() # ChatGPT: Campo magnetico teorico |B| (RMS) ao longo de um transecto ortogonal a 2 fios paralelos. # Cabo 1 (fase) em x=0; cabo 2 (neutro) em x=d; correntes opostas (+I e -I). # Soma vectorial (Bx,Bz) e so depois calcula |B| = sqrt(Bx^2 + Bz^2). B_transect_uT <- function(x_m, d_m = 0.15, P_W = 3000, V_V = 230, pf = 1.0, h_m = 0.0065, # altura efectiva (m): 0.005 + isolamento (+ raio, se quiseres) B0_uT = 0) { # corrente RMS (aprox. resistivo se pf=1) I_rms <- P_W / (V_V * pf) # constante mu0/(2*pi) k <- 2e-7 # fios x1 <- 0 x2 <- d_m I1 <- I_rms I2 <- -I_rms # ponto do sensor: z = h_m z <- h_m dx1 <- x_m - x1 dx2 <- x_m - x2 r1_sq <- dx1^2 + z^2 r2_sq <- dx2^2 + z^2 # componentes do campo de cada fio (fio ao longo de y) Bx_T <- k * ( I1 * (-z) / r1_sq + I2 * (-z) / r2_sq ) Bz_T <- k * ( I1 * (dx1) / r1_sq + I2 * (dx2) / r2_sq ) # modulo e conversao para microtesla B_uT <- sqrt(Bx_T^2 + Bz_T^2) * 1e6 + B0_uT return(B_uT) } plot_B_transect <- function(n = 2001, d_m = 0.15, P_W = 3000, V_V = 230, pf = 1.0, h_m = 0.0065, B0_uT = 0) { # transecto: 15 cm antes do 1.o cabo, 15 cm entre, 15 cm depois do 2.o # cabo1 em 0; cabo2 em +15 cm => x em [-15cm, +30cm] x_m <- seq(-0.15, 0.30, length.out = n) B_uT <- B_transect_uT(x_m, d_m, P_W, V_V, pf, h_m, B0_uT) plot(x_m * 100, B_uT, type = "l", xlab = "Posicao no transecto (cm) [cabo 1: 0 cm; cabo 2: 15 cm]", ylab = "|B| (microtesla, RMS)", main = "|B| teorico (soma vectorial)") abline(v = 0, lty = 2) # cabo 1 abline(v = d_m * 100, lty = 2) # cabo 2 invisible(data.frame(x_m = x_m, x_cm = x_m * 100, B_uT = B_uT)) } # Exemplo: plot_B_transect() rad_juntos2<-read.table("~/Downloads/RADPi02.dat", header=TRUE, sep="\t", dec=".") X2 <- rad_juntos2[c("rep1","rep2","rep3")] rad_juntos2$rep_mean <- rowMeans(X2, na.rm=TRUE) rad_juntos2$dist<-sort(rad_juntos2$dist+10) rad_juntos3<-read.table("~/Downloads/RADPi03.dat", header=TRUE, sep="\t", dec=".") # Para já usar este: ## --- janela e eixo --- win <- subset(rad_juntos, dist >= 25 & dist <= 35) x_m <- win$dist/100 # cm -> m ## --- médias (3 replicados), em Tesla (assumindo nT) --- Bmed <- rowMeans(win[c("potmed1","potmed2","potmed3")], na.rm=TRUE) * 1e-9 Bmax <- rowMeans(win[c("potmax1","potmax2","potmax3")], na.rm=TRUE) * 1e-9 SDmed <- apply(win[c("potmed1","potmed2","potmed3")], 1, sd, na.rm=TRUE) * 1e-9 SDmax <- apply(win[c("potmax1","potmax2","potmax3")], 1, sd, na.rm=TRUE) * 1e-9 dat <- data.frame( x = c(x_m, x_m), B = c(Bmed, Bmax), sd = c(SDmed, SDmax), cond = factor(rep(c("med","max"), each=length(x_m))) ) ## --- correntes RMS fixas (a tua potência média por regime) --- V <- 230 Pmed <- 420.9333 Pmax <- 1388.667 Imed <- Pmed/V Imax <- Pmax/V dat$I <- ifelse(dat$cond=="med", Imed, Imax) ## --- pesos (replicados) --- dat$w <- 1/(dat$sd^2 + 1e-30) ## --- palpites --- x0_start <- x_m[which.max(Bmax)] B0_start <- min(dat$B) b2_start <- (0.005)^2 A_start <- max((max(Bmed)-B0_start) * b2_start / Imed, 0) fit_dip <- nls( B ~ B0 + (A*I)/((x - x0)^2 + b2), data = dat, weights = w, start = list(B0=B0_start, x0=x0_start, b2=b2_start, A=A_start), algorithm="port", lower = c(B0=-Inf, x0=min(x_m), b2=1e-10, A=0) ) print(summary(fit_dip)) p <- coef(fit_dip); print(p) ## --- previsões (nT) --- xgrid <- seq(min(x_m), max(x_m), length.out=400) Bfit_med <- (p["B0"] + (p["A"]*Imed)/((xgrid - p["x0"])^2 + p["b2"])) * 1e9 Bfit_max <- (p["B0"] + (p["A"]*Imax)/((xgrid - p["x0"])^2 + p["b2"])) * 1e9 ## --- plot --- o <- order(win$dist) ylim <- range(c(Bmed*1e9, Bmax*1e9, Bfit_med, Bfit_max), na.rm=TRUE) plot(win$dist[o], (Bmed*1e9)[o], type="n", xlab="dist (cm)", ylab="|B| (nT)", ylim=ylim, main="Cabos juntos: modelo dipolar físico (linhas) vs dados (pontos)") points(win$dist[o], (Bmed*1e9)[o], pch=17) points(win$dist[o], (Bmax*1e9)[o], pch=15) lines(xgrid*100, Bfit_med, lty=2) lines(xgrid*100, Bfit_max, lty=3) legend("topright", legend=c("med dados","max dados","med dipolo","max dipolo"), pch=c(17,15,NA,NA), lty=c(NA,NA,2,3), bty="n") ## --- verificação física: escala em I --- ratio_I <- Imax/Imed ratio_obs <- (max(Bmax)-p["B0"])/(max(Bmed)-p["B0"]) cat(sprintf("\nImax/Imed = %.4f ; (Bmax-B0)/(Bmed-B0) = %.4f\n", ratio_I, ratio_obs)) rad_juntos4<-read.table("~/Downloads/RADPi03bom.dat", header=TRUE, sep="\t", dec=".") X4 <- rad_juntos4[c("potmax1uT","potmax2uT","potmax3uT")] rad_juntos4$rep_mean <- rowMeans(X4, na.rm=TRUE) rad_juntos4$refnT <- (rad_juntos4$refnT)/1000 #plot(rad_juntos4$dist,rad_juntos4$rep_mean, ylim=c(0,100)) #lines(rad_juntos4$dist,rad_juntos4$refnT, lty=3) # ============================================================= # Ajuste do modelo teórico de dois cabos com correntes opostas # ============================================================= # Dados #dist <- c(0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20, # 21,22,23,24,25,26,27,28,29,30,31,32,33,34) dist <- rad_juntos4$dist #B_obs <- c(3.871333,4.580667,5.156000,5.959333,7.269333,8.642000, # 10.390000,13.003333,16.773333,22.403333,32.610000,46.770000, # 74.600000,80.946667,63.086667,49.223333,43.926667,40.680000, # 40.680000,42.953333,48.213333,59.820000,83.516667,71.906667, # 45.656667,30.950000,22.453333,16.353333,12.770000,10.100000, # 8.269333,6.908667,5.860333,5.008333,4.350333) B_obs <- rad_juntos4$rep_mean # Parâmetros físicos V <- 230 P <- 1388.667 I <- P / V mu0 <- 4 * pi * 1e-7 # ============================================================= # Modelo teórico # x1, x2 : posições dos cabos ao longo do transecto (cm) # h : altura do sensor acima dos cabos (cm) # Devolve |B| em µT # ============================================================= B_model <- function(x_cm, x1_cm, x2_cm, h_cm) { x <- x_cm / 100 x1 <- x1_cm / 100 x2 <- x2_cm / 100 h <- h_cm / 100 r1sq <- (x - x1)^2 + h^2 r2sq <- (x - x2)^2 + h^2 Bx <- (mu0 * I / (2 * pi)) * ( h / r1sq - h / r2sq) Bz <- (mu0 * I / (2 * pi)) * ((x - x1) / r1sq - (x - x2) / r2sq) sqrt(Bx^2 + Bz^2) * 1e6 } # ============================================================= # Função RSS — usa penalização em vez de Inf para L-BFGS-B # ============================================================= rss <- function(p) { x1 <- p[1]; x2 <- p[2]; h <- p[3] if (x2 - x1 < 0.5 || h <= 0) return(1e12) B_pred <- tryCatch( B_model(dist, x1, x2, h), error = function(e) rep(NA, length(dist)) ) if (any(is.na(B_pred) | !is.finite(B_pred))) return(1e12) sum((B_obs - B_pred)^2) } # ============================================================= # Grelha de valores iniciais (garante x2 > x1) # ============================================================= grid <- expand.grid( x1 = seq(10, 16, by = 1), x2 = seq(20, 26, by = 1), h = seq(0.5, 5, by = 0.5) ) grid <- grid[grid$x2 > grid$x1, ] rss_grid <- apply(grid, 1, function(p) rss(p)) best <- as.numeric(grid[which.min(rss_grid), ]) cat(sprintf("Melhor ponto inicial: x1=%.1f x2=%.1f h=%.1f RSS=%.2f\n", best[1], best[2], best[3], min(rss_grid))) # ============================================================= # Optimização com L-BFGS-B # ============================================================= fit <- optim( par = best, fn = rss, method = "L-BFGS-B", lower = c(0, 5, 0.1), upper = c(29, 34, 20.0), control = list(maxit = 10000, factr = 1e-14) ) cat(sprintf("Convergência: %s (código %d)\n", ifelse(fit$convergence == 0, "OK", "verificar"), fit$convergence)) # ============================================================= # Parâmetros ajustados # ============================================================= x1_fit <- fit$par[1] x2_fit <- fit$par[2] h_fit <- fit$par[3] cat("\n--- Parâmetros ajustados ---\n") cat(sprintf("Posição cabo 1 (x1): %.3f cm\n", x1_fit)) cat(sprintf("Posição cabo 2 (x2): %.3f cm\n", x2_fit)) cat(sprintf("Separação real: %.3f cm\n", x2_fit - x1_fit)) cat(sprintf("Altura sensor (h): %.3f cm\n", h_fit)) cat(sprintf("Corrente usada (I): %.4f A\n", I)) # ============================================================= # Qualidade do ajuste # ============================================================= B_pred <- B_model(dist, x1_fit, x2_fit, h_fit) ss_res <- sum((B_obs - B_pred)^2) ss_tot <- sum((B_obs - mean(B_obs))^2) r2 <- 1 - ss_res / ss_tot rmse <- sqrt(ss_res / length(B_obs)) cat(sprintf("R²: %.6f\n", r2)) cat(sprintf("RMSE: %.4f µT\n", rmse)) # ============================================================= # Gráfico # ============================================================= x_seq <- seq(0, 34, by = 0.1) B_seq <- B_model(x_seq, x1_fit, x2_fit, h_fit) plot(dist, B_obs, xlab = "distância (cm)", ylab = "|B| (µT)", main = sprintf("ajuste do modelo\nR² = %.4f | RMSE = %.3f µT | h = %.2f cm", r2, rmse, h_fit), ylim = c(0, max(B_obs) * 1.15), type="n") abline(v = x1_fit, lty = 3, col = "blue", lwd = 1.5) abline(v = x2_fit, lty = 3, col = "black", lwd = 1.5) lines(x_seq, B_seq, col = "grey", lwd = 2) lines(rad_juntos4$dist, rad_juntos4$refnT, lwd=2, lty=3, col="grey") points(dist, B_obs, pch = 21, bg="red") box(lwd=3) #legend("topright", # legend = c("Dados observados", # "Modelo teórico", # sprintf("Cabo 1 (x1 = %.1f cm)", x1_fit), # sprintf("Cabo 2 (x2 = %.1f cm)", x2_fit)), # col = c("steelblue", "firebrick", "darkgreen", "darkgreen"), # pch = c(19, NA, NA, NA), # lty = c(NA, 1, 2, 2), # lwd = 2) library(ggplot2) # ============================================================= # Dados experimentais (cabos a 10 cm de separação) # ============================================================= dados <- data.frame( d_cm = c(0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20, 21,22,23,24,25,26,27,28,29,30,31,32,33,34), B_uT = c(3.871333,4.580667,5.156000,5.959333,7.269333,8.642000, 10.390000,13.003333,16.773333,22.403333,32.610000,46.770000, 74.600000,80.946667,63.086667,49.223333,43.926667,40.680000, 40.680000,42.953333,48.213333,59.820000,83.516667,71.906667, 45.656667,30.950000,22.453333,16.353333,12.770000,10.100000, 8.269333,6.908667,5.860333,5.008333,4.350333) ) # Posições dos cabos ao longo do transecto (cm) d1 <- 12.424 d2 <- 22.537 # Altura do sensor acima da mesa (cm) h_sensor <- 1.52 # ============================================================= # Campo de um fio infinito no plano 2D # x : posição ao longo do transecto # y : altura acima/abaixo da mesa # ============================================================= B_wire <- function(x, y, x0, I = 1, k = 1, eps = 1.5) { dx <- x - x0 dy <- y r2 <- dx^2 + dy^2 + eps^2 Bx <- k * I * dy / r2 # componente ao longo do transecto By <- -k * I * dx / r2 # componente vertical list(Bx = Bx, By = By) } B_pair <- function(x, y, x1, x2, I = 1, k = 1, eps = 1.5) { b1 <- B_wire(x, y, x0 = x1, I = I, k = k, eps = eps) b2 <- B_wire(x, y, x0 = x2, I = -I, k = k, eps = eps) list(Bx = b1$Bx + b2$Bx, By = b1$By + b2$By) } # ============================================================= # Grelha 2D: # x : posição ao longo do transecto (0 a 34 cm) # y : altura acima/abaixo da mesa (-15 a 15 cm) # ============================================================= x_seq <- seq(0, 34, by = 1.5) y_seq <- seq(-15, 15, by = 1.5) grid <- expand.grid(x = x_seq, y = y_seq) b <- B_pair(grid$x, grid$y, x1 = d1, x2 = d2, I = 1, k = 1, eps = 1.5) grid$Bx <- b$Bx grid$By <- b$By grid$Bmag <- sqrt(grid$Bx^2 + grid$By^2) # ============================================================= # Escalar para µT usando os dados experimentais # O sensor estava em y = h_sensor # ============================================================= linha <- data.frame(x = dados$d_cm, y = h_sensor) bl <- B_pair(linha$x, linha$y, x1 = d1, x2 = d2, I = 1, k = 1, eps = 1.5) linha$Bmag_model <- sqrt(bl$Bx^2 + bl$By^2) scale_fac <- sum(linha$Bmag_model * dados$B_uT) / sum(linha$Bmag_model^2) grid$Bx_uT <- scale_fac * grid$Bx grid$By_uT <- scale_fac * grid$By grid$Bmag_uT <- scale_fac * grid$Bmag # Comprimento das setas arrow_fac <- 0.06 grid$xend <- grid$x + arrow_fac * grid$Bx_uT grid$yend <- grid$y + arrow_fac * grid$By_uT # ============================================================= # Objectos auxiliares # ============================================================= fios <- data.frame( x = c(d1, d2), y = 0, fio = c(sprintf("Cabo 1 (%.1f cm)", d1), sprintf("Cabo 2 (%.1f cm)", d2)) ) pts <- data.frame(x = dados$d_cm, y = h_sensor, B_uT = dados$B_uT) # ============================================================= # Gráfico # ============================================================= ggplot() + # Fundo cromático do módulo do campo geom_tile( data = grid, aes(x = x, y = y, fill = Bmag_uT), alpha = 1.0 ) + scale_fill_viridis_c(name = "|B| (µT)", option = "plasma") + # Setas vectoriais a preto geom_segment( data = grid, aes(x = x, y = y, xend = xend, yend = yend), arrow = arrow(length = unit(1.2, "mm"), type = "closed"), linewidth = 0.25, colour = "lightgrey", alpha = 0.8 ) + # Superfície da mesa (y = 0) geom_hline(yintercept = 0, linewidth = 0.8, linetype = "solid", colour = "tan4") + # Linha do transecto do sensor (y = h_sensor) geom_hline(yintercept = h_sensor, linewidth = 0.5, linetype = "dashed", colour = "steelblue") + # Cabos geom_point(data = fios, aes(x = x, y = y), size = 4, shape = 21, fill = "yellow", colour = "black") + geom_text(data = fios, aes(x = x, y = y - 1.5, label = fio), hjust = 0.5, size = 3, colour = "white") + # Pontos experimentais geom_point(data = pts, aes(x = x, y = y, size = B_uT), colour = "cyan", alpha = 0.85) + scale_size_continuous(name = "B medido (µT)", range = c(1, 7)) + # Anotações annotate("text", x = 1, y = 0.5, label = "Mesa", hjust = 0, colour = "tan4", size = 3) + annotate("text", x = 1, y = h_sensor + 0.8, label = sprintf("Transecto do sensor (y = %.2f cm)", h_sensor), hjust = 0, colour = "steelblue", size = 3) + coord_equal() + labs( x = "x (cm) — posição ao longo do transecto", y = "y (cm) — altura acima e abaixo da mesa", title = "Campo magnético vectorial — plano perpendicular aos cabos", subtitle = sprintf( "Cabos em x = %.1f cm e x = %.1f cm | Sensor em y = %.2f cm", d1, d2, h_sensor) ) + theme_bw(base_size = 11) + theme( legend.position = "right", panel.grid.major = element_line(colour = "grey90"), panel.grid.minor = element_blank(), panel.background = element_rect(fill = "white"), plot.background = element_rect(fill = "white") ) rad_juntos5<-read.table("~/Downloads/RADPi04.dat", header=TRUE, sep="\t", dec=".") X5 <- rad_juntos5[c("potmax1","potmax2","potmax3")] rad_juntos5$rep_mean <- rowMeans(X5, na.rm=TRUE) #plot(rad_juntos5$dist, rad_juntos5$rep_mean, ylim=c(0,100)) # ============================================================= # Ajuste do modelo teórico de dois cabos com correntes opostas # Cabos afastados 15 cm # ============================================================= # Dados #dist <- c(0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20, # 21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39) dist <- rad_juntos5$dist #B_obs <- c(4.664667,5.432000,6.082667,6.958667,8.567000,10.113333,11.833333,15.500000, # 19.746667,26.283333,37.986667,51.796667,87.556667,73.263333,52.333333,41.696667, # 36.030000,30.610000,28.206667,27.006667,26.786667,28.120000,30.766667,32.310000, # 37.630000,45.770000,67.250000,88.453333,55.243333,39.383333,28.913333,21.790000, # 15.253333,12.373333,10.040000,8.186667,6.927000,6.394333,5.557333,4.881000) B_obs <- rad_juntos5$rep_mean # Parâmetros físicos V <- 230 P <- 1388.667 I <- P / V mu0 <- 4 * pi * 1e-7 # ============================================================= # Modelo teórico # x1, x2 : posições dos cabos ao longo do transecto (cm) # h : altura do sensor acima dos cabos (cm) # Devolve |B| em µT # ============================================================= B_model <- function(x_cm, x1_cm, x2_cm, h_cm) { x <- x_cm / 100 x1 <- x1_cm / 100 x2 <- x2_cm / 100 h <- h_cm / 100 r1sq <- (x - x1)^2 + h^2 r2sq <- (x - x2)^2 + h^2 Bx <- (mu0 * I / (2 * pi)) * ( h / r1sq - h / r2sq) Bz <- (mu0 * I / (2 * pi)) * ((x - x1) / r1sq - (x - x2) / r2sq) sqrt(Bx^2 + Bz^2) * 1e6 } # ============================================================= # Função RSS — usa penalização em vez de Inf para L-BFGS-B # ============================================================= rss <- function(p) { x1 <- p[1]; x2 <- p[2]; h <- p[3] if (x2 - x1 < 0.5 || h <= 0) return(1e12) B_pred <- tryCatch( B_model(dist, x1, x2, h), error = function(e) rep(NA, length(dist)) ) if (any(is.na(B_pred) | !is.finite(B_pred))) return(1e12) sum((B_obs - B_pred)^2) } # ============================================================= # Grelha de valores iniciais (picos observados ~12 cm e ~27 cm) # ============================================================= grid <- expand.grid( x1 = seq(10, 16, by = 1), x2 = seq(24, 31, by = 1), h = seq(0.5, 5, by = 0.5) ) grid <- grid[grid$x2 > grid$x1, ] rss_grid <- apply(grid, 1, function(p) rss(p)) best <- as.numeric(grid[which.min(rss_grid), ]) cat(sprintf("Melhor ponto inicial: x1=%.1f x2=%.1f h=%.1f RSS=%.2f\n", best[1], best[2], best[3], min(rss_grid))) # ============================================================= # Optimização com L-BFGS-B # ============================================================= fit <- optim( par = best, fn = rss, method = "L-BFGS-B", lower = c(0, 5, 0.1), upper = c(34, 39, 20.0), control = list(maxit = 10000, factr = 1e-14) ) cat(sprintf("Convergência: %s (código %d)\n", ifelse(fit$convergence == 0, "OK", "verificar"), fit$convergence)) # ============================================================= # Parâmetros ajustados # ============================================================= x1_fit <- fit$par[1] x2_fit <- fit$par[2] h_fit <- fit$par[3] cat("\n--- Parâmetros ajustados (cabos a 15 cm) ---\n") cat(sprintf("Posição cabo 1 (x1): %.3f cm\n", x1_fit)) cat(sprintf("Posição cabo 2 (x2): %.3f cm\n", x2_fit)) cat(sprintf("Separação real: %.3f cm\n", x2_fit - x1_fit)) cat(sprintf("Altura sensor (h): %.3f cm\n", h_fit)) cat(sprintf("Corrente usada (I): %.4f A\n", I)) # ============================================================= # Qualidade do ajuste # ============================================================= B_pred <- B_model(dist, x1_fit, x2_fit, h_fit) ss_res <- sum((B_obs - B_pred)^2) ss_tot <- sum((B_obs - mean(B_obs))^2) r2 <- 1 - ss_res / ss_tot rmse <- sqrt(ss_res / length(B_obs)) cat(sprintf("R²: %.6f\n", r2)) cat(sprintf("RMSE: %.4f µT\n", rmse)) # ============================================================= # Gráfico # ============================================================= x_seq <- seq(0, 39, by = 0.1) B_seq <- B_model(x_seq, x1_fit, x2_fit, h_fit) plot(dist, B_obs, xlab = "distância (cm)", ylab = "|B| (µT)", main = sprintf("ajuste do modelo\nR² = %.4f | RMSE = %.3f µT | h = %.2f cm", r2, rmse, h_fit), ylim = c(0, max(B_obs) * 1.15), type = "n") abline(v = x1_fit, lty = 3, col = "blue", lwd = 1.5) abline(v = x2_fit, lty = 3, col = "black", lwd = 1.5) lines(x_seq, B_seq, col = "grey", lwd = 2) points(dist, B_obs, pch = 21, bg = "red") box(lwd = 3)