Collimator.cpp 16.1 KB
Newer Older
gsell's avatar
gsell committed
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21
// ------------------------------------------------------------------------
// $RCSfile: Collimator.cpp,v $
// ------------------------------------------------------------------------
// $Revision: 1.1.1.1 $
// ------------------------------------------------------------------------
// Copyright: see Copyright.readme
// ------------------------------------------------------------------------
//
// Class: Collimator
//   Defines the abstract interface for a beam Collimator.
//
// ------------------------------------------------------------------------
// Class category: AbsBeamline
// ------------------------------------------------------------------------
//
// $Date: 2000/03/27 09:32:31 $
// $Author: fci $
//
// ------------------------------------------------------------------------

#include "AbsBeamline/Collimator.h"
kraus's avatar
kraus committed
22
#include "Physics/Physics.h"
23
#include "Algorithms/PartBunch.h"
gsell's avatar
gsell committed
24
#include "AbsBeamline/BeamlineVisitor.h"
25
#include "Fields/Fieldmap.h"
gsell's avatar
gsell committed
26
#include "Structure/LossDataSink.h"
kraus's avatar
kraus committed
27
#include "Utilities/Options.h"
kraus's avatar
kraus committed
28
#include "Solvers/SurfacePhysicsHandler.hh"
29
#include <memory>
gsell's avatar
gsell committed
30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55

extern Inform *gmsg;

using namespace std;

// Class Collimator
// ------------------------------------------------------------------------

Collimator::Collimator():
    Component(),
    filename_m(""),
    plane_m(OFF),
    position_m(0.0),
    PosX_m(0),
    PosY_m(0),
    PosZ_m(0),
    MomentumX_m(0),
    MomentumY_m(0),
    MomentumZ_m(0),
    time_m(0),
    id_m(0),
    informed_m(false),
    a_m(0.0),
    b_m(0.0),
    x0_m(0.0),
    y0_m(0.0),
56 57 58 59 60
    xstart_m(0.0),
    xend_m(0.0),
    ystart_m(0.0),
    yend_m(0.0),
    width_m(0.0),
gsell's avatar
gsell committed
61 62 63 64 65 66 67 68
    isAPepperPot_m(false),
    isASlit_m(false),
    isARColl_m(false),
    isACColl_m(false),
    isAWire_m(false),
    rHole_m(0.0),
    nHolesX_m(0),
    nHolesY_m(0),
69
    pitch_m(0.0),
70
    losses_m(0),
71 72
    lossDs_m(nullptr),
    sphys_m(NULL)
gsell's avatar
gsell committed
73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93
{}


Collimator::Collimator(const Collimator &right):
    Component(right),
    filename_m(right.filename_m),
    plane_m(right.plane_m),
    position_m(right.position_m),
    PosX_m(right.PosX_m),
    PosY_m(right.PosY_m),
    PosZ_m(right.PosZ_m),
    MomentumX_m(right.MomentumX_m),
    MomentumY_m(right.MomentumY_m),
    MomentumZ_m(right.MomentumZ_m),
    time_m(right.time_m),
    id_m(right.id_m),
    informed_m(right.informed_m),
    a_m(right.a_m),
    b_m(right.b_m),
    x0_m(right.x0_m),
    y0_m(right.y0_m),
94 95 96 97
    xstart_m(right.xstart_m),
    xend_m(right.xend_m),
    ystart_m(right.ystart_m),
    yend_m(right.yend_m),
98 99
    zstart_m(right.zstart_m),
    zend_m(right.zend_m),
100
    width_m(right.width_m),
gsell's avatar
gsell committed
101 102 103 104 105 106 107 108
    isAPepperPot_m(right.isAPepperPot_m),
    isASlit_m(right.isASlit_m),
    isARColl_m(right.isARColl_m),
    isACColl_m(right.isACColl_m),
    isAWire_m(right.isAWire_m),
    rHole_m(right.rHole_m),
    nHolesX_m(right.nHolesX_m),
    nHolesY_m(right.nHolesY_m),
109
    pitch_m(right.pitch_m),
110
    losses_m(0),
111 112
    lossDs_m(nullptr),
    sphys_m(NULL)
gsell's avatar
gsell committed
113
{
kraus's avatar
kraus committed
114
    setGeom();
gsell's avatar
gsell committed
115 116 117
}


118
Collimator::Collimator(const std::string &name):
gsell's avatar
gsell committed
119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135
    Component(name),
    filename_m(""),
    plane_m(OFF),
    position_m(0.0),
    PosX_m(0),
    PosY_m(0),
    PosZ_m(0),
    MomentumX_m(0),
    MomentumY_m(0),
    MomentumZ_m(0),
    time_m(0),
    id_m(0),
    informed_m(false),
    a_m(0.0),
    b_m(0.0),
    x0_m(0.0),
    y0_m(0.0),
136 137 138 139
    xstart_m(0.0),
    xend_m(0.0),
    ystart_m(0.0),
    yend_m(0.0),
140 141
    zstart_m(0.0),
    zend_m(0.0),
142
    width_m(0.0),
gsell's avatar
gsell committed
143 144 145 146 147 148 149 150
    isAPepperPot_m(false),
    isASlit_m(false),
    isARColl_m(false),
    isACColl_m(false),
    isAWire_m(false),
    rHole_m(0.0),
    nHolesX_m(0),
    nHolesY_m(0),
151
    pitch_m(0.0),
152
    losses_m(0),
153 154
    lossDs_m(nullptr),
    sphys_m(NULL)
gsell's avatar
gsell committed
155 156 157 158
{}


Collimator::~Collimator() {
kraus's avatar
kraus committed
159 160
    if(online_m)
        goOffline();
gsell's avatar
gsell committed
161 162 163 164 165 166 167
}


void Collimator::accept(BeamlineVisitor &visitor) const {
    visitor.visitCollimator(*this);
}

168 169 170

inline bool Collimator::isInColl(Vector_t R, Vector_t P, double recpgamma) {
    /**
kraus's avatar
kraus committed
171 172
       check if we are in the longitudinal
       range of the collimator
173
    */
kraus's avatar
kraus committed
174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210
    const double z = R(2) + P(2) * recpgamma;

    if((z > position_m) && (z <= position_m + getElementLength())) {
        if(isAPepperPot_m) {

            /**
               ------------
               |(0)|  |(0)|
               ----   -----
               |    a)    |
               |          |
               ----   -----
               |(0)|  |(0)|
               yL------------
               xL
               |---| d
               |--| pitch
               Observation: the area in a) is much larger than the
               area(s) (0). In a) particles are lost in (0)
               particles they are not lost.

            */
            const double h  =   pitch_m;
            const double xL = - 0.5 * h * (nHolesX_m - 1);
            const double yL = - 0.5 * h * (nHolesY_m - 1);
            bool alive = false;

            for(unsigned int m = 0; (m < nHolesX_m && (!alive)); m++) {
                for(unsigned int n = 0; (n < nHolesY_m && (!alive)); n++) {
                    double x_m = xL  + (m * h);
                    double y_m = yL  + (n * h);
                    /** are we in a) ? */
                    double rr = std::pow((R(0) - x_m) / rHole_m, 2) + std::pow((R(1) - y_m) / rHole_m, 2);
                    alive = (rr < 1.0);
                }
            }
            return !alive;
211 212
        } else if(isASlit_m || isARColl_m) {
            return (std::abs(R(0)) >= getXsize() || std::abs(R(1)) >= getYsize());
kraus's avatar
kraus committed
213 214 215 216
        } else if(isAWire_m) {
            ERRORMSG("Not yet implemented");
        } else {
            // case of an elliptic collimator
217
            return (std::pow(R(0) / getXsize(), 2.0) + std::pow(R(1) / getYsize(), 2.0)) >= 1.0;
kraus's avatar
kraus committed
218
        }
219
    }
kraus's avatar
kraus committed
220
    return false;
221 222
}

223
bool Collimator::apply(const size_t &i, const double &t, double E[], double B[]) {
gsell's avatar
gsell committed
224 225 226 227
    Vector_t Ev(0, 0, 0), Bv(0, 0, 0);
    return apply(i, t, Ev, Bv);
}

228
bool Collimator::apply(const size_t &i, const double &t, Vector_t &E, Vector_t &B) {
gsell's avatar
gsell committed
229 230 231 232
    const Vector_t &R = RefPartBunch_m->R[i] - Vector_t(dx_m, dy_m, ds_m); // including the missaligment
    const Vector_t &P = RefPartBunch_m->P[i];
    const double recpgamma = Physics::c * RefPartBunch_m->getdT() / sqrt(1.0  + dot(P, P));
    bool pdead = false;
233 234 235
    pdead = isInColl(R,P,recpgamma);

    if(pdead) {
kraus's avatar
kraus committed
236 237 238
        double frac = (R(2) - position_m) / P(2) * recpgamma;
        lossDs_m->addParticle(R,P, RefPartBunch_m->ID[i], t + frac * RefPartBunch_m->getdT(), 0);
        ++losses_m;
kraus's avatar
kraus committed
239
    }
gsell's avatar
gsell committed
240 241 242 243 244 245 246
    return pdead;
}

bool Collimator::apply(const Vector_t &R, const Vector_t &centroid, const double &t, Vector_t &E, Vector_t &B) {
    return false;
}

247 248
bool Collimator::checkCollimator(Vector_t r, Vector_t rmin, Vector_t rmax) {

kraus's avatar
kraus committed
249 250 251 252 253 254 255 256 257 258 259
    double r_start = sqrt(xstart_m * xstart_m + ystart_m * ystart_m);
    double r_end = sqrt(xend_m * xend_m + yend_m * yend_m);
    double r1 = sqrt(rmax(0) * rmax(0) + rmax(1) * rmax(1));
    bool isDead = false;
    if(rmax(2) >= zstart_m && rmin(2) <= zend_m) {
        if( r1 > r_start - 10.0 && r1 < r_end + 10.0 ){
            if(r(2) < zend_m && r(2) > zstart_m ) {
                int pflag = checkPoint(r(0), r(1));
                isDead = (pflag != 0);
            }
        }
260
    }
kraus's avatar
kraus committed
261
    return isDead;
262 263
}

gsell's avatar
gsell committed
264 265

// rectangle collimators in cyclotron cyclindral coordiantes
266
// without surfacephysics, the particle hitting collimator is deleted directly
gsell's avatar
gsell committed
267 268 269 270
bool Collimator::checkCollimator(PartBunch &bunch, const int turnnumber, const double t, const double tstep) {

    bool flagNeedUpdate = false;
    Vector_t rmin, rmax;
271

gsell's avatar
gsell committed
272
    bunch.get_bounds(rmin, rmax);
273 274
    double r_start = sqrt(xstart_m * xstart_m + ystart_m * ystart_m);
    double r_end = sqrt(xend_m * xend_m + yend_m * yend_m);
gsell's avatar
gsell committed
275
    double r1 = sqrt(rmax(0) * rmax(0) + rmax(1) * rmax(1));
276 277

    if(rmax(2) >= zstart_m && rmin(2) <= zend_m) {
kraus's avatar
kraus committed
278 279 280 281 282 283 284 285 286 287 288 289 290 291 292
        if( r1 > r_start - 10.0 && r1 < r_end + 10.0 ){
            size_t tempnum = bunch.getLocalNum();
            int pflag = 0;
            for(unsigned int i = 0; i < tempnum; ++i) {
                if(bunch.PType[i] == ParticleType::REGULAR && bunch.R[i](2) < zend_m && bunch.R[i](2) > zstart_m ) {
                    pflag = checkPoint(bunch.R[i](0), bunch.R[i](1));
                    if(pflag != 0) {
                        if (!sphys_m)
                            lossDs_m->addParticle(bunch.R[i], bunch.P[i], bunch.ID[i]);
                        bunch.Bin[i] = -1;
                        flagNeedUpdate = true;
                    }
                }
            }
        }
293 294 295
    }
    reduce(&flagNeedUpdate, &flagNeedUpdate + 1, &flagNeedUpdate, OpBitwiseOrAssign());
    if (flagNeedUpdate && sphys_m) {
kraus's avatar
kraus committed
296
        sphys_m->apply(bunch);
gsell's avatar
gsell committed
297
    }
298
    return flagNeedUpdate;
gsell's avatar
gsell committed
299 300 301 302 303 304
}

void Collimator::initialise(PartBunch *bunch, double &startField, double &endField, const double &scaleFactor) {
    RefPartBunch_m = bunch;
    position_m = startField;
    endField = position_m + getElementLength();
305

kraus's avatar
kraus committed
306 307
    sphys_m = getSurfacePhysics();

adelmann's avatar
adelmann committed
308
    if (!sphys_m) {
kraus's avatar
kraus committed
309 310 311 312
        if (filename_m == std::string(""))
            lossDs_m = std::unique_ptr<LossDataSink>(new LossDataSink(getName(), !Options::asciidump));
        else
            lossDs_m = std::unique_ptr<LossDataSink>(new LossDataSink(filename_m.substr(0, filename_m.rfind(".")), !Options::asciidump));
adelmann's avatar
adelmann committed
313
    }
kraus's avatar
kraus committed
314

kraus's avatar
kraus committed
315
    goOnline(-1e6);
gsell's avatar
gsell committed
316 317 318 319
}

void Collimator::initialise(PartBunch *bunch, const double &scaleFactor) {
    RefPartBunch_m = bunch;
320 321

    sphys_m = getSurfacePhysics();
adelmann's avatar
adelmann committed
322 323

    if (!sphys_m) {
kraus's avatar
kraus committed
324 325 326 327
        if (filename_m == std::string(""))
            lossDs_m = std::unique_ptr<LossDataSink>(new LossDataSink(getName(), !Options::asciidump));
        else
            lossDs_m = std::unique_ptr<LossDataSink>(new LossDataSink(filename_m.substr(0, filename_m.rfind(".")), !Options::asciidump));
adelmann's avatar
adelmann committed
328
    }
kraus's avatar
kraus committed
329

kraus's avatar
kraus committed
330
    goOnline(-1e6);
gsell's avatar
gsell committed
331 332 333 334
}

void Collimator::finalise()
{
kraus's avatar
kraus committed
335 336 337
    if(online_m)
        goOffline();
    *gmsg << "* Finalize probe" << endl;
gsell's avatar
gsell committed
338 339
}

kraus's avatar
kraus committed
340
void Collimator::goOnline(const double &) {
341
    print();
gsell's avatar
gsell committed
342 343 344 345 346 347 348 349 350 351 352 353

    PosX_m.reserve((int)(1.1 * RefPartBunch_m->getLocalNum()));
    PosY_m.reserve((int)(1.1 * RefPartBunch_m->getLocalNum()));
    PosZ_m.reserve((int)(1.1 * RefPartBunch_m->getLocalNum()));
    MomentumX_m.reserve((int)(1.1 * RefPartBunch_m->getLocalNum()));
    MomentumY_m.reserve((int)(1.1 * RefPartBunch_m->getLocalNum()));
    MomentumZ_m.reserve((int)(1.1 * RefPartBunch_m->getLocalNum()));
    time_m.reserve((int)(1.1 * RefPartBunch_m->getLocalNum()));
    id_m.reserve((int)(1.1 * RefPartBunch_m->getLocalNum()));
    online_m = true;
}

adelmann's avatar
adelmann committed
354 355 356
void Collimator::print() {
    if(RefPartBunch_m == NULL) {
        if(!informed_m) {
357
            std::string errormsg = Fieldmap::typeset_msg("BUNCH SIZE NOT SET", "warning");
358
            ERRORMSG(errormsg << endl);
adelmann's avatar
adelmann committed
359 360 361 362 363 364 365 366 367 368
            if(Ippl::myNode() == 0) {
                ofstream omsg("errormsg.txt", ios_base::app);
                omsg << errormsg << endl;
                omsg.close();
            }
            informed_m = true;
        }
        return;
    }

369
    *gmsg << level3;
adelmann's avatar
adelmann committed
370
    if(isAPepperPot_m)
371 372
        *gmsg << "* Pepperpot x= " << a_m << " y= " << b_m << " r= " << rHole_m
              << " nx= " << nHolesX_m << " ny= " << nHolesY_m << " pitch= " << pitch_m << endl;
adelmann's avatar
adelmann committed
373
    else if(isASlit_m)
374 375
        *gmsg << "* Slit x= " << getXsize() << " Slit y= " << getYsize()
              << " start= " << position_m << " fn= " << filename_m << endl;
adelmann's avatar
adelmann committed
376
    else if(isARColl_m)
377 378 379
        *gmsg << "* RCollimator a= " << getXsize() << " b= " << getYsize()
              << " start= " << position_m << " fn= " << filename_m
              << " ny= " << nHolesY_m << " pitch= " << pitch_m << endl;
adelmann's avatar
adelmann committed
380
    else if(isACColl_m)
381 382
        *gmsg << "* CCollimator angle start " << xstart_m << " (Deg) angle end " << ystart_m << " (Deg) "
              << "R start " << xend_m << " (mm) R rend " << yend_m << " (mm)" << endl;
adelmann's avatar
adelmann committed
383
    else if(isAWire_m)
384
        *gmsg << "* Wire x= " << x0_m << " y= " << y0_m << endl;
adelmann's avatar
adelmann committed
385
    else
386 387
        *gmsg << "* ECollimator a= " << getXsize() << " b= " << b_m << " start= " << position_m
              << " fn= " << filename_m << " ny= " << nHolesY_m << " pitch= " << pitch_m << endl;
adelmann's avatar
adelmann committed
388 389
}

gsell's avatar
gsell committed
390
void Collimator::goOffline() {
kraus's avatar
kraus committed
391 392 393
    if (online_m && lossDs_m)
        lossDs_m->save();
    online_m = false;
gsell's avatar
gsell committed
394 395 396 397 398 399
}

bool Collimator::bends() const {
    return false;
}

400
void Collimator::setOutputFN(std::string fn) {
gsell's avatar
gsell committed
401 402 403 404
    filename_m = fn;
}

string Collimator::getOutputFN() {
adelmann's avatar
adelmann committed
405 406
    if (filename_m == std::string(""))
        return getName();
kraus's avatar
kraus committed
407
    else
adelmann's avatar
adelmann committed
408
        return filename_m.substr(0, filename_m.rfind("."));
gsell's avatar
gsell committed
409 410
}

411 412 413 414
unsigned int Collimator::getLosses() const {
    return losses_m;
}

gsell's avatar
gsell committed
415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448
void Collimator::setXsize(double a) {
    a_m = a;
}

void Collimator::setYsize(double b) {
    b_m = b;
}

void Collimator::setXpos(double x0) {
    x0_m = x0;
}

void Collimator::setYpos(double y0) {
    y0_m = y0;
}


double Collimator::getXsize(double a) {
    return a_m;
}

double Collimator::getYsize(double b) {
    return b_m;
}

double Collimator::getXpos() {
    return x0_m;
}

double Collimator::getYpos() {
    return y0_m;

    // --------Cyclotron collimator
}
449 450
void Collimator::setXStart(double xstart) {
    xstart_m = xstart;
gsell's avatar
gsell committed
451 452
}

453 454
void Collimator::setXEnd(double xend) {
    xend_m = xend;
gsell's avatar
gsell committed
455 456
}

457 458
void Collimator::setYStart(double ystart) {
    ystart_m = ystart;
gsell's avatar
gsell committed
459 460
}

461 462
void Collimator::setYEnd(double yend) {
    yend_m = yend;
gsell's avatar
gsell committed
463 464
}

465 466 467 468 469 470 471 472
void Collimator::setZStart(double zstart) {
    zstart_m = zstart;
}

void Collimator::setZEnd(double zend) {
    zend_m = zend;
}

473 474
void Collimator::setWidth(double width) {
    width_m = width;
gsell's avatar
gsell committed
475 476
}

477 478
double Collimator::getXStart() {
    return xstart_m;
gsell's avatar
gsell committed
479 480
}

481 482
double Collimator::getXEnd() {
    return xend_m;
gsell's avatar
gsell committed
483 484
}

485 486
double Collimator::getYStart() {
    return ystart_m;
gsell's avatar
gsell committed
487 488
}

489 490
double Collimator::getYEnd() {
    return yend_m;
gsell's avatar
gsell committed
491 492
}

493 494 495 496 497 498 499 500
double Collimator::getZStart() {
    return zstart_m;
}

double Collimator::getZEnd() {
    return zend_m;
}

gsell's avatar
gsell committed
501
double Collimator::getWidth() {
502
    return width_m;
gsell's avatar
gsell committed
503 504 505 506 507 508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525 526 527 528 529 530 531 532 533 534 535 536 537 538 539 540 541 542 543 544 545
}

//-------------------------------

void Collimator::setRHole(double r) {
    rHole_m = r;
}
void Collimator::setNHoles(unsigned int nx, unsigned int ny) {
    nHolesX_m = nx;
    nHolesY_m = ny;
}
void Collimator::setPitch(double p) {
    pitch_m = p;
}


void Collimator::setPepperPot() {
    isAPepperPot_m = true;
}
void Collimator::setSlit() {
    isASlit_m = true;
}

void Collimator::setRColl() {
    isARColl_m = true;
}

void Collimator::setCColl() {
    isACColl_m = true;
}

void Collimator::setWire() {
    isAWire_m = true;
}
void Collimator::getDimensions(double &zBegin, double &zEnd) const {
    zBegin = position_m;
    zEnd = position_m + getElementLength();

    // zBegin = position_m - 0.005;
    //  zEnd = position_m + 0.005;

}

546 547
ElementBase::ElementType Collimator::getType() const {
    return COLLIMATOR;
gsell's avatar
gsell committed
548 549 550 551 552 553 554 555 556 557 558 559 560 561 562 563 564 565 566
}

string Collimator::getCollimatorShape() {
    if(isAPepperPot_m)
        return "PeperPot";
    else if(isASlit_m)
        return "Slit";
    else if(isARColl_m)
        return "RCollimator";
    else if(isACColl_m)
        return "CCollimator";
    else if(isAWire_m)
        return "Wire";
    else
        return "ECollimator";
}

void Collimator::setGeom() {

567 568
    double slope;
    if (xend_m == xstart_m)
kraus's avatar
kraus committed
569
        slope = 1.0e12;
570
    else
kraus's avatar
kraus committed
571
        slope = (yend_m - ystart_m) / (xend_m - xstart_m);
572 573 574 575 576 577

    double coeff2 = sqrt(1 + slope * slope);
    double coeff1 = slope / coeff2;
    double halfdist = width_m / 2.0;
    geom_m[0].x = xstart_m - halfdist * coeff1;
    geom_m[0].y = ystart_m + halfdist / coeff2;
gsell's avatar
gsell committed
578

579 580
    geom_m[1].x = xstart_m + halfdist * coeff1;
    geom_m[1].y = ystart_m - halfdist / coeff2;
gsell's avatar
gsell committed
581

582 583
    geom_m[2].x = xend_m + halfdist * coeff1;
    geom_m[2].y = yend_m - halfdist  / coeff2;
gsell's avatar
gsell committed
584

585 586
    geom_m[3].x = xend_m - halfdist * coeff1;
    geom_m[3].y = yend_m + halfdist / coeff2;
gsell's avatar
gsell committed
587 588 589 590

    geom_m[4].x = geom_m[0].x;
    geom_m[4].y = geom_m[0].y;

591
    if (zstart_m > zend_m){
kraus's avatar
kraus committed
592 593 594 595
        double tempz = 0.0;
        tempz = zstart_m;
        zstart_m = zend_m;
        zend_m = tempz;
596
    }
gsell's avatar
gsell committed
597 598 599 600 601 602 603 604 605 606 607 608 609 610 611 612
}


int Collimator::checkPoint(const double &x, const double &y) {
    int    cn = 0;

    for(int i = 0; i < 4; i++) {
        if(((geom_m[i].y <= y) && (geom_m[i+1].y > y))
           || ((geom_m[i].y > y) && (geom_m[i+1].y <= y))) {

            float vt = (float)(y - geom_m[i].y) / (geom_m[i+1].y - geom_m[i].y);
            if(x < geom_m[i].x + vt * (geom_m[i+1].x - geom_m[i].x))
                ++cn;
        }
    }
    return (cn & 1);  // 0 if even (out), and 1 if odd (in)
613
}