/* rectmagnet.c - Rectangular magnet with fringe fields
 *
 * Author: Bas van der Geer, Pulsar Physics
 * Note  : Corrected version with continuously differenteial fields everywhere
 */

/* With special thanks to Bruno Muratory from Daresbury Laboratory
 * for providing us with the fringe field expressions
 */

#include <stdio.h>
#include <math.h>
#include "elem.h"

struct rectmagnet_info
{
  double a,b,Bo ;
  double dl,b1,b2 ;
} ;

static int rectmagnet_sim(gptpar *par,double t,struct rectmagnet_info *info) ;
static int rectmagnet_overlap(gptpar *par,double t,struct rectmagnet_info *info) ;
static int rectmagnet_simple(gptpar *par,double t,struct rectmagnet_info *info) ;

#define NGAP 10 /* Transient region is NGAP/b1 */

void rectmagnet_fix_init(gptinit *init)
{
  struct rectmagnet_info *info ;

  gptbuildECS( init ) ;

  if( gptgetargnum(init)!=6 )
    gpterror( "Syntax: %s(ECS,a,b,Bfield,dl,b1,b2)\n", gptgetname(init) ) ;

  info = (struct rectmagnet_info *)gptmalloc( sizeof(struct rectmagnet_info) ) ;

  info->a  = gptgetargdouble(init,1)/2 ;
  info->b  = gptgetargdouble(init,2)/2 ;
  info->Bo = gptgetargdouble(init,3) ;
  info->dl = gptgetargdouble(init,4) ;
  info->b1 = gptgetargdouble(init,5) ;
  info->b2 = gptgetargdouble(init,6) ;

  if( info->b1<0 )
    gpterror( "%s: b1 must be nonnegative", gptgetname(init) ) ;

  if( info->b1==0 )
    gptaddEBelement( init, rectmagnet_simple, gptfree, GPTELEM_LOCAL, info ) ;
  else if( (info->b+info->dl) * info->b1 < NGAP )
    gptaddEBelement( init, rectmagnet_overlap, gptfree, GPTELEM_GLOBAL, info ) ;
  else
    gptaddEBelement( init, rectmagnet_sim, gptfree, GPTELEM_GLOBAL, info ) ;
}

// Normal simulation function for 'top-hat' with smooth corners
static int rectmagnet_sim(gptpar *par,double t,struct rectmagnet_info *info)
{
  double d ;    	/* Distance from edge */
  double f, h, Bn ;

  if( fabs(X)>info->a )
    return 0 ;

  d = fabs(Z)-info->b-info->dl ;

  if( d*info->b1 > NGAP ) return 0 ;

  if( fabs(Y)*info->b1>=gpt_pi && d<=0)
  {
    gptremoveparticle(par) ;
    return 1 ;
  }

  f = info->b1*d + info->b2*(d*d-Y*Y) ;
  h = Y*(info->b1+2*info->b2*d) ;

  Bn = info->Bo/(1+2*exp(f)*cos(h)+exp(2*f)) ;
  BY = Bn*(1+exp(f)*cos(h)) ;
  BZ = -Bn*exp(f)*sin(h) ;

  if(Z<0) BZ=-BZ ;

  return 1 ;
}

// Overlapping fringes
static int rectmagnet_overlap(gptpar *par,double t,struct rectmagnet_info *info)
{
  if( fabs(X)>info->a )
    return 0 ;

  double dl = -Z-info->b-info->dl ;
  double dr =  Z-info->b-info->dl ;

  if( dl*info->b1>NGAP || dr*info->b1>NGAP )
    return 0 ;

  if( fabs(Y)*info->b1>=gpt_pi && dl<=0 && dr<=0)
  {
    gptremoveparticle(par) ;
    return 1 ;
  }

  double fl = info->b1*dl + info->b2*(dl*dl-Y*Y) ;
  double fr = info->b1*dr + info->b2*(dr*dr-Y*Y) ;
  double hl = Y*(info->b1+2*info->b2*dl) ;
  double hr = Y*(info->b1+2*info->b2*dr) ;

  double Bnl = info->Bo/(1+2*exp(fl)*cos(hl)+exp(2*fl)) ;
  double Bnr = info->Bo/(1+2*exp(fr)*cos(hr)+exp(2*fr)) ;

  BY = Bnl*(1+exp(fl)*cos(hl)) + Bnr*(1+exp(fr)*cos(hr)) - info->Bo ;
  BZ = Bnl*exp(fl)*sin(hl)     - Bnr*exp(fr)*sin(hr) ;

  return 1 ;
}

// No fringes
static int rectmagnet_simple(gptpar *par,double t,struct rectmagnet_info *info)
{
  if( fabs(X)>info->a || fabs(Z)>info->b )
    return 0 ;

  BY = info->Bo ;

  return 1 ;
}