/* rocot_check_mva.f -- translated by f2c (version 20100827).
   You must link the resulting object file with libf2c:
	on Microsoft Windows system, link with libf2c.lib;
	on Linux or Unix systems, link with .../path/to/libf2c.a -lm
	or, if you install libf2c.a in a standard place, with -lf2c -lm
	-- in that order, at the end of the command line, as in
		cc *.o -lf2c -lm
	Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,

		http://www.netlib.org/f2c/libf2c.zip
*/

#include "f2c.h"

/* Common Block Declarations */

struct {
    real grtube, pari, gjtube, tettu, phitu, gjxtu, gjytu, gjztu;
} tube_;

#define tube_1 tube_

struct {
    char modj[1];
} repaj_;

#define repaj_1 repaj_

/* Table of constant values */

static integer c__1 = 1;
static integer c__1000 = 1000;
static integer c__6 = 6;
static integer c__4 = 4;

/* Main program */ int MAIN__(void)
{
    /* Format strings */
    static char fmt_200[] = "(a)";
    static char fmt_300[] = "(i8)";
    static char fmt_100[] = "(3(1pe14.5))";

    /* System generated locals */
    olist o__1;
    cllist cl__1;

    /* Builtin functions */
    integer s_wsfe(cilist *), e_wsfe(void), do_fio(integer *, char *, ftnlen),
	     f_open(olist *), f_clos(cllist *);
    /* Subroutine */ int s_stop(char *, ftnlen);

    /* Local variables */
    static integer i__;
    static real bx[1000], by[1000], bz[1000], var[9]	/* was [3][3] */, 
	    vecp[9]	/* was [3][3] */, valp[3];
    static integer icrep;
    extern /* Subroutine */ int mat_cp_varmin__(integer *, real *, real *, 
	    real *, integer *, integer *, real *, real *, real *), 
	    writ_para_tube__(integer *), sigtest_(real *, real *, real *, 
	    integer *);

    /* Fortran I/O blocks */
    static cilist io___1 = { 0, 6, 0, fmt_200, 0 };
    static cilist io___2 = { 0, 6, 0, fmt_200, 0 };
    static cilist io___3 = { 0, 6, 0, fmt_200, 0 };
    static cilist io___4 = { 0, 6, 0, fmt_200, 0 };
    static cilist io___5 = { 0, 6, 0, fmt_200, 0 };
    static cilist io___9 = { 0, 1, 0, "(a)", 0 };
    static cilist io___10 = { 0, 1, 0, fmt_300, 0 };
    static cilist io___11 = { 0, 1, 0, fmt_100, 0 };
    static cilist io___17 = { 0, 1, 0, "(a)", 0 };
    static cilist io___18 = { 0, 1, 0, "(a)", 0 };
    static cilist io___19 = { 0, 1, 0, fmt_100, 0 };
    static cilist io___20 = { 0, 6, 0, fmt_200, 0 };
    static cilist io___21 = { 0, 6, 0, fmt_200, 0 };
    static cilist io___22 = { 0, 6, 0, fmt_200, 0 };
    static cilist io___23 = { 0, 6, 0, fmt_200, 0 };
    static cilist io___24 = { 0, 6, 0, fmt_200, 0 };
    static cilist io___25 = { 0, 6, 0, fmt_200, 0 };



/*     ***************************************************************0** */
/* *   Rocot_Check_MVA: program de test des modules varmin de le rocotlib */

/*     calcul la base (vecteurs propres) de la matrice de variance */
/*     d'un signal bx,by,bz */

/*     P. Robert, aout 2001 */
/*     ***************************************************************0** */



    s_wsfe(&io___1);
    e_wsfe();
    s_wsfe(&io___2);
    do_fio(&c__1, "=======================", (ftnlen)23);
    e_wsfe();
    s_wsfe(&io___3);
    do_fio(&c__1, "rocot_check_mva program", (ftnlen)23);
    e_wsfe();
    s_wsfe(&io___4);
    do_fio(&c__1, "=======================", (ftnlen)23);
    e_wsfe();
    s_wsfe(&io___5);
    e_wsfe();

/* *** chargement d'un signal test */

    sigtest_(bx, by, bz, &c__1000);

/* *** ouverture du fichier resu et chargement du signal d'entree */

    o__1.oerr = 0;
    o__1.ounit = 1;
    o__1.ofnmlen = 20;
    o__1.ofnm = "rocot_check_mva.resu";
    o__1.orl = 0;
    o__1.osta = 0;
    o__1.oacc = 0;
    o__1.ofm = 0;
    o__1.oblnk = 0;
    f_open(&o__1);

    writ_para_tube__(&c__1);

    s_wsfe(&io___9);
    do_fio(&c__1, "Input signal:", (ftnlen)13);
    e_wsfe();
    s_wsfe(&io___10);
    do_fio(&c__1, (char *)&c__1000, (ftnlen)sizeof(integer));
    e_wsfe();
    s_wsfe(&io___11);
    for (i__ = 1; i__ <= 1000; ++i__) {
	do_fio(&c__1, (char *)&bx[i__ - 1], (ftnlen)sizeof(real));
	do_fio(&c__1, (char *)&by[i__ - 1], (ftnlen)sizeof(real));
	do_fio(&c__1, (char *)&bz[i__ - 1], (ftnlen)sizeof(real));
    }
    e_wsfe();

/* *** calcul de variance du signal: */
/*     calcul de la matrice de variance, des  valeurs propres et */
/*     des vecteurs propres ; verification de l'orthogonalite */
/*     et du sens de la base trouvee, et passage des donnees */
/*     dans cette nouvelle base si icrep NE 0. */

    icrep = 1;

    mat_cp_varmin__(&c__6, bx, by, bz, &c__1000, &icrep, var, valp, vecp);

/* *** ecriture du signal dans la base de la variance */

    s_wsfe(&io___17);
    do_fio(&c__1, "-------------------------------------", (ftnlen)37);
    e_wsfe();
    s_wsfe(&io___18);
    do_fio(&c__1, "MVA signal:", (ftnlen)11);
    e_wsfe();
    s_wsfe(&io___19);
    for (i__ = 1; i__ <= 1000; ++i__) {
	do_fio(&c__1, (char *)&bx[i__ - 1], (ftnlen)sizeof(real));
	do_fio(&c__1, (char *)&by[i__ - 1], (ftnlen)sizeof(real));
	do_fio(&c__1, (char *)&bz[i__ - 1], (ftnlen)sizeof(real));
    }
    e_wsfe();

/* *** fermeture du fichier resultat */

    cl__1.cerr = 0;
    cl__1.cunit = 1;
    cl__1.csta = 0;
    f_clos(&cl__1);

    s_wsfe(&io___20);
    do_fio(&c__1, " ", (ftnlen)1);
    e_wsfe();
    s_wsfe(&io___21);
    do_fio(&c__1, "Input signal and MVA coord. syst. signal are in rocot_che"
	    "ck_mva.resu file", (ftnlen)73);
    e_wsfe();
    s_wsfe(&io___22);
    e_wsfe();
    s_wsfe(&io___23);
    do_fio(&c__1, "=======================", (ftnlen)23);
    e_wsfe();
    s_wsfe(&io___24);
    do_fio(&c__1, "end of rocot_check_mva ", (ftnlen)23);
    e_wsfe();
    s_wsfe(&io___25);
    do_fio(&c__1, "=======================", (ftnlen)23);
    e_wsfe();


    s_stop("rocot_check_mva.exe: normal termination", (ftnlen)39);
    return 0;
} /* MAIN__ */


/*     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX */

/* Subroutine */ int sigtest_(real *bx, real *by, real *bz, integer *nbp)
{
    /* Format strings */
    static char fmt_200[] = "(a,i9)";

    /* System generated locals */
    integer i__1;

    /* Builtin functions */
    integer s_wsfe(cilist *), e_wsfe(void), do_fio(integer *, char *, ftnlen);

    /* Local variables */
    static integer i__;
    static real dx;
    extern /* Subroutine */ int cal_b_tube__(real *, real *, real *, real *, 
	    real *, real *, real *);
    static real sx, sy, sz, bmoy;
    extern /* Subroutine */ int read_para_tube__(void), writ_para_tube__(
	    integer *);

    /* Fortran I/O blocks */
    static cilist io___27 = { 0, 6, 0, fmt_200, 0 };
    static cilist io___28 = { 0, 6, 0, fmt_200, 0 };
    static cilist io___29 = { 0, 6, 0, fmt_200, 0 };
    static cilist io___30 = { 0, 6, 0, fmt_200, 0 };
    static cilist io___31 = { 0, 6, 0, fmt_200, 0 };





    /* Parameter adjustments */
    --bz;
    --by;
    --bx;

    /* Function Body */
    read_para_tube__();

    dx = tube_1.grtube / 60.f;

    s_wsfe(&io___27);
    e_wsfe();
    s_wsfe(&io___28);
    do_fio(&c__1, "Computation of B vector along the trajectory", (ftnlen)44);
    e_wsfe();
    s_wsfe(&io___29);
    do_fio(&c__1, "--------------------------------------------", (ftnlen)44);
    e_wsfe();
    s_wsfe(&io___30);
    e_wsfe();

    s_wsfe(&io___31);
    do_fio(&c__1, "nb. of points=", (ftnlen)14);
    do_fio(&c__1, (char *)&(*nbp), (ftnlen)sizeof(integer));
    e_wsfe();

    writ_para_tube__(&c__6);

    i__1 = *nbp;
    for (i__ = 1; i__ <= i__1; ++i__) {

	sx = (real) (i__ - 1) * dx - tube_1.pari * 10.f;
	sy = -tube_1.pari;
	sz = 0.f;

	cal_b_tube__(&sx, &sy, &sz, &bx[i__], &by[i__], &bz[i__], &bmoy);

    }
    return 0;
} /* sigtest_ */


/*     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX */

/*     magsructlib: "magnetic structure libarary" */
/*                  bibliotheque de modules permettant le calcul */
/*                  du champ et de la densite de courant produit */
/*                  par une structure de courant telle un tube de courant */
/*                  a densite uniforme ou gaussienne. */

/*                  Biblioteque developpee pour les simulation CLUSTER */
/*                  Patrick ROBERT, CRPE, 1990-1995 */

/*     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX */

/* Subroutine */ int read_para_tube__(void)
{
    /* Format strings */
    static char fmt_200[] = "(a)";
    static char fmt_300[] = "(f10.3)";
    static char fmt_100[] = "(a1)";

    /* Builtin functions */
    integer s_wsfe(cilist *), do_fio(integer *, char *, ftnlen), e_wsfe(void),
	     s_rsle(cilist *), do_lio(integer *, integer *, char *, ftnlen), 
	    e_rsle(void), s_rsfe(cilist *), e_rsfe(void);

    /* Local variables */
    extern /* Subroutine */ int carsphe_(real *, real *, real *, real *, real 
	    *, real *);

    /* Fortran I/O blocks */
    static cilist io___37 = { 0, 6, 0, fmt_200, 0 };
    static cilist io___38 = { 0, 6, 0, fmt_200, 0 };
    static cilist io___39 = { 0, 6, 0, fmt_200, 0 };
    static cilist io___40 = { 0, 6, 0, fmt_200, 0 };
    static cilist io___41 = { 0, 6, 0, fmt_200, 0 };
    static cilist io___42 = { 0, 5, 0, 0, 0 };
    static cilist io___43 = { 0, 6, 0, fmt_300, 0 };
    static cilist io___44 = { 0, 6, 0, fmt_200, 0 };
    static cilist io___45 = { 0, 6, 0, fmt_200, 0 };
    static cilist io___46 = { 0, 5, 0, 0, 0 };
    static cilist io___47 = { 0, 6, 0, fmt_300, 0 };
    static cilist io___48 = { 0, 6, 0, fmt_200, 0 };
    static cilist io___49 = { 0, 6, 0, fmt_200, 0 };
    static cilist io___50 = { 0, 5, 0, 0, 0 };
    static cilist io___51 = { 0, 6, 0, fmt_300, 0 };
    static cilist io___52 = { 0, 6, 0, fmt_200, 0 };
    static cilist io___53 = { 0, 6, 0, fmt_200, 0 };
    static cilist io___54 = { 0, 6, 0, fmt_200, 0 };
    static cilist io___55 = { 0, 5, 0, 0, 0 };
    static cilist io___56 = { 0, 6, 0, fmt_300, 0 };
    static cilist io___57 = { 0, 6, 0, fmt_200, 0 };
    static cilist io___58 = { 0, 6, 0, fmt_200, 0 };
    static cilist io___59 = { 0, 5, 0, 0, 0 };
    static cilist io___60 = { 0, 6, 0, fmt_300, 0 };
    static cilist io___61 = { 0, 6, 0, fmt_200, 0 };
    static cilist io___62 = { 0, 6, 0, fmt_200, 0 };
    static cilist io___63 = { 0, 5, 0, fmt_100, 0 };
    static cilist io___64 = { 0, 6, 0, fmt_100, 0 };
    static cilist io___65 = { 0, 6, 0, fmt_200, 0 };





/*     *********************************************************** */
/*     lecture des parametres du tube */
/*     *********************************************************** */

    s_wsfe(&io___37);
    do_fio(&c__1, " ", (ftnlen)1);
    e_wsfe();
    s_wsfe(&io___38);
    do_fio(&c__1, "Current tube characteristics:", (ftnlen)29);
    e_wsfe();
    s_wsfe(&io___39);
    do_fio(&c__1, "-----------------------------", (ftnlen)29);
    e_wsfe();

    s_wsfe(&io___40);
    do_fio(&c__1, " ", (ftnlen)1);
    e_wsfe();
    s_wsfe(&io___41);
    do_fio(&c__1, "Radius of the current tube, in km", (ftnlen)33);
    e_wsfe();
    s_rsle(&io___42);
    do_lio(&c__4, &c__1, (char *)&tube_1.grtube, (ftnlen)sizeof(real));
    e_rsle();
    s_wsfe(&io___43);
    do_fio(&c__1, (char *)&tube_1.grtube, (ftnlen)sizeof(real));
    e_wsfe();

    s_wsfe(&io___44);
    do_fio(&c__1, " ", (ftnlen)1);
    e_wsfe();
    s_wsfe(&io___45);
    do_fio(&c__1, "impact parameter (min dist to the curr. tube), in km", (
	    ftnlen)52);
    e_wsfe();
    s_rsle(&io___46);
    do_lio(&c__4, &c__1, (char *)&tube_1.pari, (ftnlen)sizeof(real));
    e_rsle();
    s_wsfe(&io___47);
    do_fio(&c__1, (char *)&tube_1.pari, (ftnlen)sizeof(real));
    e_wsfe();

    s_wsfe(&io___48);
    do_fio(&c__1, " ", (ftnlen)1);
    e_wsfe();
    s_wsfe(&io___49);
    do_fio(&c__1, "Current density, in A/km2", (ftnlen)25);
    e_wsfe();
    s_rsle(&io___50);
    do_lio(&c__4, &c__1, (char *)&tube_1.gjtube, (ftnlen)sizeof(real));
    e_rsle();
    s_wsfe(&io___51);
    do_fio(&c__1, (char *)&tube_1.gjtube, (ftnlen)sizeof(real));
    e_wsfe();

    s_wsfe(&io___52);
    do_fio(&c__1, " ", (ftnlen)1);
    e_wsfe();
    s_wsfe(&io___53);
    do_fio(&c__1, "Theta angle of the current direction/ given frame (d)", (
	    ftnlen)53);
    e_wsfe();
    s_wsfe(&io___54);
    do_fio(&c__1, "teta ?", (ftnlen)6);
    e_wsfe();
    s_rsle(&io___55);
    do_lio(&c__4, &c__1, (char *)&tube_1.tettu, (ftnlen)sizeof(real));
    e_rsle();
    s_wsfe(&io___56);
    do_fio(&c__1, (char *)&tube_1.tettu, (ftnlen)sizeof(real));
    e_wsfe();

    s_wsfe(&io___57);
    do_fio(&c__1, " ", (ftnlen)1);
    e_wsfe();
    s_wsfe(&io___58);
    do_fio(&c__1, "Phi   angle of the current direction/ given frame (d)", (
	    ftnlen)53);
    e_wsfe();
    s_rsle(&io___59);
    do_lio(&c__4, &c__1, (char *)&tube_1.phitu, (ftnlen)sizeof(real));
    e_rsle();
    s_wsfe(&io___60);
    do_fio(&c__1, (char *)&tube_1.phitu, (ftnlen)sizeof(real));
    e_wsfe();

    carsphe_(&tube_1.gjtube, &tube_1.tettu, &tube_1.phitu, &tube_1.gjxtu, &
	    tube_1.gjytu, &tube_1.gjztu);


    s_wsfe(&io___61);
    do_fio(&c__1, " ", (ftnlen)1);
    e_wsfe();
    s_wsfe(&io___62);
    do_fio(&c__1, "Current density model (\"u\" or \"g\"/ uniform or gauss.)",
	     (ftnlen)53);
    e_wsfe();
    s_rsfe(&io___63);
    do_fio(&c__1, repaj_1.modj, (ftnlen)1);
    e_rsfe();
    s_wsfe(&io___64);
    do_fio(&c__1, repaj_1.modj, (ftnlen)1);
    e_wsfe();
    s_wsfe(&io___65);
    e_wsfe();


    return 0;
} /* read_para_tube__ */


/*     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX */

/* Subroutine */ int writ_para_tube__(integer *ifc)
{
    /* Format strings */
    static char fmt_100[] = "(\002 R=\002,f8.1,\002 km  IP=\002,f8.1,\002 km"
	    " J=\002,f6.3,\002 A/km2  tet=\002,f5.1,\002 d.  phi=\002,f5.1"
	    ",\002 d. mod=\002,a1)";

    /* Builtin functions */
    integer s_wsfe(cilist *), do_fio(integer *, char *, ftnlen), e_wsfe(void);

    /* Fortran I/O blocks */
    static cilist io___66 = { 0, 0, 0, fmt_100, 0 };





/*     *********************************************************** */
/*     ecriture des parametres du tube */
/*     *********************************************************** */

    io___66.ciunit = *ifc;
    s_wsfe(&io___66);
    do_fio(&c__1, (char *)&tube_1.grtube, (ftnlen)sizeof(real));
    do_fio(&c__1, (char *)&tube_1.pari, (ftnlen)sizeof(real));
    do_fio(&c__1, (char *)&tube_1.gjtube, (ftnlen)sizeof(real));
    do_fio(&c__1, (char *)&tube_1.tettu, (ftnlen)sizeof(real));
    do_fio(&c__1, (char *)&tube_1.phitu, (ftnlen)sizeof(real));
    do_fio(&c__1, repaj_1.modj, (ftnlen)1);
    e_wsfe();


    return 0;
} /* writ_para_tube__ */


/*     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX */

/* Subroutine */ int cal_b_tube__(real *xs, real *ys, real *zs, real *bxs, 
	real *bys, real *bzs, real *bmod)
{
    /* Builtin functions */
    double acos(doublereal);
    /* Subroutine */ int s_stop(char *, ftnlen);
    double exp(doublereal), sqrt(doublereal);

    /* Local variables */
    static real c__, x, y, z__, r2, pi, bx, by, bz, gr2, rsr2;
    extern /* Subroutine */ int xybvdh_(real *, real *, real *, real *, real *
	    , real *, real *, real *, real *), vdhxyb_(real *, real *, real *,
	     real *, real *, real *, real *, real *, real *);




/*     *********************************************************** */
/*     calcul de B dans le repere structure */
/*     repere structure=repere ou le tube est oblique */
/*     common /tube/=axe et module de J dans le repere structure */
/*     *********************************************************** */

    pi = acos(-1.f);
/*     valeur de Mu0 en nT*km/A */
/*     rmu=4.*pi*(1.e-7)*(1.e6) */
/*     c=Mu0*J/2 */
    c__ = pi * .2f * tube_1.gjtube;


/* *** passage au repere tube incline par rapport au repere structure */

    xybvdh_(xs, ys, zs, &tube_1.gjxtu, &tube_1.gjytu, &tube_1.gjztu, &x, &y, &
	    z__);

    r2 = x * x + y * y;
    gr2 = tube_1.grtube * tube_1.grtube;


    if (*(unsigned char *)repaj_1.modj != 'u' && *(unsigned char *)
	    repaj_1.modj != 'g') {
	s_stop("CAL_B_TUBE     *** ABORTED ! UNDEFINED MODEL OF J", (ftnlen)
		49);
    }
    if (*(unsigned char *)repaj_1.modj == 'u') {

/* *** densite uniforme */

	if (r2 >= gr2) {
	    bx = -c__ * gr2 * y / r2;
	    by = c__ * gr2 * x / r2;
	    bz = 0.f;

	} else {
	    bx = -c__ * y;
	    by = c__ * x;
	    bz = 0.f;
	}

    } else {

/* *** densite gaussienne */

	rsr2 = r2 / gr2;

	bx = -c__ * (1.f / rsr2) * (1.f - exp(-rsr2)) * y;
	by = c__ * (1.f / rsr2) * (1.f - exp(-rsr2)) * x;
	bz = 0.f;

/*     um=1.06 */
/*     rm=sqrt(um)*grtube */
/*     bmax=c*rm*(1.-exp(-um))/um */

    }

    *bmod = sqrt(bx * bx + by * by + bz * bz);

/* *** passage au repere structure */

    vdhxyb_(&bx, &by, &bz, &tube_1.gjxtu, &tube_1.gjytu, &tube_1.gjztu, bxs, 
	    bys, bzs);

    return 0;
} /* cal_b_tube__ */


/*     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX */

/* Subroutine */ int cal_j_tube__(real *xs, real *ys, real *zs, real *xsj, 
	real *ysj, real *zsj, real *rsj)
{
    /* System generated locals */
    real r__1;

    /* Builtin functions */
    double sqrt(doublereal);
    /* Subroutine */ int s_stop(char *, ftnlen);
    double exp(doublereal);

    /* Local variables */
    static real r__, x, y, z__, xj, yj, zj;
    extern /* Subroutine */ int xybvdh_(real *, real *, real *, real *, real *
	    , real *, real *, real *, real *), vdhxyb_(real *, real *, real *,
	     real *, real *, real *, real *, real *, real *);




/*     *********************************************************** */
/*     calcul de J dans le repere structure */
/*     repere structure=repere ou le tube est oblique */
/*     common /tube/=axe et module de J dans le repere structure */
/*     *********************************************************** */


/* *** passage au repere tube */

    xybvdh_(xs, ys, zs, &tube_1.gjxtu, &tube_1.gjytu, &tube_1.gjztu, &x, &y, &
	    z__);

    r__ = sqrt(x * x + y * y);

    if (*(unsigned char *)repaj_1.modj != 'u' && *(unsigned char *)
	    repaj_1.modj != 'g') {
	s_stop("CAL_B_TUBE     *** ABORTED ! UNDEFINED MODEL OF J", (ftnlen)
		49);
    }
    if (*(unsigned char *)repaj_1.modj == 'u') {

/* *** modele uniforme */

	if (r__ > tube_1.grtube) {
	    xj = 0.f;
	    yj = 0.f;
	    zj = 0.f;

	} else {
	    xj = 0.f;
	    yj = 0.f;
	    zj = tube_1.gjtube;
	}

    } else {

/* *** modele gaussien */

	xj = 0.f;
	yj = 0.f;
/* Computing 2nd power */
	r__1 = r__ / tube_1.grtube;
	zj = tube_1.gjtube * exp(-(r__1 * r__1));

    }
    *rsj = sqrt(xj * xj + yj * yj + zj * zj);

/* *** passage au repere structure */

    vdhxyb_(&xj, &yj, &zj, &tube_1.gjxtu, &tube_1.gjytu, &tube_1.gjztu, xsj, 
	    ysj, zsj);

    return 0;
} /* cal_j_tube__ */


/*     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX */

/* Subroutine */ int carsphe_(real *r__, real *tet, real *phi, real *x, real *
	y, real *z__)
{
    /* Builtin functions */
    double sin(doublereal), cos(doublereal);

    /* Local variables */
    static real pp;


/*     ****************************************************************** */
/*     auteur   :p. robert */
/*     categorie:astronomie et changement de base */
/*     objet    :conversion de coordonnees spheriques en cartesiennes */
/*     ****************************************************************** */

    pp = .017453292222222222f;

    *x = *r__ * sin(*tet * pp) * cos(*phi * pp);
    *y = *r__ * sin(*tet * pp) * sin(*phi * pp);
    *z__ = *r__ * cos(*tet * pp);

    return 0;
} /* carsphe_ */


/*     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX */

/* Subroutine */ int sincosp_(real *x, real *y, real *z__, real *st, real *ct,
	 real *sp, real *cp)
{
    /* Builtin functions */
    double sqrt(doublereal);
    /* Subroutine */ int s_stop(char *, ftnlen);

    /* Local variables */
    static real r__, rp;


/*     ****************************************************************** */
/*     auteur   :p. robert */
/*     categorie:astronomie et changement de base */
/*     objet    :sin et cos des angles polaires d"un vecteur cartesien */
/*     ****************************************************************** */

    r__ = sqrt(*x * *x + *y * *y + *z__ * *z__);

    *ct = *z__ / r__;
    if (dabs(*ct) > 1.f) {
	s_stop("SINCOSP        *** ABORTED ! ERROR 1", (ftnlen)36);
    }
    *st = sqrt(1.f - *ct * *ct);
    if (dabs(*st) > 1.f) {
	s_stop("SINCOSP        *** ABORTED ! ERROR 2", (ftnlen)36);
    }

    *sp = 0.f;
    *cp = 1.f;
    rp = sqrt(*x * *x + *y * *y);

    if (rp < 1e-30f) {
	return 0;
    }
    *sp = *y / rp;
    if (dabs(*sp) > 1.f) {
	s_stop("SINCOSP        *** ABORTED ! ERROR 3", (ftnlen)36);
    }
    *cp = *x / rp;
    if (dabs(*cp) > 1.f) {
	s_stop("SINCOSP        *** ABORTED ! ERROR 4", (ftnlen)36);
    }

    return 0;
} /* sincosp_ */


/*     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX */

/* Subroutine */ int vdhxyb_(real *x, real *y, real *b, real *bv, real *bd, 
	real *bh, real *v, real *d__, real *h__)
{
    /* Builtin functions */
    double sqrt(doublereal);

    /* Local variables */
    static real b1, b2, b3, cp, ct, sp, st, cpv, spv, stv;
    extern /* Subroutine */ int sincosp_(real *, real *, real *, real *, real 
	    *, real *, real *);


/*     ****************************************************************** */
/*     auteur   :p. robert */
/*     categorie:astronomie et changement de base */
/*     objet    :passage du repere ligne de force au repere vdh */
/*     ****************************************************************** */

/*      input=x,y,b avec axe x dans le plan du meridien magnetique */
/*      output=v,d,h */

    sincosp_(bv, bd, bh, &st, &ct, &sp, &cp);

    stv = sqrt(1.f - cp * cp * st * st);
    if (stv > 1e-30f) {
	spv = -sp / stv;
	cpv = cp * ct / stv;

    } else {
	spv = 0.f;
	cpv = 1.f;
    }

    b1 = cpv * *x - spv * *y;
    b2 = spv * *x + cpv * *y;
    b3 = *b;

    *v = cp * ct * b1 - sp * b2 + cp * st * b3;
    *d__ = sp * ct * b1 + cp * b2 + sp * st * b3;
    *h__ = -st * b1 + ct * b3;

    return 0;
} /* vdhxyb_ */


/*     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX */

/* Subroutine */ int xybvdh_(real *v, real *d__, real *h__, real *bv, real *
	bd, real *bh, real *x, real *y, real *b)
{
    /* Builtin functions */
    double sqrt(doublereal);

    /* Local variables */
    static real b1, b2, b3, cp, ct, sp, st, cpv, spv, stv;
    extern /* Subroutine */ int sincosp_(real *, real *, real *, real *, real 
	    *, real *, real *);


/*     ****************************************************************** */
/*     auteur   :p. robert */
/*     categorie:astronomie et changement de base */
/*     objet    :passage du repere vdh au repere ligne de force */
/*     ****************************************************************** */

/*      input=v,d,h */
/*      output=x,y,b */
/*             x est dans le plan du meridien magnetique */


    sincosp_(bv, bd, bh, &st, &ct, &sp, &cp);

    b1 = ct * cp * *v + ct * sp * *d__ - st * *h__;
    b2 = -sp * *v + cp * *d__;
    b3 = st * cp * *v + st * sp * *d__ + ct * *h__;

    stv = sqrt(1.f - cp * cp * st * st);
    if (stv > 1e-30f) {
	spv = -sp / stv;
	cpv = cp * ct / stv;

    } else {
	spv = 0.f;
	cpv = 1.f;
    }

    *x = cpv * b1 + spv * b2;
    *y = -spv * b1 + cpv * b2;
    *b = b3;

    return 0;
} /* xybvdh_ */

/* Main program alias */ int rocot_check_mva__ () { MAIN__ (); return 0; }
