hext_kq.c

Go to the documentation of this file.

/*
    This file is part of magpar.

    Copyright (C) 2009 Josephat Kalezhi
    Copyright (C) 2002-2010 Werner Scholz

    www:   http://www.magpar.net/
    email: magpar(at)magpar.net

    magpar is free software; you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation; either version 2 of the License, or
    (at your option) any later version.

    magpar is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with magpar; if not, write to the Free Software
    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
*/

/*
Karlqvist type ring head model

Reference:

C. Denis Mee, Eric D. Daniel,
"Magnetic Recording Technology"
2nd ed. 1990
Sec. 2.1.2, p. 2.3
http://books.google.co.uk/books?id=extQfrcDveYC&pg=SA2-PA1&dq=recording+and+reproduce+theory+middleton&cd=1#v=onepage&q=recording%20and%20reproduce%20theory%20middleton&f=true

Hx = (Hg/pi)*(atan((g/2 + x)/y) + atan((g/2 - x)/y) )
Hy = -(Hg/2*pi)*log(((g/2 + x)2 + y2)/((g/2 -x)2 + y2))

Hx/Hy represent the horizontal/vertical component of the field.
Hg is the head gap field
g is the gapsize
x is the horizontal offset
y is the vertical offset
*/

static char const Id[] = "$Id: hext_kq.c 3024 2010-03-28 17:17:35Z scholz $\n\n";
static char const Td[] = "$Today: " __FILE__ "  " __DATE__ "  "  __TIME__  " $\n\n";

#include "griddata.h"

static int       doinit=1;
static int       fieldon=0;
static Vec       VHthis;
static PetscReal hextamp; /* Field strength in Tesla */

int Hext_kq_Init(GridData *gdata)
{
  MagparFunctionLogBegin;

  PetscTruth flg;
  ierr = PetscOptionsGetInt(PETSC_NULL,"-hext_kq",(PetscInt*)&fieldon,&flg);CHKERRQ(ierr);
  if (flg!=PETSC_TRUE) {
    fieldon=0;
    PetscPrintf(PETSC_COMM_WORLD,
    "Option -hext_kq not found, using default value: %i (field off)\n",
    fieldon
    );
  }

  if (!fieldon) {
    PetscPrintf(PETSC_COMM_WORLD,"Hext_kq off\n");
    MagparFunctionLogReturn(0);
  }
  PetscPrintf(PETSC_COMM_WORLD,"Hext_kq on\n");

  PetscReal hgapini;
  ierr = PetscOptionsGetReal(PETSC_NULL,"-hext_kq_gapini",&hgapini,&flg);CHKERRQ(ierr);
  if (flg!=PETSC_TRUE) {
    hgapini=50.0;
    PetscPrintf(PETSC_COMM_WORLD,
    "Option -hext_kq_gapini not found. Using default value %g\n",
    hgapini
    );
  }
  PetscPrintf(PETSC_COMM_WORLD,"hext_kq_gapini %g\n",hgapini);

  PetscReal gapsize;
  ierr = PetscOptionsGetReal(PETSC_NULL,"-hext_kq_gapsize",&gapsize,&flg);CHKERRQ(ierr);
  if (flg!=PETSC_TRUE) {
    gapsize=25.0;
    PetscPrintf(PETSC_COMM_WORLD,
    "Option -hext_kq_gapsize not found. Using default value %g\n",
    gapsize
    );
  }
  PetscPrintf(PETSC_COMM_WORLD,"kq_gap_size %g\n",gapsize);

  PetscReal x_offset;
  ierr = PetscOptionsGetReal(PETSC_NULL,"-hext_kq_x_offset",&x_offset,&flg);CHKERRQ(ierr);
  if (flg!=PETSC_TRUE) {
    x_offset=0.0;
    PetscPrintf(PETSC_COMM_WORLD,
    "Option -hext_kq_x_offset not found. Using default value %g\n",
    x_offset
    );
  }
  PetscPrintf(PETSC_COMM_WORLD,"kq_value of x_offset %g\n",x_offset);

  PetscReal z_offset;
  ierr = PetscOptionsGetReal(PETSC_NULL,"-hext_kq_z_offset",&z_offset,&flg);CHKERRQ(ierr);
  if (flg!=PETSC_TRUE) {
    z_offset=0.0;
    PetscPrintf(PETSC_COMM_WORLD,
    "Option -hext_kq_z_offset not found. Using default value %g\n",
    z_offset
    );
  }
  PetscPrintf(PETSC_COMM_WORLD,"kq_z_offset %g\n",z_offset);

  ierr = VecDuplicate(gdata->M,&VHthis);CHKERRQ(ierr);

  /* this sets the external field according to Karlqvist */

  /* make the vector of the external field accessible */
  PetscReal *ta_Hext;
  ierr = VecGetArray(VHthis,&ta_Hext);CHKERRQ(ierr);

  hgapini *=1000.0; /* in A/m */
  hgapini /= gdata->hscale/MU0; /* scaled */
  hextamp=hgapini;

  for (int j=0; j<gdata->ln_vert; j++) {
    PetscReal x,y,z;
    x = gdata->vertxyz[ND*gdata->vertl2g[j]+0];
    y = gdata->vertxyz[ND*gdata->vertl2g[j]+1];
    z = gdata->vertxyz[ND*gdata->vertl2g[j]+2];

    /* Using a ring head */

    PetscReal x_coord, z_coord;
    /* horizontal separation between head centre and point in question */
    x_coord = x_offset*1e-9 + x*gdata->lenscale;
    /* vertical separation between head centre and point in question */
    z_coord = z_offset*1e-9 + z*gdata->lenscale;

    PetscReal hx,hy,hz;
    hx = (hgapini/PETSC_PI)*
         (atan(((gapsize*1e-9)/2.0 + x_coord)/z_coord) +
          atan(((gapsize*1e-9)/2.0 - x_coord)/z_coord));
    hy = 0;
    hz = -(hgapini/(2*PETSC_PI))*
         log((((gapsize*1e-9)/2.0 + x_coord)*((gapsize*1e-9)/2.0 + x_coord) + z_coord*z_coord)/
             (((gapsize*1e-9)/2.0 - x_coord)*((gapsize*1e-9)/2.0 - x_coord) + z_coord*z_coord)
            );

    PetscPrintf(PETSC_COMM_WORLD,"Coordinates: (%g,%g,%g)\n",x,y,z);
    PetscPrintf(PETSC_COMM_WORLD,"Field:       (%g,%g,%g)\n",hx,hy,hz);

    ta_Hext[ND*j+0]= hx;
    ta_Hext[ND*j+1]= hy;
    ta_Hext[ND*j+2]= hz;
  }
  ierr = VecRestoreArray(VHthis,&ta_Hext);CHKERRQ(ierr);

  ierr = VecScale(gdata->VHext,hextamp);CHKERRQ(ierr);

  MagparFunctionLogReturn(0);
}


int Hext_kq(GridData *gdata,Vec VHextsum,PetscReal *hext)
{
  MagparFunctionInfoBegin;

  if (doinit) {
    ierr = Hext_kq_Init(gdata);CHKERRQ(ierr);
    doinit=0;
  }
  if (!fieldon) {
    MagparFunctionInfoReturn(0);
  }

  ierr = VecAXPY(VHextsum,1.0,VHthis);CHKERRQ(ierr);

  *hext=hextamp;

  MagparFunctionProfReturn(0);
}