In chrono35/ArMag: ArMag

knitr::opts_chunk$set(echo = TRUE)

#shinyShortcut(shinyDirectory = "~/Documents/AM/R_Projects/ArMag/Visu_AM", OS = "unix",   gitIgnore = FALSE)

Memo

Quand on modifie le code, il faut penser à faire la documentation avec:

setwd("/Users/dufresne/Documents/R_Project/ArMag") devtools::document()

Chargement des fonctions d'ArMag via GitHub

if (!require(devtools))
 {install.packages("devtools")}

devtools::install_github("chrono35/ArMag", force = TRUE)
library("ArMag") #, lib.loc="/Library/Frameworks/R.framework/Versions/3.5/Resources/library")

Chargement des fonctions d'ArMag via source

Si besoin

# source('R/ArMag.R', echo=FALSE)

Chargement fichier test avec fonction read.Pal

Lecture fichier pal

#file.AM <- "~/Documents/AM/GDRE1_INT.AMD"
# file.AM <- "/Volumes/Partage_Mesures_Archeomag/Mesures/22061_Glomel/22061A.AMD"
# file.AM <- "~/Documents/AM/22061_Glomel/22061A.AMD"
#file.AM <- "~/Documents/AM/53130_Laval/53130_Mesures/53130B1_2019.AMD"
#file.Pal <- "~/Documents/AM/Agnes&Maria/Portugal.txt"
file.Pal <- "examples/Portugal.txt"

mes <- NULL
mes <- read.Pal.mesures(file.Pal)
mes.info <- read.Pal.info (file.Pal)

demag(mes)

# export en csv
#write.csv(mes, file = "~/Documents/AM/Agnes&Maria/Portugal.csv")

# calcul de l'anisotropie moyenne
ani.moyen <- anisotropy.mean.eigen.tensor(mes, mes.info$number[1:20], step.value = 480, step.code = c("Z+", "Z-", "X+", "X-", "Y+", "Y-", "ZB"))

# calcul et stockage des anisotropies individuelles
ani.all.spe <- anisotropy.eigen.tensors.all(mes, mes.info$number[1:20], step.value = 480, step.code = c("Z+", "Z-", "X+", "X-", "Y+", "Y-", "ZB"))
ani.all.spe.mane <- cbind(mes.info$name[1:20], ani.all.spe)

#png(paste(TeX.floder,"22061A_ANI.png", sep = "") , bg = "transparent", width = 1000) # Directement dans le dossier TeX
#par(mfrow = c(1, 2)) # pour exportation

par(mfrow = c(1, 2), cex.axis = 0.5)

lambert.ID.tensors(ani.all.spe)
lambert.ID.tensors(ani.moyen, pt.col = "red", new = FALSE)

flinn(ani.all.spe)
flinn(ani.moyen, pt.col = "red", new = FALSE)

# dev.off() # fin sauvegarde-exportation

Chargement fichier test avec fonction read.AM

Test execution

# file.AM <- "~/Documents/AM/GDRE1_INT.AMD"
# file.AM <- "/Volumes/Partage_Mesures_Archeomag/Mesures/22061_Glomel/22061A.AMD"
# file.AM <- "~/Documents/AM/22061_Glomel/22061A.AMD"
 file.AM <- "examples/14039C.AMP"
# file.AM <- "~/Documents/AM/Aire/40001A.AMD"

mes <- NULL
mes <- read.AM.mesures(file.AM)
mes.info <- read.AM.info (file.AM)

select <- NULL
select$I <- mes$I[mes$step.value== 0]
select$D <- mes$I[mes$step.value== 0]

stat.mcFadden(select$I[1:8], select$D[1:8])


#lambert(mes, inc.lim = c(0,90))

Codage correction anisotropie

# calcul les nouveaux paramètres d'un echantillon avec un tenseur donnée
#procedure Calcul_Aniso_Specimen(var R_Cal_echt:Res_Cal_xyz;Tens_ani:TENSEUR_ANISO );
mes.sel <- NULL
mes.sel1 <- extract.mesures.specimen.name("14P2", mes)

volume <- mes.info$vol[which(mes.info$name == trimws(mes.sel1$name[1]))]
TH <- mes.info$TH[which(mes.info$name == trimws(mes.sel1$name[1]))]
mes.sel1$X <- mes.sel1$X/volume * 1E6
mes.sel1$Y <- mes.sel1$Y/volume * 1E6
mes.sel1$Z <- mes.sel1$Z/volume * 1E6
ani1 <- anisotropy.matrix.symetric(mes.sel1, step.value = 550, TH = TH, volume = volume )

mes.sel2 <- extract.mesures.specimen.name("16P1", mes)
volume <- mes.info$vol[which(mes.info$name == trimws(mes.sel2$name[1]))]
TH <- mes.info$TH[which(mes.info$name == trimws(mes.sel2$name[1]))]
ani2 <- anisotropy.matrix.symetric(mes.sel2, step.value = 550, TH = TH, volume = volume )

mes.sel3 <- extract.mesures.specimen.name("24P1", mes)
volume <- mes.info$vol[which(mes.info$name == trimws(mes.sel3$name[1]))]
TH <- mes.info$TH[which(mes.info$name == trimws(mes.sel3$name[1]))]
ani3 <- anisotropy.matrix.symetric(mes.sel3, step.value = 550, TH = TH, volume = volume )


mes.sel4 <- extract.mesures.specimen.name("15P1", mes)
volume <- mes.info$vol[which(mes.info$name == trimws(mes.sel4$name[1]))]
TH <- mes.info$TH[which(mes.info$name == trimws(mes.sel4$name[1]))]
ani4 <- anisotropy.matrix.symetric(mes.sel4, step.value = 550, TH = TH, volume = volume )


mes.sel.ARN <- mes[mes$step.value == 0,]
mes.sel.ARN$X <- mes.sel.ARN$X /volume * 1E6
mes.sel.ARN$Y <- mes.sel.ARN$Y/volume * 1E6
mes.sel.ARN$Z <- mes.sel.ARN$Z/volume * 1E6

ani <- anisotropy.matrix.symetric(mes, step.value = 550)
ani.tens <- anisotropy.eigen.tensor(mes, step.value = 550 ) # eigen.tensor

# matrice moyenne
  ani.mean <- (ani2+ ani3 + ani4)/3
  eigen(ani.mean)

test correction.anisotropie

mes.sel <- NULL
mes.sel <- extract.mesures.specimen.name("24P1", mes)
volume <- mes.info$vol[which(mes.info$name == trimws(mes.sel$name[1]))]
TH <- mes.info$TH[which(mes.info$name == trimws(mes.sel$name[1]))]

mes.sel.ARN <- mes.sel[mes.sel$step.value == 0,]
mes.sel.ARN$X <- mes.sel.ARN$X /volume * 1E6
mes.sel.ARN$Y <- mes.sel.ARN$Y /volume * 1E6
mes.sel.ARN$Z <- mes.sel.ARN$Z /volume * 1E6

ani <- anisotropy.matrix.symetric(mes.sel, step.value = 550, TH = TH, volume = volume )

# normalisation des mesures en A/m
mes.sel$X <- mes.sel$X /volume * 1E6
mes.sel$Y <- mes.sel$Y /volume * 1E6
mes.sel$Z <- mes.sel$Z /volume * 1E6
correction.anisotropy(mes.sel, ani)

Tracer

mes.ech <- NULL
mes.ech <- extract.mesures.specimen.name("10P1", mes)
  mes.ech <- remove.step(mes.ech, step.value = 550, step.code = c("??"))
   mes.ech <- remove.step(mes.ech, step.value = 400, step.code = c("??"))
#mes.sel2 <- extract.mesures.specimen.number(2, mes2)

Tracer de Zijderveld

par(pty="s") # force une figure carré
zijderveld1(mes.ech$X, mes.ech$Y, mes.ech$Z) #, pt.name = mes.ech$Etap)
par(mfrow = c(1,2), pty="m", cex.lab = 0.5, cex.axis = 0.6) # séparation en 2 colonnes et rétablie une figure taille maximale
# cex.axis défini la taille du texte des échelles
# cex.lab défini la taille du texte des étapes 
zijderveld1(mes.ech, legend.pos = "topright")
zijderveld2(mes.ech, pt.names = NULL)

Nous devons obtenir la même chose. Chacune des fonction utilisant des paramètres différents

lambert.XYZ.specimen(mes.ech)
lambert.ID.specimen(mes.ech)

Calcul sur matrice

volume <- mes.info$vol[which(mes.info$name == trimws(mes.sel$name[1]))]
TH <- mes.info$TH[which(mes.info$name == trimws(mes.sel$name[1]))]

ani <- anisotropy.matrix.symetric(mes.sel, step.value = 550, TH = TH, volume = volume )
ani

TH <- 60
mes.sel <- NULL

mes.sel1 <- extract.mesures.specimen.number(1, mes)
mes.sel2 <- extract.mesures.specimen.number(2, mes)
mes.sel3 <- extract.mesures.specimen.number(3, mes)

ani1 <- anisotropy.eigen.tensor(mes.sel1, step.value = 550)
ani2 <- anisotropy.eigen.tensor(mes.sel2, step.value = 550 )
ani3 <- anisotropy.eigen.tensor(mes.sel3, step.value = 550 )

Data <- rbind(ani1, ani2, ani3)
lambert.ID.tensors(Data, point.col = c("blue", "red", "green"))

numbers <- c(7)
Data <- anisotropy.eigen.tensors.numbers(numbers, mes, mes.info, 400)
lambert.ID.tensors(Data, point.col = c("blue", "red", "green"))


numbers.selected <- NULL
  for (i in 1: length(mes$number) ) {
    if ( length(which( numbers.selected == mes$number[i] )) == 0 )
      numbers.selected <- c(numbers.selected, mes$number[i])
  }

flinn(Data)#, points.col = c("red", "blue"), pt.names = c("titi", "toto"))

mes.seltt <- NULL
mes.seltt <- rbind(mes.sel1, mes.sel2)

mes.sel2 <- extract.mesures.specimen.name("1P1", mes)
mes.sel3 <- extract.mesures.specimen.name("16P1", mes)
demag(mes.sel2, etap.J0 = NULL)
demag(mes, etap.J0 = 0, pt.col = rainbow(length(mes.info$name)))
demag(rbind(mes.sel2, mes.sel3), normalize = TRUE)

partial.component(mes.sel2$X, mes.sel2$Y, mes.sel2$Z)

par(mfrow = c(2, 2), cex.lab = 0.7, cex.axis = .7, cex = 0.7, cex.main = 1, cex.sub = 0.1,
    mai = c(0.5, 0.5, 0.7, 0.3), oma = c(0, 1, 1, 1))#, pty ="s" )
zijderveld1.T1T2(mes.sel2)
zijderveld2.T1T2(mes.sel2)

# suppression des étapes d'anisotropie
mes.sel2 <- remove.step(mes.sel2) 
lambert(mes.sel2, inc.lim = c(0, 90))
demag(mes.sel2, etap.J0 = NULL)

file.AM <- "examples/test.txt"
list.ech <- c("1T",   "17T",  "35T",  "47T",  "50T",  "59T",  "83T",  "87T",  "89T",  "94T",  "100T", "102T", "103T" )
genere.AMD(file.AM, list.ech)

Calcul vecteur partiel

sp.res <- NULL
for (sp in 1: length(mes.info$number)){
  sp.mes <- extract.mesures.specimen.number(mes.info$number[sp], mes)


  # correction d'anisotropie
  volume <- mes.info$vol[which(mes.info$name == trimws(sp.mes$name[1]))]
  TH <- mes.info$TH[which(mes.info$name == trimws(sp.mes$name[1]))]

  # normalisation des mesures en A/m
  sp.mes$X <- sp.mes$X /volume * 1E6
  sp.mes$Y <- sp.mes$Y /volume * 1E6
  sp.mes$Z <- sp.mes$Z /volume * 1E6

  # récupération de l'ARN
  #sp.mes.ARN <- sp.mes[sp.mes$step.value == 0,]


  # évaluation de l'anisotropie
  sp.ani <- anisotropy.matrix.symetric(sp.mes, step.value = 550, TH = TH, volume = volume )
 # controle de l'existance de l'anisotropie
  if (!is.na(sum(sp.ani))) {
    # suppression des étapes d'anisotropie
    sp.mes <- remove.step(sp.mes) #, step.value = "0", step.code = "N") 

    # correction de l'anisotropie  
    sp.cor <- correction.anisotropy(sp.mes, sp.ani)

  } else {
    warning("No anisotropy step \n")
    sp.cor <- sp.mes
  }

  # determination de la composante
  sp.cp <- partial.component(sp.cor$X, sp.cor$Y, sp.cor$Z)

  # repliement
  sp.cp <- repliement.auto(sp.cp$I, sp.cp$D, name = sp.mes$name)
  sp.res <- rbind.data.frame(sp.res, cbind(name =mes.info$name[sp], sp.cp) )
}


stat.mcFadden(sp.res$I, sp.res$D)


lambert(sp.res, inc.lim = c(0,90))

Calcul azimut du soleil

JulianDay <- function(day, month, year, UT=12)
  # valid from 1900/3/1 to 2100/2/28 # http://jgiesen.de/elevazmoon/basics/meeus.htm
  # equivalent à la fonction dans R : julian(Sys.Date(), -2440588)
  # ISOdate(year, month, day, hour = 12, min = 0, sec = 0, tz = "GMT")
{
  #  if (month<=2) {
   #   month=month+12
   #   year=year-1
   #   }
    #return ( trunc(365.25*year) + trunc(30.6001*(month+1)) - 15 + 1720996.5 + day + UT/24.0 )
  return ( 367*year-7*trunc( (year+trunc((month+9)/12 ) )/4)+275*trunc(month/9)+day+1721014 ) # https://cral-perso.univ-lyon1.fr/labo/fc/cdroms/docu_astro/jour_julien/jour_julien.pdf
} 

# Le jour julien 0 commence le 24 novembre -4713 (4712 BC) à 12h
#' Le nombre de jour julien pour les calculs astronomiques
#' @seealso \code{\link{https://codes-sources.commentcamarche.net/source/31774-calcul-de-la-position-du-soleil-declinaison-angle-horaire-altitude-et-azimut-altaz-solaire}}
#' @export
jour.julien <- function( jour,   mois,   annee,   heure,   minute,   seconde)
{
        day <- jour + heure/24 + minute/1440 + seconde/86400
        year <- annee
        month <- mois

        if (month == 1 || month == 2)
        {
                year <- year-1
                month <- month+12
        }

        a <- trunc(year/100)
        b <- 2 - a + trunc(a/4)

        jour_julien <- trunc(365.25*(year+4716)) + trunc(30.6001*(month+1)) + day + b - 1524.5

        return (as.numeric(jour_julien))
}

juju <- function( jour,   mois,   annee)
{
            jour_julien <- (1461 * (annee + 4800 + (mois - 14)/12))/4 +(367 * (mois - 2 - 12 * ((mois - 14)/12)))/12- (3 * ((annee + 4900 + (mois - 14)/12)/100))/4 + jour - 32075.0
        return (as.numeric(jour_julien) )
}

# https://codes-sources.commentcamarche.net/source/31774-calcul-de-la-position-du-soleil-declinaison-angle-horaire-altitude-et-azimut-altaz-solaire

# Pour controler
# https://www.sunearthtools.com/dp/tools/pos_sun.php?lang=fr
# https://fr.planetcalc.com/320/

jour2= 28 ; mois2= 04; annee2= 1994; heure2= 12; minute2= 51; seconde2= 10
longitude= 20; longmin= 20; longsec= 46    # Est est positif, et Ouest est négatif
latitude= 50; latmin= 41; latsec= 53; deviation = 0
tz = "UTC"
sun.azimuth(jour2, mois2, annee2, heure2, minute2, seconde2, longitude, longmin, longsec, latitude, latmin, latsec)
360-sun.azimuth(jour2, mois2, annee2, heure2, minute2, seconde2, longitude, longmin, longsec, latitude, latmin, latsec)

sun.azimuth (10, 11, 2019, 11, 19, seconde=0, 48, longmin=0, longsec=0, 45, latmin=0, latsec=0)

test igrf13syn

isv <- 1
date <- 2000
itype <- 1
alt <- 50
colat <- 90 - 45
elong <- 10

igrf13syn(isv=isv, date=date, itype=itype , alt=alt, colat= colat, elong=elong)

Vieux fichier

file.AM <- "examples/old-type.AMD"
oldType.file.info <- read.AM.oldType.info(file.AM, encoding = "macroman")
oldType.file <- read.AM.oldType.mesures(file.AM, encoding = "macroman")

test codage diagramme ARAI

file.INT <- "examples/INT_example.AMD"

mesINT <- NULL
mesINT <- read.AM.mesures(file.INT)
mesINT.info <- read.AM.info (file.INT)
mes.sel<- extract.mesures.specimen.name("40001B_11B1", mesINT)

# reference 
# Coe 1978 : DOI: 10.1029/JB083iB04p01740
# Prévost et Al. 1985 DOI: 10.1029/JB090iB12p10417
par(pty="s", "xaxp")
relative = FALSE
verbose = TRUE
show.plot = TRUE
vol=10.8
TH = 60
aim.coef = 1E-10*1E6/vol #1E6
show.step.value = FALSE
R.mark = 'R'  # Positive pTRM
V.mark = 'V'  # Negative pTRM
P.mark = 'P'  # pTRM check
L.mark = "L"  # sLow cooling 
Q.mark = "Q"  # Quick cooling
pt.col = "blue"
loop.col = "forestgreen"
step.J0 = "20N0" # ou NULL

begin.step.value = 0
end.step.value = 700

par(pty="s", "xaxp")
arai(mes.sel, begin.step.value = 250, end.step.value = 700, aim.coef = 1E-10*1E6/vol)

chrono35/ArMag documentation built on June 19, 2024, 6:38 p.m.