diff --git a/maglib/src/dconjrv.f b/maglib/src/dconjrv.f new file mode 100755 index 0000000..77d47fb --- /dev/null +++ b/maglib/src/dconjrv.f @@ -0,0 +1,420 @@ + subroutine dconjrv (magin,year,magout,indgm,indval,tilt,rb, + > rggsm,rgsmg,dir,rre,theto, + > phio,rfin,rmax,np,tr,tthet,tphi,ier) +c* +c*********************************************************************** +c* +c* "Copyright [c] CNES 98 - tous droits reserves" +c* ********************************************** +c* +c*PRO MAGLIB +c* +c*VER 99.03.31 - V 1.0 +c*VER 01.05.30 - V 2.0 +c*VER 03.01.06 - V 2.1 +c*VER 01.01.07 - V 3.0 +c*VER 13.10.10 - V 4.0 +c*VER 17.02.23 - V 5.0 +c*VER 20.07.15 - V 6.0 +c* +c*AUT spec. CNES - JC KOSIK - juin 1995 +c*AUT port. CISI +c*AUT adapt. AKKA +c* +c*ROL Theme : Calculs de geophysique +c*ROL Calcul du point conjugue d'un point donne en tenant +c*ROL compte d'une combinaison de champs interns et externe. +c* +c*PAR magin (I) : type de champ magnetique interne +c* +c*PAR year (I) : annee fractionnaire +c* +c*PAR magout (I) : type de champ magnetique externe +c* +c*PAR indgm (I) : type d'indice geomagnetique Kp ou Ae +c*PAR indval (I) : indice geomagnetique de 1 a 6 +c* +c*PAR tilt (I) : angle de tilt (radians) +c* +c*PAR rb (I) : distance subsolaire normalisee de la magnetopause +c*PAR : (rayons terrestres) rb = 10. * re +c* +c*PAR rggsm (I) : matrice (3,3) de passage du repere geographique au +c*PAR : repere solaire magnetospherique +c*PAR rgsmg (I) : matrice (3,3) de passage du repere solaire +c*PAR : magnetospherique au repere geographique +c* +c*PAR dir (I) : direction du trace des lignes de champ +c*PAR : dir = -1 trace vers le meme hemisphere +c*PAR : dir = +1 trace vers l'hemisphere oppose +c* +c*PAR rre (I) : distance radiale geocentrique (rayons terrestres) +c*PAR theto (I) : colatitude geocentrique (radians) +c*PAR phio (I) : longitude geocentrique (radians) +c* +c*PAR rfin (I) : distance geocentrique d'arret (rayons terrestres) +c* +c*PAR rmax (I) : distance geocentrique maximum autorisee +c*PAR : (rayons terrestres) +c* +c*PAR np (O) : nombre de points calcules +c* +c*PAR tr (O) : tableau des distances geocentriques des points de la +c*PAR : ligne de champ (rayons terrestres) +c*PAR tthet (O) : tableau des colatitudes des points de la +c*PAR : ligne de champ (radians) +c*PAR tphi (O) : tableau des longitudes des points de la +c*PAR : ligne de champ (radians) +c* +c*PAR ier (O) : code de retour +c* +c*NOT magin : 1 = dipole + g11 + h11 +c*NOT magin : 2 = champ DGRF 1945-1970 +c*NOT magin : 3 = champ DGRF 1970-1995 +c*NOT magin : 3 = champ DGRF 1995-2000 +c* +c*NOT year : year doit etre >= 1945. +c* +c*NOT magout : 0 = pas de champ externe +c*NOT magout : 1 = Tsyganenko 87 +c*NOT magout : 2 = Tsyganenko 89 +c*NOT magout : 3 = Kosik 97 +c* +c*NOT indgm : 1 ---> indices geomagnetique Kp +c*NOT indval : 1 ---> Kp = 0 , 0+ +c*NOT indval : 2 ---> Kp = 1- , 1 , 1+ +c*NOT indval : 3 ---> Kp = 2- , 2 , 2+ +c*NOT indval : 4 ---> Kp = 3- , 3 , 3+ +c*NOT indval : 5 ---> Kp = 4- , 4 , 4+ +c*NOT indval : 6 ---> Kp > 5- +c* +c*NOT indgm : 2 ---> geomagnetique Ae +c*NOT indval : 1 ---> Ae = 0 - 50 +c*NOT indval : 2 ---> Ae = 50 - 100 +c*NOT indval : 3 ---> Ae = 100 - 150 +c*NOT indval : 4 ---> Ae = 150 - 250 +c*NOT indval : 5 ---> Ae = 250 - 400 +c*NOT indval : 6 ---> Ae >= 400 +c* +c*NOT dir : +1 = trace vers les altitudes plus elevees +c*NOT : (depuis la surface vers l hemisphere oppose) +c*NOT : -1 = trace vers les altitudes plus basses +c*NOT : (vers la surface du meme hemisphere +c* +c*NOT rfin : 1 = dernier point sur Terre +c* +c*NOT rmax : point d arret des lignes de champ ouvertes +c* +c*NOT ier : 1 = l > 500 ou r > rmax ou (rre-1) <= 0.01 +c*NOT : 2 = sortie de magnetopause +c* +c*NOT common : util +c* +c*INF utilise : magtot, spcar, mpause, step, carsp +c* +c*HST version 1.0 - 99.03.31 - creation de la maglib au CDPP +c*HST version 2.0 - 01.05.30 - Ajout indice Kosik 97 +c*HST version 2.1 - 03.01.06 - corrections en compilation avec g77 +c*HST version 3.0 - 01.01.07 - mise a jour des modeles de champ interne +c*HST version 4.0 - 13.10.10 - mise a jour des modeles de champ interne +c*HST (DGRF05 et IGRF10) +c*HST version 5.0 - 17.02.23 - mise a jour des modeles de champ interne +c*HST (DGRF10 et IGRF15) +c*HST version 5.0 - 20.07.15 - mise a jour des modeles de champ interne +c*HST (DGRF15 et IGRF20) +c* +c*********************************************************************** +c* + implicit none +c +c --------------------------------- +c*FON Declaration identificateur rcs_id +c --------------------------------- +c + character rcs_id*100 +c +c -------------------------- +c*FON Declaration des parametres +c -------------------------- +c + integer magin + double precision year + integer magout, indgm, indval + double precision tilt, rb + double precision rggsm(3,3), rgsmg(3,3) + double precision dir + double precision rre, theto, phio + double precision rfin, rmax + integer np + double precision tr(500), tthet(500), tphi(500) + integer ier +c +c ---------------------------------- +c*FON Declaration des variables communes +c ---------------------------------- +c + double precision pi,dpi,rad,deg,pid,xmu,rayt +c +c*COM pi : constante pi (obtenue a partir de acos(-1.)) +c*COM dpi : constante 2 * pi +c*COM pid : constante pi / 2 +c*COM rad : facteur de conversion degres ----> radians +c*COM deg : facteur de conversion radians ----> degres +c*COM xmu : constante de gravitation terrestre (km**3/sec**2) +c*COM rayt : rayon equatorial terrestre (km) +c + common/util/pi,dpi,rad,deg,pid,xmu,rayt +c +c --------------------------------- +c*FON Declaration des variables locales +c --------------------------------- +c + integer ipop +c*LOC ipop : indicateur de positionnement par rapport a la magnetosphere +c* + integer ier1,ier2,ier3,ier4,ier5,ier6,ier7 +c*LOC ier1,ier2,ier3,ier4,ier5,ier6,ier7 : codes retour des modules appeles +c + double precision br,bt,bp +c*LOC br,bt,bp : composantes radiale, tangentielle et azimuthale du champ +c + double precision x1,x2,x3 + double precision xxpp,xxp,xx,yypp,yyp,yy,zzpp + double precision zzp,zz +c*LOC x1,x2,x3 : coordonnees cartesiennes x, y et z du point courant +c + double precision rr,thetr,phir +c*LOC rr,thetr,phir : coordonnees spheriques +c + double precision errin +c*LOC errin: erreur maximale permise lors de l'integration +c + double precision tx(500),ty(500),tz(500) +c*LOC tx,ty,tz : tableau des coordonnees cartesiennes x, y et z +c + double precision a1,a2,a3,a4,a5,a6,a7 +c*LOC a1,a2,a3,a4,a5,a6,a7 : termes des equations +c + double precision altgc +c*LOC altgc : altitude en km du 1er point de la ligne de champ +c + double precision epsr +c*LOC epsr : epsilon pour validite de l'altitude +c + double precision ds0,rpp,r,b,sgn,ds,rp +c*LOC ds0,rpp,r,b,sgn,ds,rp : variables pour la determination du ds +c + double precision px,py,pz +c*LOC px,py,pz : coordonees cartesiennes +c +c ----------------------------------------- +c*FON Declaration de la fonction interne pinter +c ----------------------------------------- +c + double precision pinter +c + SAVE +c +c ---------------------------------------- +c*FON Definition de la fonction interne pinter +c ---------------------------------------- +c + pinter(a1,a2,a3,a4,a5,a6,a7)= + > ((a2 - a3) * (a7 - a2) * (a7 - a3) * a4 + > - (a1 - a3) * (a7 - a1) * (a7 - a3) * a5 + > + (a1 - a2) * (a7 - a2) * (a7 - a1) * a6) + > / ((a1 - a2) * (a2 - a3) * (a1 - a3)) +c +c --------------------------------- +c*FON Affectation identificateur rcs_id +c --------------------------------- +c + data rcs_id /" + >$Id$"/ +c +c ****************** +c Debut de programme +c ****************** +c +c --------------- +c*FON Initialisations +c --------------- +c + ier = 0 + ier1 = 0 + ier2 = 0 + ier3 = 0 + ier4 = 0 + ier5 = 0 + ier6 = 0 + ier7 = 0 + ipop = 0 +c + rr = rre + thetr = theto + phir = phio +c +c ------------------- +c*FON Determination du ds +c ------------------- +c + altgc = (rr - 1.0d0) + epsr = 0.001d0 + ds0 = 0.5d0 + if (dir .lt. 0.0d0) then + if (altgc .gt. epsr .and. altgc .lt. 2.0d0) then + ds0 = altgc / 20.0d0 + else if (altgc .le. epsr) then + ier = 0 + np = 1 + tr(1) = rr + tthet(1) = thetr + tphi(1) = phir + endif + endif + ds = ds0 +c +c --------------------------------------------------------------- +c*FON Hemisphere Nord Br < 0, dir < 0, sgn > 0, trace dans le sens du +c*FON champ, donc vers l'hemisphere Nord. Si dir > 0, sgn <0, trace +c*FON dans le sens oppose au champ, donc vers l'hemisphere Sud. +c*FON Hemisphere Sud : Br > 0, dir < 0, sgn < 0, trace dans le sens +c*FON oppose au champ, donc vers l'hemisphere Sud. Si dir > 0, sgn > 0, +c*FON trace dans le sens du champ, donc vers l'hemisphere Nord. +c --------------------------------------------------------------- +c + call magtotv(magin,year,magout,indgm,indval,tilt,rggsm,rgsmg,rr, + > thetr,phir,br,bt,bp,b,ier1) +c + errin = 0.0005d0 + sgn = sign(1.d0, br * dir) + ds = ds0 * sgn + np = 1 +c +c ---------------------------------------------------- +c*FON Calcul de la position cartesienne du point de depart +c ---------------------------------------------------- +c + call spcar(rr,thetr,phir,x1,x2,x3,ier2) +c + tr(1) = rr + tthet(1) = thetr + tphi(1) = phir +c + tx(np) = x1 + ty(np) = x2 + tz(np) = x3 +c +c ------------------------------------------------------------------- +c*FON Calcul et test d'appartenance a la magnetosphere de type Shabansky, +c*FON si ipop = 2 point en dehors ipop = 0, point a l'interieur +c*FON Si le point de depart est a l'exterieur on ne fait pas les calculs +c*FON et on met 999 dans les valeurs en 25 +c ----------------------------------------------- +c + call mpause(rb,rggsm,rr,thetr,phir,ipop,ier3) +c + if (ipop .eq. 2) go to 10 +c +40 continue + np = np + 1 +c +c ------------------------------------------------------- +c*FON Test de sortie de la magnetosphere, si ipop=2 on +c*FON arrete les calculs sans mettre 999 aux valeurs +c*FON et on se debranche sur le return, on a trace un bout de +c*FON ligne de champ c'est tout +c ------------------------------------------------------- +c + call mpause(rb,rggsm,rr,thetr,phir,ipop,ier4) +c + if (ipop .eq. 2) go to 30 +c +c -------------------------------------------------- +c*FON Integration pas a pas par appels successifs a step +c -------------------------------------------------- +c + call stepv(magin,magout,indgm,indval,year,tilt, + > rggsm,rgsmg,x1,x2,x3,ds,errin,ier5) +c + tx(np) = x1 + ty(np) = x2 + tz(np) = x3 +c +c ------------------------------------------------ +c*FON Conversion des coordonnees cartesiennes du point +c*FON en coordonnees spheriques +c ------------------------------------------------ +c + call carsp(x1,x2,x3,rr,thetr,phir,ier6) +c + tr(np) = rr + tthet(np) = thetr + tphi(np) = phir +c + if (rr .gt. rfin .and. np .lt. 500 .and. rr .lt. rmax) go to 40 +c + if (np .ge. 500 .or. rr .ge. rmax) then + ier = 1 + tr(np) = 999.0d0 + tthet(np) = 999.0d0 * rad + tphi(np) = 999.0d0 * rad + endif +c +c ---------------------------------------------------------- +c*FON La limite rfin est depassee, on interpole pour obtenir les +c*FON coordonnees exactes du point en cartesiennes converties +c*FON en coordonnees spheriques apres calcul +c ---------------------------------------------------------- +c + if (rr .lt. rfin) then +c + rpp = tr(np-2) + rp = tr(np-1) + r = tr(np) + xxpp = tx(np-2) + xxp = tx(np-1) + xx = tx(np) + yypp = ty(np-2) + yyp = ty(np-1) + yy = ty(np) + zzpp = tz(np-2) + zzp = tz(np-1) + zz = tz(np) + px = pinter(rpp,rp,r,xxpp,xxp,xx,rfin) + py = pinter(rpp,rp,r,yypp,yyp,yy,rfin) + pz = pinter(rpp,rp,r,zzpp,zzp,zz,rfin) +c + call carsp (px,py,pz,rr,thetr,phir,ier7) +c + tr(np) = rr + tthet(np) = thetr + tphi(np) = phir +c + endif +c +c -------------------------------------------------------------- +c*FON Si la distance limite rmax ou si le nbr de points depasse 500 +c*FON ou si on arrive en dehors de la magnetosphere de Shabansky on +c*FON remplit les coordonnes par des 999 ou 999*rad pour obtenir des +c*FON 999 en degres apres la conversion radians degres +c -------------------------------------------------------------- +c +10 continue +c + if (ipop .eq. 2) then +c + tr(np) = 999.0d0 + tthet(np) = 999.0d0 * rad + tphi(np) = 999.0d0 * rad + ier = 2 +c + endif +30 continue +c +c **************** +c Fin de programme +c **************** +c + return + end diff --git a/maglib/src/derivv.f b/maglib/src/derivv.f new file mode 100755 index 0000000..560b70f --- /dev/null +++ b/maglib/src/derivv.f @@ -0,0 +1,228 @@ + subroutine derivv (magin,magout,indgm,indval,year,tilt, + > rggsm,rgsmg,ds3,xgc,ygc,zgc,rxp1,rxp2, + > rxp3,ier) +c* +c*********************************************************************** +c* +c* "Copyright [c] CNES 98 - tous droits reserves" +c* ********************************************** +c* +c*pro MAGLIB +c* +c*VER 99.03.31 - V 1.0 +c*VER 01.05.30 - V 2.0 +c*VER 03.01.06 - V 2.1 +c*VER 01.01.07 - V 3.0 +c*VER 13.10.10 - V 4.0 +c*VER 17.02.23 - V 5.0 +c*VER 20.07.15 - V 6.0 +c* +c*AUT spec. CNES - JC KOSIK - avril 1996 +c*AUT port. CISI +c*AUT adpat. AKKA +c* +c*ROL Theme : Calculs de geophysique +c*ROL Calcul des pas elementaires rxp1,rxp2,rxp3 au point +c*ROL courant xgc,ygc,zgc. +c* +c*PAR magin (I) : type de champ magnetique interne +c*PAR magout (I) : type de champ magnetique externe +c* +c*PAR indgm (I) : type d'indice geomagnetique Kp ou Ae +c*PAR indval (I) : indice geomagnetique : niveau d'amplitude du champ +c* +c*PAR year (I) : annee fractionnaire (>= 2000.) +c* +c*PAR tilt (I) : angle de tilt (radians) +c* +c*PAR rggsm (I) : matrice de passage du repere geographique au +c*PAR : repere solaire magnetospherique +c*PAR rgsmg (I) : matrice de passage du repere solaire +c*PAR : magnetospherique au repere geographique +c* +c*PAR ds3 (I) : pas elementaire d'integration +c* +c*PAR xgc (I) : coordonnee geocentrique en x +c*PAR ygc (I) : coordonnee geocentrique en y +c*PAR zgc (I) : coordonnee geocentrique en z +c* +c*PAR rxp1 (O) : pas elementaire sur x (au point xgc) +c*PAR rxp2 (O) : pas elementaire sur y (au point ygc) +c*PAR rxp3 (O) : pas elementaire sur z (au point zgc) +c* +c*PAR ier (O) : code de retour +c* +c*NOT magin : 1 = dipole + g11 + h11 +c*NOT magin : 2 = champ DGRF 1945-1970 +c*NOT magin : 3 = champ DGRF 1970-1995 +c* +c*NOT magout : 0 = pas de champ externe +c*NOT magout : 1 = Tsyganenko 87 +c*NOT magout : 2 = Tsyganenko 89 +c*NOT magout : 3 = Kosik 97 +c* +c*NOT indgm : 1 ---> indice geomagnetique Kp +c*NOT indval : 1 ---> Kp = 0 , 0+ +c*NOT indval : 2 ---> Kp = 1- , 1 , 1+ +c*NOT indval : 3 ---> Kp = 2- , 2 , 2+ +c*NOT indval : 4 ---> Kp = 3- , 3 , 3+ +c*NOT indval : 5 ---> Kp = 4- , 4 , 4+ +c*NOT indval : 6 ---> Kp > 5- +c* +c*NOT indgm : 2 ---> indice geomagnetique Ae +c*NOT indval : 1 ---> Ae = 0 - 50 +c*NOT indval : 2 ---> Ae = 50 - 100 +c*NOT indval : 3 ---> Ae = 100 - 150 +c*NOT indval : 4 ---> Ae = 150 - 250 +c*NOT indval : 5 ---> Ae = 250 - 400 +c*NOT indval : 6 ---> Ae >= 400 +c* +c*NOT On utilise la formule generale dx/ds = bx/b +c*NOT ou ds est le pas elementaire le long de la ligne de champ, b le +c*NOT champ total, dx est la composante de ds suivant l'axe x et bx +c*NOT la composante du champ suivant ce meme axe. +c* +c*INF utilise : carsp, inmag, outma1, vspvcar +c* +c*HST version 1.0 - 99.03.31 - creation de la maglib au CDPP +c*HST version 2.0 - 01.05.30 - correction de commentaires de code +c*HST version 2.1 - 03.01.06 - corrections en compilation avec g77 +c*HST version 3.0 - 01.01.07 - mise a jour des modeles de champ interne +c*HST version 4.0 - 13.10.10 - mise a jour des modeles de champ interne +c*HST (DGRF05 et IGRF10) +c*HST version 5.0 - 17.02.23 - mise a jour des modeles de champ interne +c*HST (DGRF10 et IGRF15) +c*HST version 6.0 - 20.07.15 - mise a jour des modeles de champ interne +c*HST (DGRF15 et IGRF20) +c* +c*********************************************************************** +c* + implicit none +c +c --------------------------------- +c*FON Declaration identificateur rcs_id +c --------------------------------- +c + character rcs_id*100 +c +c -------------------------- +c*FON Declaration des parametres +c -------------------------- +c + integer indgm,indval,magin,magout + double precision year + double precision tilt + double precision rggsm(3,3), rgsmg(3,3) + double precision ds3 + double precision xgc, ygc, zgc + double precision rxp1, rxp2, rxp3 + integer ier +c +c --------------------------------- +c*FON Declaration des variables locales +c --------------------------------- +c + integer ier1,ier2,ier3,ier4 +c*LOC ier1,ier2,ier3,ier4 : codes retour des modules appeles +c + double precision rre,thetr,phir +c*LOC rre,thetr,phir : coordonnees spheriques +c + double precision bre,bte,bpe +c*LOC bre,bte,bpe : composantes radiale, tangentielle et azimuthale du champ +c*LOC : d'origine externe +c + double precision bri,bti,bpi +c*LOC bri,bti,bpi : composantes radiale, tangentielle et azimuthale du champ +c*LOC : d'origine interne +c + double precision br,bt,bp +c*LOC br,bt,bp : composantes radiale, tangentielle et azimuthale du champ total +c + double precision bi +c*LOC bi : module du champ +c + double precision bxt,byt,bzt +c*LOC bxt,byt,bzt : coordonnees solaires magnetospheriques x,y et z du champ +c + SAVE +c +c --------------------------------- +c*FON Affectation identificateur rcs_id +c --------------------------------- +c + data rcs_id /" + >$Id$"/ +c +c ****************** +c Debut de programme +c ****************** +c +c ----------------------------------- +c*FON Initialisations des codes de retour +c ----------------------------------- +c + ier = 0 + ier1 = 0 + ier2 = 0 + ier3 = 0 + ier4 = 0 +c +c --------------------------------------------------- +c*FON Initialisations des composantes du champ magnetique +c --------------------------------------------------- +c + bre = 0.d0 + bte = 0.d0 + bpe = 0.d0 +c + bri = 0.d0 + bti = 0.d0 + bpi = 0.d0 +c +c ---------------------------------------------------------------------- +c*FON Transformation des coordonnees geocentriques en coordonnees spheriques +c ---------------------------------------------------------------------- +c + call carsp (xgc,ygc,zgc,rre,thetr,phir,ier1) +c +c -------------------------------------------- +c*FON Calcul du champ magnetique d'origine interne +c -------------------------------------------- +c + call inmagv(magin,year,rre,thetr,phir,bri,bti,bpi,bi,ier2) +c +c -------------------------------------------- +c*FON Calcul du champ magnetique d'origine externe +c -------------------------------------------- +c + call outma1(magout,indgm,indval,tilt,rggsm,rgsmg, + > rre,thetr,phir,bre,bte,bpe,ier3) +c +c --------------------------------------------------------------- +c*FON Calcul des composantes du champ magnetique total en coordonnees +c*FON spheriques puis conversion en coordonnees cartesiennes. +c --------------------------------------------------------------- +c + br = bri + bre + bt = bti + bte + bp = bpi + bpe +c + call vspvcar(thetr,phir,br,bt,bp,bxt,byt,bzt,ier4) +c + bt = sqrt(bxt*bxt + byt*byt + bzt*bzt) +c +c ------------------------------------- +c*FON Calcul des pas elementaires sur x,y,z +c ------------------------------------- +c + rxp1 = ds3 * bxt / bt + rxp2 = ds3 * byt / bt + rxp3 = ds3 * bzt / bt +c +c **************** +c Fin de programme +c **************** +c + return + end diff --git a/maglib/src/getmag.f b/maglib/src/getmag.f index 2125c39..e95b8ab 100755 --- a/maglib/src/getmag.f +++ b/maglib/src/getmag.f @@ -100,21 +100,34 @@ c*FON by the geomagnetic calculations for the given year c -------------------------------------------------- c - call inigeom(iyear,imonth,iday,ihour,imin,isec,year,alfag, + if (iyear .ge. 2000) then + call inigeom(iyear,imonth,iday,ihour,imin,isec,year,alfag, > tetdip,phidip,alfas,deltas,rig,rgi,rgdip,rdipg, > rgsm,rsmg,tilt,rggsm,rgsmg,rgse,rseg,rigsm,rgsmi, > ifail) -c ------------------------------------------------------------- -c*FON Calculation of the position of the point in the magnetosphere -c ------------------------------------------------------------- -c -c call posmag(magout,isatex,year,rrmag,thetr,phir,alfag, > alfas,deltas,tilt,rgsm,rggsm,rgsmg,rgdip,rgse, > tetdip,phidip,xgsm,ygsm,zgsm,xgse,ygse,zgse,tgl, > flg,xlamb,tglc,hsl,clatgmr,clongmr,iposmg,ifail) + else + call inigeomv(iyear,imonth,iday,ihour,imin,isec,year,alfag, + > tetdip,phidip,alfas,deltas,rig,rgi,rgdip,rdipg, + > rgsm,rsmg,tilt,rggsm,rgsmg,rgse,rseg,rigsm,rgsmi, + > ifail) + + call posmagv(magout,isatex,year,rrmag,thetr,phir,alfag, + > alfas,deltas,tilt,rgsm,rggsm,rgsmg,rgdip,rgse, + > tetdip,phidip,xgsm,ygsm,zgsm,xgse,ygse,zgse,tgl, + > flg,xlamb,tglc,hsl,clatgmr,clongmr,iposmg,ifail) + + endif + +c ------------------------------------------------------------- +c*FON Calculation of the position of the point in the magnetosphere +c ------------------------------------------------------------- +c c return diff --git a/maglib/src/inmagv.f b/maglib/src/inmagv.f new file mode 100755 index 0000000..c2ee1de --- /dev/null +++ b/maglib/src/inmagv.f @@ -0,0 +1,144 @@ + subroutine inmagv (magin,year,rre,thet,phi,bri,bti,bpi,bi,ier) +c* +c*********************************************************************** +c* +c* "Copyright [c] CNES 98 - tous droits reserves" +c* ********************************************** +c* +c*PRO MAGLIB +c* +c*VER 99.03.31 - V 1.0 +c*VER 01.05.30 - V 2.0 +c*VER 03.01.06 - V 2.1 +c*VER 01.01.07 - V 3.0 +c*VER 13.10.10 - V 4.0 +c*VER 17.02.24 - V 5.0 +c*VER 20.07.15 - V 6.0 +c* +c*AUT spec. CNES - JC KOSIK - juin 1995 +c*AUT port. CISI +c*AUT adapt. AKKA +c* +c*ROL Theme : Calculs de geophysique +c*ROL Calcul du champ magnetique d'origine interne. +c* +c*PAR magin (I) : type de champ magnetique interne +c* +c*PAR year (I) : annee fractionnaire (>= 1945. pour magin = 2 ; +c*PAR : >= 1970. pour magin = 3) +c* +c*PAR rre (I) : distance radiale geocentrique (rayons terrestres) +c*PAR thet (I) : colatitude geocentrique (radians) +c*PAR phi (I) : longitude geocentrique (radians) +c* +c*PAR bri (O) : composante radiale du champ magnetique le long du +c* : meridien positive vers l'exterieur (gauss) +c*PAR bti (O) : composante tangentielle du champ magnetique le long +c*PAR : du meridien positive vers le sud (gauss) +c*PAR bp (O) : composante azimuthale du champ magnetique, positive +c*PAR : vers l'est (gauss) +c* +c*PAR bi (O) : module du champ (gauss) +c* +c*PAR ier (O) : code de retour +c* +c*NOT magin : 1 = dipole + g11 + h11 +c*NOT magin : 2 = champ DGRF 1945-1970 +c*NOT magin : 3 = champ DGRF 1970-1995 +c* +c*NOT ier : sans objet +c* +c*INF utilise : dipol, dgrf45_70, dgrf70_95 +c* +c*HST version 1.0 - 99.03.31 - creation de la maglib au CDPP +c*HST version 2.0 - 01.05.30 - correction de commentaires de code +c*HST version 2.1 - 03.01.06 - corrections en compilation avec g77 +c*HST version 3.0 - 01.01.07 - utilisation du champ 2005 +c*HST version 4.0 - 13.10.10 - Introduction du champ igrf10 et +c*HST passage de igrf05 a dgrf05 +c*HST version 5.0 - 17.02.24 - Introduction du champ igrf15 et +c*HST passage de igrf10 a dgrf10 +c*HST version 6.0 - 20.07.15 - Introduction du champ igrf20 et +c*HST passage de igrf15 a dgrf15 +c* +c*********************************************************************** +c* + implicit none +c +c --------------------------------- +c*FON Declaration identificateur rcs_id +c --------------------------------- +c + character rcs_id*100 +c +c -------------------------- +c*FON Declaration des parametres +c -------------------------- +c + integer magin + double precision year + double precision rre, thet, phi + double precision bri, bti, bpi, bi + integer ier +c +c --------------------------------- +c*FON Declaration des variables locales +c --------------------------------- +c + integer ier1,ier2,ier3 +c*LOC ier1,ier2,ier3 : codes retour des modules appeles +c + SAVE +c +c --------------------------------- +c*FON Affectation identificateur rcs_id +c --------------------------------- +c + data rcs_id /" + >$Id$"/ +c +c ****************** +c Debut de programme +c ****************** +c + ier = 0 + ier1 = 0 + ier2 = 0 + ier3 = 0 +c +c ------------------------------------ +c*FON Calcul du champ (dipole + g11 + h11) +c ------------------------------------ +c + if (magin .eq. 1) then +c + call dipol(year,rre,thet,phi,bri,bti,bpi,bi,ier1) +c +c ------------------------- +c*FON Calcul du champ DGRF 2000 +c ------------------------- +c + else if (magin .eq. 2) then +c + call dgrf45_70(year,rre,thet,phi,bri,bti,bpi,bi,ier2) +c +c ------------------------- +c*FON Calcul du champ DGRF 2005 +c ------------------------- +c + else if (magin .eq. 3) then +c + call dgrf70_95(year,rre,thet,phi,bri,bti,bpi,bi,ier2) +c + else if(magin .eq. 4) then +c + call dgrf95(year,rre,thet,phi,bri,bti,bpi,bi,ier2) +c + endif +c +c **************** +c Fin de programme +c **************** +c + return + end diff --git a/maglib/src/magtotv.f b/maglib/src/magtotv.f new file mode 100755 index 0000000..328b41c --- /dev/null +++ b/maglib/src/magtotv.f @@ -0,0 +1,199 @@ + subroutine magtotv (magin,year,magout,indgm,indval,tilt,rggsm, + > rgsmg,rre,thet,phi,br,bt,bp,bb,ier) +c* +c*********************************************************************** +c* +c* "Copyright [c] CNES 98 - tous droits reserves" +c* ********************************************** +c* +c*PRO MAGLIB +c* +c*VER 99.03.31 - V 1.0 +c*VER 01.05.30 - V 2.0 +c*VER 03.01.06 - V 2.1 +c*VER 01.01.07 - V 3.0 +c*VER 13.10.10 - V 4.0 +c*VER 17.02.24 - V 5.0 +c*VER 20.07.16 - V 6.0 +c* +c*AUT spec. CNES - JC KOSIK - juin 1995 +c*AUT port. CISI +c*AUT adapt. AKKA +c* +c*ROL Theme : Calculs de geophysique +c*ROL Calcul du champ magnetique total en un point. +c*ROL Ce champ est la combinaison possible d'un champ dipole +c*ROL ou complet et d'un champ externe, ici Tsyganenko. +c* +c*PAR magin (I) : type de champ magnetique interne +c* +c*PAR year (I) : annee fractionnaire (>= 1945.) +c* +c*PAR magout (I) : type de champ magnetique externe +c* +c*PAR indgm (I) : type d'indice geomagnetique Kp ou Ae +c*PAR indval (I) : indice geomagnetique : niveau d'amplitude du champ +c* +c*PAR tilt (I) : angle de tilt (radians) +c* +c*PAR rggsm (I) : matrice (3,3) de passage du repere geocentrique +c*PAR : au repere magnetospherique +c* +c*PAR rgsmg (I) : matrice (3,3) de passage du repere solaire +c*PAR : magnetospherique au repere geocentrique +c* +c*PAR rre (I) : distance radiale geocentrique (rayons terrestres) +c*PAR thet (I) : colatitude geocentrique (radians) +c*PAR phi (I) : longitude geocentrique (radians) +c* +c*PAR br (O) : composante radiale du champ total (gauss) +c*PAR bt (O) : composante tangentielle du champ total (gauss) +c*PAR bp (O) : composante azimuthale du champ total (gauss) +c* +c*PAR bb (O) : module du champ magnetique total (gauss) +c*PAR ier (O) : code de retour +c* +c*NOT magin : 1 = dipole + g11 + h11 +c*NOT magin : 2 = champ DGRF 1945-1970 +c*NOT magin : 3 = champ DGRF 1970-1995 +c* +c*NOT magout : 0 = pas de champ externe +c*NOT magout : 1 = Tsyganenko 87 +c*NOT magout : 2 = Tsyganenko 89 +c*NOT magout : 3 = Kosik 97 +c* +c*NOT indgm : 1 ---> indice geomagnetique Kp +c*NOT indval : 1 ---> Kp = 0 , 0+ +c*NOT indval : 2 ---> Kp = 1- , 1 , 1+ +c*NOT indval : 3 ---> Kp = 2- , 2 , 2+ +c*NOT indval : 4 ---> Kp = 3- , 3 , 3+ +c*NOT indval : 5 ---> Kp = 4- , 4 , 4+ +c*NOT indval : 6 ---> Kp > 5- +c* +c*NOT indgm : 2 ---> indice geomagnetique Ae +c*NOT indval : 1 ---> Ae = 0 - 50 +c*NOT indval : 2 ---> Ae = 50 - 100 +c*NOT indval : 3 ---> Ae = 100 - 150 +c*NOT indval : 4 ---> Ae = 150 - 250 +c*NOT indval : 5 ---> Ae = 250 - 400 +c*NOT indval : 6 ---> Ae >= 400 +c* +c*NOT ier : sans objet +c* +c*INF utilise : inmag, outma1 +c* +c*HST version 1.0 - 99.03.31 - creation de la maglib au CDPP +c*HST version 2.0 - 01.05.30 - Ajout indice Kosik 97 +c*HST version 2.1 - 03.01.06 - corrections en compilation avec g77 +c*HST version 3.0 - 01.01.07 - utilisation du champ 2005 +c*HST version 4.0 - 13.10.10 - Introductiondu champ IGRF10 et +c*HST passage de IGRF05 a DGRF05 +c*HST version 5.0 - 17.02.24 - Introduction du champ IGRF15 et +c*HST passage de IGRF10 a DGRF10 +c*HST version 5.0 - 20.07.16 - Introduction du champ IGRF20 et +c*HST passage de IGRF15 a DGRF15 +* +c*********************************************************************** +c* + implicit none +c +c --------------------------------- +c*FON Declaration identificateur rcs_id +c --------------------------------- +c + character rcs_id*100 +c +c -------------------------- +c*FON Declaration des parametres +c -------------------------- +c + integer magin + double precision year + integer magout, indgm, indval + double precision tilt + double precision rggsm(3,3), rgsmg(3,3) + double precision rre, thet, phi + double precision br, bt, bp, bb + integer ier +c +c --------------------------------- +c*FON Declaration des variables locales +c --------------------------------- +c + integer ier1,ier2 +c*LOC ier1,ier2 : codes retour des modules appeles +c + double precision bpe,bre,bte +c*LOC bpe,bre,bte : composantes radiale, tangentielle et azimuthale du champ +c* : externe +c + double precision bpi,bri,bti +c*LOC bpi,bri,bti : composantes radiale, tangentielle et azimuthale du champ +c* : interne +c + double precision bi +c*LOC bi : module du champ +c + SAVE +c +c --------------------------------- +c*FON Affectation identificateur rcs_id +c --------------------------------- +c + data rcs_id /" + >$Id$"/ +c +c ****************** +c Debut de programme +c ****************** +c + ier = 0 + ier1 = 0 + ier2 = 0 +c +c -------------------------------------------- +c*FON Initialisation a 0. des composantes du champ +c -------------------------------------------- +c + bre = 0.d0 + bte = 0.d0 + bpe = 0.d0 +c + bri = 0.d0 + bti = 0.d0 + bpi = 0.d0 +c + if (magin .ne. 0) then +c +c -------------------------------------------- +c*FON Calcul du champ magnetique d'origine interne +c -------------------------------------------- +c + call inmagv(magin,year,rre,thet,phi,bri,bti,bpi,bi,ier1) +c + endif +c +c -------------------------------------------- +c*FON Calcul du champ magnetique d'origine externe +c -------------------------------------------- +c + if (magout .ne. 0) then + call outma1(magout,indgm,indval,tilt,rggsm,rgsmg, + > rre,thet,phi,bre,bte,bpe,ier2) + endif +c +c -------------------------------- +c*FON Calcul du champ magnetique total +c -------------------------------- +c + br = bri + bre + bt = bti + bte + bp = bpi + bpe + bb = sqrt(br * br + bt * bt + bp * bp) +c +c **************** +c Fin de programme +c **************** +c + return + end diff --git a/maglib/src/posmag.f b/maglib/src/posmag.f index b47a830..16868ba 100755 --- a/maglib/src/posmag.f +++ b/maglib/src/posmag.f @@ -30,7 +30,7 @@ c*PAR magout (I) : type de champ magnetique externe c* c*PAR isatex (I) : type de satellite c* -c*PAR year (I) : annee fractionnaire >= 2005 +c*PAR year (I) : annee fractionnaire >= 2000 c* c*PAR rre (I) : distance radiale geocentrique (rayons terrestres) c*PAR thetr (I) : colatitude geocentrique (radians) @@ -325,10 +325,39 @@ c c c ------------------------------------------------------ c*FON Calcul du point conjugue situe dans le meme hemisphere -c*FON magin = 3 = champ IGRF 2005 +c*NOT magin : 2 = champ DGRF 2000 +c*NOT magin : 3 = champ DGRF 2005 +c*NOT magin : 4 = champ DGRF 2010 +c*NOT magin : 5 = champ DGRF 2015 +c*NOT magin : 6 = champ IGRF 2020 + c ------------------------------------------------------ c - magin = 3 + if (year .ge. 2000 .and. year .lt. 2005) then +c + magin = 2 +c + else if (year .ge. 2005 .and. year .lt. 2010) then +c + magin = 3 +c + else if (year .ge. 2010 .and. year .lt. 2015) then +c + magin = 4 +c + else if (year .ge. 2015 .and. year .lt. 2020) then +c + magin = 5 +c + else if (year .ge. 2020) then +c + magin = 6 +c + else + ier = 1 + endif +c +c indgm = 1 indval = 3 rb = 11.d0 diff --git a/maglib/src/posmagv.f b/maglib/src/posmagv.f new file mode 100755 index 0000000..f09c7b9 --- /dev/null +++ b/maglib/src/posmagv.f @@ -0,0 +1,481 @@ + subroutine posmagv (magout,isatex,year,rre,thetr,phir,alfag,alfas, + > deltas,tilt,rgsm,rggsm,rgsmg,rgdip,rgse, + > tetdip,phidip,xgsm,ygsm,zgsm,xgse,ygse,zgse, + > tgl,flg,xlamb,tglc,hsl,clatgm,clongm, + > iposmg,ier) +c* +c*********************************************************************** +c* +c* "Copyright [c] CNES 98 - tous droits reserves" +c* ********************************************** +c* +c*PRO MAGLIB +c* +c*VER 99.03.31 - V 1.0 +c*VER 01.05.30 - V 2.0 +c*VER 03.01.06 - V 2.1 +c*VER 01.01.07 - V 3.0 +c* +c*AUT spec. CNES - JC KOSIK - janvier 1991 +c*AUT port. CISI +c* +c*ROL Theme : Frontieres et regions +c*ROL Calcul de la position d'un satellite dans la magnetosphere +c*ROL par rapport aux differentres regions. +c*ROL Calcul aussi des parametres essentiels tgl, L et mlt, +c*ROL et des coordonnees dans divers systemes d axes. +c*ROL L'indicateur de positionnement iposmg est renseigne. +c* +c*PAR magout (I) : type de champ magnetique externe +c* +c*PAR isatex (I) : type de satellite +c* +c*PAR year (I) : annee fractionnaire >= 2005 +c* +c*PAR rre (I) : distance radiale geocentrique (rayons terrestres) +c*PAR thetr (I) : colatitude geocentrique (radians) +c*PAR phir (I) : longitude geocentrique (radians) +c* +c*PAR alfag (I) : ascension droite de Greenwich (radians) +c*PAR alfas (I) : ascension droite du soleil (radians) +c*PAR deltas (I) : declinaison du soleil (radians) +c* +c*PAR tilt (I) : angle de tilt (radians) +c* +c*PAR rgsm (I) : matrice de passage du repere geographique au +c*PAR : repere solaire magnetique +c* +c*PAR rggsm (I) : matrice de passage du repere geographique au +c*PAR : repere solaire magnetospherique +c*PAR rgsmg (I) : matrice de passage du repere magnetospherique +c*PAR : au repere geographique +c* +c*PAR rgdip (I) : matrice (3,3) de passage du repere geographique +c*PAR : au repere dipolaire +c* +c*PAR rgse (I) : matrice de passage du repere geographique +c*PAR : au repere solaire ecliptique +c* +c*PAR tetdip (I) : colatitude geocentrique du dipole (radians) +c*PAR phidip (I) : longitude geocentrique du dipole (radians) +c* +c*PAR xgsm (O) : coordonnee solaire magnetospherique en x (rayons terrestres) +c*PAR ygsm (O) : coordonnee solaire magnetospherique en y (rayons terrestres) +c*PAR zgsm (O) : coordonnee solaire magnetospherique en z (rayons terrestres) +c* +c*PAR xgse (O) : coordonnee solaire ecliptique en x (rayons terrestres) +c*PAR ygse (O) : coordonnee solaire ecliptique en y (rayons terrestres) +c*PAR zgse (O) : coordonnee solaire ecliptique en z (rayons terrestres) +c* +c*PAR tgl (O) : temps geomagnetique local du satellite +c* : (heures fractionnaires) +c*PAR flg (O) : parametre L de Galperin +c*PAR xlamb (O) : latitude invariante (radians) +c* +c*PAR tglc (O) : temps geomagnetique local du point conjugue +c* : (heures fractionnaires) +c* +c*PAR hsl (O) : hauteur de l'ombre (kilometres) +c* +c*PAR clatgm (O) : latitude geomagnetique (radians) +c*PAR clongm (O) : longitude geomagnetique (radians) +c*PAR iposmg (O) : tableau (dim=15) des indicateurs de positions +c* +c*PAR ier (O) : code de retour +c* +c*NOT magout : 0 = pas de champ externe +c*NOT magout : 1 = Tsyganenko 87 +c*NOT magout : 2 = Tsyganenko 89 +c*NOT magout : 3 = kosik 97 +c* +c*NOT isatex : 0 = satellite proche (apogee < 8. r. terre.) +c*NOT isatex : 1 = satellite excentrique (apogee > 8. r. terre.) +c* +c*NOT iposmg : 1 = le satellite appartient a la region +c*NOT iposmg : 0 = le satellite n'appartient pas a la region +c* +c*NOT iposmg(1) : calotte polaire +c*NOT iposmg(2) : ovale auroral +c*NOT iposmg(3) : cusp +c*NOT iposmg(4) : zone auroral diffuse +c*NOT iposmg(5) : ceinture de Van Allen pour xbelt = 6. +c*NOT iposmg(6) : plasmasphere +c*NOT iposmg(7) : couche neutre +c*NOT iposmg(8) : couche de plasma +c*NOT iposmg(9) : magnetosphere +c*NOT iposmg(10) : magnetopause +c*NOT iposmg(11) : magnetogaine +c*NOT iposmg(12) : onde de choc +c*NOT iposmg(13) : vent solaire +c*NOT iposmg(14) : ombre de la terre (satellite auroral) +c*NOT iposmg(15) : ceinture de Van Allen pour xbelt = 3.5 +c* +c*NOT ier : sans objet +c* +c*NOT common : util +c* +c*INF utilise : spcar, geogsm, tgml, geose, bwshff, mpsib, msheath +c*INF utilise : dconjr, invar, invlat, ggeom, oval, calpol, cusp +c*INF utilise : gussen, rbelt, chapel, cahsl, posnsh, pospsh +c* +c*HST version 1.0 - 99.03.31 - creation de la maglib au CDPP +c*HST version 2.0 - 01.05.30 - correction de commentaires de code +c*HST version 2.1 - 03.01.06 - corrections en compilation avec g77 +c*HST version 3.0 - 01.01.07 - utilisation du champ IGRF 2005 +c* +c*********************************************************************** +c* + implicit none +c +c --------------------------------- +c*FON Declaration identificateur rcs_id +c --------------------------------- +c + character rcs_id*100 +c +c -------------------------- +c*FON Declaration des parametres +c -------------------------- +c + integer magout, isatex + double precision year + double precision rre, thetr, phir + double precision alfag, alfas, deltas, tilt + double precision rgsm(3,3), rggsm(3,3), rgsmg(3,3) + double precision rgdip(3,3), rgse(3,3) + double precision tetdip, phidip + double precision tr(500), tthet(500), tphi(500) + double precision xgsm, ygsm, zgsm + double precision xgse, ygse, zgse + double precision tgl, flg, xlamb, tglc, hsl + double precision clatgm, clongm + integer iposmg(15) + double precision znsh + integer ier +c +c ---------------------------------- +c*FON Declaration des variables communes +c ---------------------------------- +c + double precision pi,dpi,rad,deg,pid,xmu,rayt +c +c*COM pi : constante pi (obtenue a partir de acos(-1.)) +c*COM dpi : constante 2 * pi +c*COM pid : constante pi / 2 +c*COM rad : facteur de conversion degres ----> radians +c*COM deg : facteur de conversion radians ----> degres +c*COM xmu : constante de gravitation terrestre (km**3/sec**2) +c*COM rayt : rayon equatorial terrestre (km) +c + common/util/pi,dpi,rad,deg,pid,xmu,rayt +c +c --------------------------------- +c*FON Declaration des variables locales +c --------------------------------- +c + integer imagn,ibwsh,imp,isheath,isolw,ioval,ipol + integer icusp,iguss,iallen,ichapp,ihsl,insh,ipsh +c*LOC Indicateurs d'appartenance aux zones +c + integer i +c*LOC i : indice de boucle +c + integer magin +c*LOC magin : type de champ magnetique interne (= 2 : DGRF45_70; 3 : DGRF70_95) +c + integer indgm +c*LOC indgm : type d'indice geomagnetique (=1 : Kp) +c + integer indval +c*LOC indval : indice geomagnetique : niveau d'amplitude du champ (1 a 6) +c + integer np +c*LOC np : nombre de points calcules +c + integer ier1,ier2,ier3,ier4,ier5,ier6,ier7,ier9 + integer ier10,ier11,ier12,ier13,ier14,ier15,ier16 + integer ier17,ier18,ier19 +c*LOC ier1 a ier19 : codes retour des modules appeles +c + double precision rmagc,thetc,phic +c*LOC rmagc,thetc,phic : coordonnees spheriques +c* + double precision xmag,ymag,zmag +c*LOC xmag,ymag,zmag : composantes cartesiennes en x, y et z +c + double precision xbelt +c*LOC xbelt : valeur de la limite externe de la ceinture +c + double precision rfin,rmax +c*LOC rfin,rmax : distances geocentriques d'arret et maximum autorisee +c + double precision dir +c*LOC dir : direction du trace des lignes de champ +c + double precision rb +c*LOC rb : distance subsolaire de la magnetopause +c + SAVE +c +c --------------------------------- +c*FON Affectation identificateur rcs_id +c --------------------------------- +c + data rcs_id /" + >$Id$"/ +c +c ****************** +c Debut de programme +c ****************** +c + ier = 0 + ier1 = 0 + ier2 = 0 + ier3 = 0 + ier4 = 0 + ier5 = 0 + ier6 = 0 + ier7 = 0 + ier9 = 0 + ier10 = 0 + ier11 = 0 + ier12 = 0 + ier13 = 0 + ier14 = 0 + ier15 = 0 + ier16 = 0 + ier17 = 0 + ier18 = 0 + ier19 = 0 +c +c -------------------------------------------------- +c*FON Mise a 0 du tableau d'indicateur de positionnement +c*FON Initialisations de valeurs par defaut +c -------------------------------------------------- +c + do 10 i = 1, 15 + iposmg(i) = 0 +10 continue +c +c dans tous les cas mais surtout pour satellites auroraux : + iposmg(9) = 1 +c + tgl = 999.0d0 + flg = 999.0d0 + xlamb = 999.0d0 * rad + tglc = 999.0d0 + clatgm = 999.0d0 * rad + clongm = 999.0d0 * rad +c +c ------------------------------------------------- +c*FON Calcul des coordonnees solaires magnetospheriques +c ------------------------------------------------- +c + call spcar(rre,thetr,phir,xmag,ymag,zmag,ier1) +c + call geogsm(rggsm,xmag,ymag,zmag,xgsm,ygsm,zgsm,ier2) +c +c --------------------------------------------------- +c*FON Calcul du temps geomagnetique local du satellite tgl +c --------------------------------------------------- +c + call tgml(rgsm,thetr,phir,tgl,ier3) +c +c ------------------------------------------ +c*FON Calcul des coordonnees solaire-ecliptiques +c ------------------------------------------ +c + call geose(rgse,xmag,ymag,zmag,xgse,ygse,zgse,ier4) +c +c -------------------------------------------------------- +c*FON Si nosat est superieur a 2, verification si le satellite +c*FON Appartient a la magnetosphere de sibeck +c -------------------------------------------------------- +c + imagn = 1 +c + if (isatex .eq. 1) then +c +c --------------------------------------------------------------- +c*FON Calcul de la position par rapport a l onde de choc de Fairfield +c*FON et la magnetopause de Sibeck +c --------------------------------------------------------------- +c + call bwshff(xgsm,ygsm,zgsm,ibwsh,ier5) +c + if (ibwsh .eq. 1) then + iposmg(12) = ibwsh + endif +c + call mpsib(xgsm,ygsm,zgsm,imp,ier6) +c + if (imp .eq. 1) then + iposmg(10) = imp + endif +c + call msheath (xgsm,ygsm,zgsm,imagn,isheath,isolw,ier7) +c + iposmg(13) = isolw + iposmg(11) = isheath + iposmg(9) = imagn +c + endif +c + if (imagn .eq. 1) then +c +c ------------------------------------------------------ +c*FON Calcul du point conjugue situe dans le meme hemisphere +c*NOT magin : 1 = dipole + g11 + h11 +c*NOT magin : 2 = champ DGRF 1945-1970 +c*NOT magin : 3 = champ DGRF 1970-1995 +c*NOT magin : 4 = champ DGRF 1995-2000 + +c ------------------------------------------------------ +c +c + if (year .ge. 1945 .and. year .lt. 1970) then +c + magin = 2 +c + else if (year .ge. 1970 .and. year .lt. 1995) then +c + magin = 3 +c + else if (year .ge. 1995 .and. year .lt. 2000) then +c + magin = 4 +c + else + ier = 1 + endif + + + indgm = 1 + indval = 3 + rb = 11.d0 + rfin = 1.d0 + rmax = 40.d0 + dir = -1.d0 +c + call dconjrv(magin,year,magout,indgm,indval,tilt,rb,rggsm, + > rgsmg,dir,rre,thetr,phir,rfin,rmax,np,tr, + > tthet,tphi,ier1) +c + if (ier1 .eq. 0) then +c + rmagc = tr(np) + thetc = tthet(np) + phic = tphi(np) +c + call invar(rmagc,thetc,phic,flg,ier2) +c + call invlat(flg,xlamb,ier3) +c +c --------------------------------------- +c*FON Calcul d appartenance a l ovale auroral +c --------------------------------------- +c +c ---------------------------------------------------------- +c*FON Calcul des coordonnees geomagnetiques et mlt (appele tglc) +c*FON du pied de la ligne de champ passant par le satellite +c ---------------------------------------------------------- +c + call ggeom(rgdip,thetc,phic,clatgm,clongm,ier9) + call tgml(rgsm,thetc,phic,tglc,ier10) + + call oval(tglc,clatgm,ioval,ier11) + call calpol(tglc,clatgm,ipol,ier12) +c +c ------------------------------------------- +c*FON Calcul de la position par rapport a la CUSP +c ------------------------------------------- +c + call cusp(clatgm,tglc,icusp,ier13) +c +c ------------------------------------------------------ +c*FON Calcul de la frontiere sud de la zone aurorale diffuse +c ------------------------------------------------------ +c + call gussen(tglc,xlamb,iguss,ier14) +c + iposmg(1) = ipol + iposmg(2) = ioval + iposmg(3) = icusp + iposmg(4) = iguss +c + endif +c +c ------------------------------------------------------- +c*FON Si le satellite est a moins de 7 rayons terrestres..... +c ------------------------------------------------------- +c + if (rre .lt. 7.d0) then +c +c ------------------------------------------------ +c*FON Calcul d appartenance a la ceinture de Van Allen +c ------------------------------------------------ +c +c Ceinture a 3.5 r. terr. +c ----------------------- +c + xbelt = 3.5d0 + call rbelt(xbelt,tetdip,phidip,rre,thetr, + + phir,iallen,ier15) + iposmg(15) = iallen +c +c Ceinture a 6.0 r. terr. +c ----------------------- +c + xbelt = 6.d0 + call rbelt(xbelt,tetdip,phidip,rre,thetr, + + phir,iallen,ier15) + iposmg(5) = iallen +c +c --------------------------------------------------- +c*FON Calcul d appartenance a la plasmasphere de Chappell +c --------------------------------------------------- +c + call chapel(tgl,flg,ichapp,ier16) + iposmg(6) = ichapp +c +c -------------------------------------------------------- +c*FON Calcul de la distance du satellite a l ombre de la terre +c -------------------------------------------------------- +c + call cahsl(alfag,alfas,deltas,rre,thetr,phir, + + hsl,ihsl,ier17) + iposmg(14) = ihsl +c + endif +c + endif +c +c ----------------------------------------------------------------- +c*FON Calcul de la distance a la couche neutre et a la couche de plasma +c*FON si satellite excentrique et xgsm inferieur a -9 rt +c ----------------------------------------------------------------- +c + if (isatex .eq. 1 .and. xgsm .le. -9.d0) then +c +c ---------------------------------------------------- +c*FON Calcul de la position par rapport a la couche neutre +c ---------------------------------------------------- +c + call posnsh(tilt,xgsm,ygsm,zgsm,insh,znsh,ier18) + iposmg(7) = insh +c +c ------------------------------------------------------- +c*FON Calcul de la position par rapport a la couche de plasma +c ------------------------------------------------------- +c + call pospsh(tilt,xgsm,ygsm,zgsm,ipsh,ier19) + iposmg(8) = ipsh +c + endif +c +c **************** +c Fin de programme +c **************** +c + return + end diff --git a/maglib/src/stepv.f b/maglib/src/stepv.f new file mode 100755 index 0000000..9cdcd3c --- /dev/null +++ b/maglib/src/stepv.f @@ -0,0 +1,236 @@ + subroutine stepv (magin,magout,indgm,indval,year,tilt, + > rggsm,rgsmg,x,y,z,ds,errin,ier) +c* +c*********************************************************************** +c* +c* "Copyright [c] CNES 98 - tous droits reserves" +c* ********************************************** +c* +c*PRO MAGLIB +c* +c*VER 99.03.31 - V 1.0 +c*VER 01.05.30 - V 2.0 +c*VER 03.01.06 - V 2.1 +c*VER 01.01.07 - V 3.0 +c* +c*AUT spec. CNES - JC KOSIK - avril 1996 +c*AUT port. CISI +c* +c*ROL Theme : Calculs de geophysique +c*ROL Calcul d'un pas elementaire de la ligne de champ par la +c*ROL methode d'integration de Runge Kutta a pas varaiable et +c*ROL les coefficients de Merson. +c* +c*PAR magin (I) : type de champ magnetique interne +c*PAR magout (I) : type de champ magnetique externe +c* +c*PAR indgm (I) : type d'indice geomagnetique Kp ou Ae +c*PAR indval (I) : indice geomagnetique : niveau d'amplitude du champ +c* +c*PAR year (I) : annee fractionnaire (>= 2000.) +c* +c*PAR tilt (I) : angle de tilt (radians) +c* +c*PAR rggsm (I) : matrice de passage du repere geographique au +c*PAR : repere solaire magnetospherique +c*PAR rgsmg (I) : matrice de passage du repere solaire +c*PAR : magnetospherique au repere geographique +c* +c*PAR x (I/O) : coordonnee geocentrique en x +c*PAR y (I/O) : coordonnee geocentrique en y +c*PAR z (I/O) : coordonnee geocentrique en z +c* +c*PAR ds (I) : pas elementaire d'integration +c* +c*PAR errin (I) : erreur maximale permise lors de l'integration +c* +c*PAR ier (O) : code retour +c* +c*NOT magin : 1 = dipole + g11 + h11 +c*NOT magin : 2 = champ DGRF 1945-1970 +c*NOT magin : 3 = champ DGRF 1970-1995 +c* +c*NOT magout : 0 = pas de champ externe +c*NOT magout : 1 = Tsyganenko 87 +c*NOT magout : 2 = Tsyganenko 89 +c*NOT magout : 3 = Kosik 97 +c* +c*NOT indgm : 1 ---> indice geomagnetique Kp +c*NOT indval : 1 ---> Kp = 0 , 0+ +c*NOT indval : 2 ---> Kp = 1- , 1 , 1+ +c*NOT indval : 3 ---> Kp = 2- , 2 , 2+ +c*NOT indval : 4 ---> Kp = 3- , 3 , 3+ +c*NOT indval : 5 ---> Kp = 4- , 4 , 4+ +c*NOT indval : 6 ---> Kp > 5- +c* +c*NOT indgm : 2 ---> indice geomagnetique Ae +c*NOT indval : 1 ---> Ae = 0 - 50 +c*NOT indval : 2 ---> Ae = 50 - 100 +c*NOT indval : 3 ---> Ae = 100 - 150 +c*NOT indval : 4 ---> Ae = 150 - 250 +c*NOT indval : 5 ---> Ae = 250 - 400 +c*NOT indval : 6 ---> Ae >= 400 +c* +c*NOT La division par 3 du pas, permet des coefficients plus simples. +c* +c*INF utilise : deriv +c* +c*HST version 1.0 - 99.03.31 - creation de la maglib au CDPP +c*HST version 2.0 - 01.05.30 - correction de commentaires de code +c*HST version 2.1 - 03.01.06 - corrections en compilation avec g77 +c*HST version 3.0 - 01.01.07 - extension au champ IGRF 2005 +c* +c*********************************************************************** +c* + implicit none +c +c --------------------------------- +c*FON Declaration identificateur rcs_id +c --------------------------------- +c + character rcs_id*100 +c +c -------------------------- +c*FON Declaration des parametres +c -------------------------- +c + integer magin, magout, indgm, indval + double precision year + double precision tilt + double precision rggsm(3,3), rgsmg(3,3) + double precision x, y, z + double precision ds + double precision errin + integer ier +c +c --------------------------------- +c*FON Declaration des variables locales +c --------------------------------- +c + double precision r11,r12,r13 +c*LOC r11,r12,r13 : pas elementaires sur x, y et z (aux points x, y et z) +c + double precision x2,y2,z2 +c*LOC x2,y2,z2 : coordonnees cartesiennes du point au 2eme pas +c + double precision r21,r22,r23 +c*LOC r21,r22,r23 : pas elementaires sur x, y et z (aux points x2, y2 et z2) +c + double precision x3,y3,z3 +c*LOC x3,y3,z3 : coordonnees cartesiennes du point au 3eme pas +c + double precision r31,r32,r33 +c*LOC r31,r32,r33 : pas elementaires sur x, y et z (aux points x3, y3 et z3) +c + double precision x4,y4,z4 +c*LOC x4,y4,z4 : coordonnees cartesiennes du point au 3eme pas +c + double precision r41,r42,r43 +c*LOC r41,r42,r43 : pas elementaires sur x, y et z (aux points x4, y4 et z4) +c + double precision x5,y5,z5 +c*LOC x5,y5,z5 : coordonnees cartesiennes du point au 5eme pas +c + double precision r51,r52,r53 +c*LOC r51,r52,r53 : pas elementaires sur x, y et z (aux points x5, y5 et z5) +c + double precision ds3,errcur +c*LOC Variables de travail intermediaires +c + SAVE +c +c --------------------------------- +c*FON Affectation identificateur rcs_id +c --------------------------------- +c + data rcs_id /" + >$Id$"/ +c +c ****************** +c Debut de programme +c ****************** +c + ier = 0 +c +10 continue +c + ds3 = ds / 3.d0 +c +c ---------------------------------------------------------- +c*FON 1er appel au sous programme de calcul des pas elementaires +c ---------------------------------------------------------- +c + call derivv(magin,magout,indgm,indval,year,tilt, + > rggsm,rgsmg,ds3,x,y,z,r11,r12,r13,ier) +c + x2 = x + r11 + y2 = y + r12 + z2 = z + r13 +c +c ------------------------------------------------------------ +c*FON 2ieme appel au sous programme de calcul des pas elementaires +c ------------------------------------------------------------ +c + call derivv(magin,magout,indgm,indval,year,tilt, + > rggsm,rgsmg,ds3,x2,y2,z2,r21,r22,r23,ier) +c + x3 = x + .5d0 * (r11 + r21) + y3 = y + .5d0 * (r12 + r22) + z3 = z + .5d0 * (r13 + r23) +c +c ------------------------------------------------------------ +c*FON 3ieme appel au sous programme de calcul des pas elementaires +c ------------------------------------------------------------ +c + call derivv(magin,magout,indgm,indval,year,tilt, + > rggsm,rgsmg,ds3,x3,y3,z3,r31,r32,r33,ier) +c + x4 = x +.375d0 * (r11 + 3.d0 * r31) + y4 = y +.375d0 * (r12 + 3.d0 * r32) + z4 = z +.375d0 * (r13 + 3.d0 * r33) +c +c ------------------------------------------------------------ +c*FON 4ieme appel au sous programme de calcul des pas elementaires +c ------------------------------------------------------------ +c + call derivv(magin,magout,indgm,indval,year,tilt, + > rggsm,rgsmg,ds3,x4,y4,z4,r41,r42,r43,ier) +c + x5 = x + 1.5d0 * (r11 -3.d0 * r31 + 4.d0 * r41) + y5 = y + 1.5d0 * (r12 -3.d0 * r32 + 4.d0 * r42) + z5 = z + 1.5d0 * (r13 -3.d0 * r33 + 4.d0 * r43) +c +c ------------------------------------------------------------ +c*FON 5ieme appel au sous programme de calcul des pas elementaires +c ------------------------------------------------------------ +c + call derivv(magin,magout,indgm,indval,year,tilt, + > rggsm,rgsmg,ds3,x5,y5,z5,r51,r52,r53,ier) +c + errcur = abs(r11 - 4.5d0 * r31 + 4.d0 * r41 -.5d0 * r51) + + > abs(r12 - 4.5d0 * r32 + 4.d0 * r42 -.5d0 * r52) + + > abs(r13 - 4.5d0 * r33 + 4.d0 * r43 -.5d0 * r53) +c +c ------------------------------- +c*FON Test sur la precision de calcul +c ------------------------------- +c + if (errcur .lt. errin) then + x = x + .5d0 * (r11 + 4.d0 * r41 + r51) + y = y + .5d0 * (r12 + 4.d0 * r42 + r52) + z = z + .5d0 * (r13 + 4.d0 * r43 + r53) + else + ds = ds * .5d0 + goto 10 + endif +c + if (errcur.lt. errin * .04d0 .and. abs(ds) .lt. 1.33d0) then + ds = ds * 1.5d0 + endif +c +c **************** +c Fin de programme +c **************** +c + return + end -- libgit2 0.21.2