[430] | 1 | !The Full Polymorphic Package |
---|
| 2 | !Copyright (C) Etienne Forest |
---|
| 3 | |
---|
| 4 | module definition |
---|
| 5 | ! use define_newda |
---|
| 6 | ! use precision_constants ! added to replace use define_newda |
---|
| 7 | ! use scratch_size |
---|
| 8 | ! use DABNEW |
---|
| 9 | ! use my_own_1D_TPSA |
---|
| 10 | use lielib_yang_berz, junk_no=>no,junk_nd=>nd,junk_nd2=>nd2,junk_ndpt=>ndpt,junk_nv=>nv |
---|
| 11 | ! use newda |
---|
| 12 | ! USE LIELIB_ETIENNE |
---|
| 13 | implicit none |
---|
| 14 | public |
---|
| 15 | logical(lp) :: newread=.false. ,newprint = .false. , first_time = .true. |
---|
| 16 | logical(lp) :: print77=.true. ,read77 = .true. |
---|
| 17 | logical(lp) :: no_ndum_check = .false. |
---|
| 18 | logical(lp),TARGET :: insane_PTC = .false. |
---|
| 19 | logical(lp),TARGET :: setknob = .false. !@1 Real part of knobs cannot set |
---|
| 20 | logical(lp),TARGET :: knob=.true. !@1 Knobs are effective |
---|
| 21 | integer, target :: npara_fpp !@1 position of last non-parameter tpsa variable |
---|
| 22 | ! complex(dp), parameter :: i_ = = cmplx(zero,one,kind=dp) |
---|
| 23 | ! complex(dp) :: i_ = (zero,one) ! cmplx(zero,one,kind=dp) |
---|
| 24 | complex(dp), parameter :: i_ = ( 0.0_dp,1.0_dp ) ! cmplx(zero,one,kind=dp) |
---|
| 25 | integer master |
---|
| 26 | ! integer,parameter::lnv=100 |
---|
| 27 | ! scratch variables |
---|
| 28 | INTEGER iassdoluser(ndumt) |
---|
| 29 | integer DUMMY,temp |
---|
| 30 | integer iass0user(ndumt) |
---|
| 31 | integer,parameter::ndim2=2*ndim |
---|
| 32 | integer,parameter::mmmmmm1=1,mmmmmm2=2,mmmmmm3=3,mmmmmm4=4 |
---|
| 33 | ! type (taylorlow) DUMMYl,templ !,DUMl(ndum) |
---|
| 34 | ! private NDC,NDC2,NDT,IREF,itu,iflow,jtune,nres,ifilt ,idpr |
---|
| 35 | ! private nplane,idsta,ista |
---|
| 36 | ! private xintex,dsta,sta,angle,rad,ps,rads,mx |
---|
| 37 | ! numerical differentiation by knobs |
---|
| 38 | logical(lp) :: knob_numerical=.false. |
---|
| 39 | real(dp) :: knob_eps(lnv)=1e-6_dp |
---|
| 40 | integer :: knob_i =0 |
---|
| 41 | INTEGER,PARAMETER::NMAX=20 |
---|
| 42 | integer,private,parameter::n_max=10 ! sagan stuff |
---|
| 43 | INTEGER, PARAMETER :: CASE1=1,CASE2=2, CASE0=0, CASEP1=-1,CASEP2=-2 |
---|
| 44 | INTEGER, PARAMETER :: CASET=3,CASETF1=4,CASETF2=5 |
---|
| 45 | INTEGER,PARAMETER :: ISPIN0R=1,ISPIN1R=3 |
---|
| 46 | logical(lp) :: doing_ac_modulation_in_ptc=.false. |
---|
| 47 | integer, target :: nb_ =0 ! global group index |
---|
| 48 | ! |
---|
| 49 | TYPE sub_taylor |
---|
| 50 | INTEGER j(lnv) |
---|
| 51 | INTEGER min,max |
---|
| 52 | END TYPE sub_taylor |
---|
| 53 | |
---|
| 54 | !!&1 |
---|
| 55 | TYPE taylor |
---|
| 56 | INTEGER I !@1 integer I is a pointer in old da-package of Berz |
---|
| 57 | ! type (taylorlow) j !@1 Taylorlow is an experimental type not supported |
---|
| 58 | END TYPE taylor |
---|
| 59 | !@2 UNIVERSAL_TAYLOR is used by Sagan in BMAD Code at Cornell |
---|
| 60 | !@2 Also used by MAD-XP |
---|
| 61 | TYPE UNIVERSAL_TAYLOR |
---|
| 62 | INTEGER, POINTER:: N,NV ! Number of coeeficients and number of variables |
---|
| 63 | REAL(DP), POINTER,dimension(:)::C ! Coefficients C(N) |
---|
| 64 | INTEGER, POINTER,dimension(:,:)::J ! Exponents of each coefficients J(N,NV) |
---|
| 65 | END TYPE UNIVERSAL_TAYLOR |
---|
| 66 | !@3 ---------------------------------------------</br> |
---|
| 67 | TYPE complextaylor |
---|
| 68 | type (taylor) r !@1 Real part |
---|
| 69 | type (taylor) i !@1 Imaginary part |
---|
| 70 | END TYPE complextaylor |
---|
| 71 | |
---|
| 72 | !&2 |
---|
| 73 | ! this is a real polymorphic type |
---|
| 74 | |
---|
| 75 | TYPE REAL_8 |
---|
| 76 | TYPE (TAYLOR) T !@1 USED IF TAYLOR |
---|
| 77 | REAL(DP) R !@1 USED IF REAL |
---|
| 78 | !&2 |
---|
| 79 | INTEGER KIND !@1 0,1,2,3 (1=REAL,2=TAYLOR,3=TAYLOR KNOB, 0=SPECIAL) |
---|
| 80 | INTEGER I !@1 USED FOR KNOBS AND SPECIAL KIND=0 |
---|
| 81 | REAL(DP) S !@1 SCALING FOR KNOBS AND SPECIAL KIND=0 |
---|
| 82 | LOGICAL(lp) :: ALLOC !@1 IF TAYLOR IS ALLOCATED IN DA-PACKAGE |
---|
| 83 | integer g,nb ! group index, number in group |
---|
| 84 | !&2 |
---|
| 85 | END TYPE REAL_8 |
---|
| 86 | !@3 ---------------------------------------------</br> |
---|
| 87 | TYPE double_complex |
---|
| 88 | type (complextaylor) t |
---|
| 89 | complex(dp) r |
---|
| 90 | logical(lp) alloc |
---|
| 91 | integer kind |
---|
| 92 | integer i,j |
---|
| 93 | complex(dp) s |
---|
| 94 | integer g,nb ! group index |
---|
| 95 | END TYPE double_complex |
---|
| 96 | |
---|
| 97 | type(taylor) varf1,varf2 |
---|
| 98 | type(complextaylor) varc1,varc2 |
---|
| 99 | |
---|
| 100 | !Radiation |
---|
| 101 | TYPE ENV_8 |
---|
| 102 | type (REAL_8) v |
---|
| 103 | type (REAL_8) e(ndim2) |
---|
| 104 | type (REAL_8) sigma0(ndim2) |
---|
| 105 | type (REAL_8) sigmaf(ndim2) |
---|
| 106 | END TYPE ENV_8 |
---|
| 107 | |
---|
| 108 | type spinor |
---|
| 109 | real(dp) x(3) |
---|
| 110 | ! real(dp) G |
---|
| 111 | end type spinor |
---|
| 112 | |
---|
| 113 | type spinor_8 |
---|
| 114 | type(real_8) x(3) |
---|
| 115 | end type spinor_8 |
---|
| 116 | |
---|
| 117 | type res_spinor_8 |
---|
| 118 | type(double_complex) x(3) |
---|
| 119 | end type res_spinor_8 |
---|
| 120 | |
---|
| 121 | ! scratch levels of DA using linked list |
---|
| 122 | |
---|
| 123 | type dascratch |
---|
| 124 | type(taylor), pointer :: t |
---|
| 125 | TYPE (dascratch),POINTER :: PREVIOUS |
---|
| 126 | TYPE (dascratch),POINTER :: NEXT |
---|
| 127 | end type dascratch |
---|
| 128 | |
---|
| 129 | TYPE dalevel |
---|
| 130 | INTEGER, POINTER :: N ! TOTAL ELEMENT IN THE CHAIN |
---|
| 131 | ! |
---|
| 132 | logical(lp),POINTER ::CLOSED |
---|
| 133 | TYPE (dascratch), POINTER :: PRESENT |
---|
| 134 | TYPE (dascratch), POINTER :: END |
---|
| 135 | TYPE (dascratch), POINTER :: START |
---|
| 136 | TYPE (dascratch), POINTER :: START_GROUND ! STORE THE GROUNDED VALUE OF START DURING CIRCULAR SCANNING |
---|
| 137 | TYPE (dascratch), POINTER :: END_GROUND ! STORE THE GROUNDED VALUE OF END DURING CIRCULAR SCANNING |
---|
| 138 | END TYPE dalevel |
---|
| 139 | |
---|
| 140 | !&1 |
---|
| 141 | TYPE DAMAP |
---|
| 142 | TYPE (TAYLOR) V(ndim2) ! Ndim2=6 but allocated to nd2=2,4,6 ! etienne_oct_2004 |
---|
| 143 | END TYPE DAMAP |
---|
| 144 | !&1 |
---|
| 145 | !@3 ---------------------------------------------</br> |
---|
| 146 | |
---|
| 147 | TYPE GMAP |
---|
| 148 | TYPE (TAYLOR) V(lnv) ! |
---|
| 149 | integer N |
---|
| 150 | END TYPE GMAP |
---|
| 151 | |
---|
| 152 | !&4 |
---|
| 153 | TYPE vecfield |
---|
| 154 | type (taylor) v(ndim2) !@1 <font face="Times New Roman">V<sub>i</sub>∂<sub>i</sub></font> Operator |
---|
| 155 | integer ifac !@1 Type of Factorization 0,1,-1 (One exponent, Dragt-Finn, Reversed Dragt-Finn) |
---|
| 156 | END TYPE vecfield |
---|
| 157 | !@3 ---------------------------------------------</br> |
---|
| 158 | TYPE pbfield |
---|
| 159 | type (taylor) h |
---|
| 160 | integer ifac |
---|
| 161 | integer nd_used |
---|
| 162 | END TYPE pbfield |
---|
| 163 | !&4 |
---|
| 164 | |
---|
| 165 | |
---|
| 166 | TYPE tree |
---|
| 167 | type (taylor) branch(ndim2) |
---|
| 168 | END TYPE tree |
---|
| 169 | |
---|
| 170 | !Radiation |
---|
| 171 | TYPE radtaylor |
---|
| 172 | type (taylor) v |
---|
| 173 | type (taylor) e(ndim2) |
---|
| 174 | END TYPE radtaylor |
---|
| 175 | |
---|
| 176 | !&5 |
---|
| 177 | TYPE DRAGTFINN |
---|
| 178 | real(dp) constant(ndim2) |
---|
| 179 | type (damap) Linear |
---|
| 180 | type (vecfield) nonlinear |
---|
| 181 | type (pbfield) pb |
---|
| 182 | END TYPE DRAGTFINN |
---|
| 183 | !@3 ---------------------------------------------</br> |
---|
| 184 | TYPE reversedragtfinn |
---|
| 185 | real(dp) CONSTANT(NDIM2) |
---|
| 186 | type (damap) Linear |
---|
| 187 | type (vecfield) nonlinear |
---|
| 188 | type (pbfield) pb |
---|
| 189 | END TYPE reversedragtfinn |
---|
| 190 | !@3 ---------------------------------------------</br> |
---|
| 191 | TYPE ONELIEEXPONENT |
---|
| 192 | real(dp) EPS |
---|
| 193 | type (vecfield) VECTOR |
---|
| 194 | type (pbfield) pb |
---|
| 195 | END TYPE ONELIEEXPONENT |
---|
| 196 | !@3 ---------------------------------------------</br> |
---|
| 197 | |
---|
| 198 | TYPE normalform |
---|
| 199 | type (damap) A_t ! Total A : A_t= A1 o A_rest |
---|
| 200 | type (damap) A1 !@1 Dispersion |
---|
| 201 | type (reversedragtfinn) A !@1 Linear A_t and nonlinear A_t |
---|
| 202 | type (dragtfinn) NORMAL !@1 Normal is the Normal Form R |
---|
| 203 | type (damap) DHDJ !@1 Contains the tunes in convenient form: extracted from NORMAL (=R) |
---|
| 204 | real(dp) TUNE(NDIM),DAMPING(NDIM) !@1 linear tune and linear damping |
---|
| 205 | integer nord,jtune !@1 nord=1 A1 first order in parameters |
---|
| 206 | integer NRES,M(NDIM,NRESO),PLANE(NDIM) !@1 NRES,M(NDIM,NRESO) -> resonances left in the map |
---|
| 207 | logical(lp) AUTO |
---|
| 208 | END TYPE normalform |
---|
| 209 | !@3 ---------------------------------------------</br> |
---|
| 210 | TYPE genfield |
---|
| 211 | type (taylor) h |
---|
| 212 | type (damap) m |
---|
| 213 | type (taylor) d(ndim,ndim) |
---|
| 214 | type (damap) linear |
---|
| 215 | type (damap) lineart |
---|
| 216 | type (damap) mt |
---|
| 217 | real(dp) constant(ndim2),eps |
---|
| 218 | integer imax !@1 imax=Maximum Number of Iteration (default=1000) |
---|
| 219 | integer ifac !@1 ifac = the map is raised to the power 1/ifac and iterated ifac times (default=1) |
---|
| 220 | logical(lp) linear_in !@1 Linear part is left in the map (default=.false.) |
---|
| 221 | integer no_cut !@1 Original map is not symplectic on and above no_cut |
---|
| 222 | END TYPE genfield |
---|
| 223 | |
---|
| 224 | !@3 ---------------------------------------------</br> |
---|
| 225 | |
---|
| 226 | |
---|
| 227 | TYPE pbresonance |
---|
| 228 | type (pbfield) cos,sin |
---|
| 229 | integer ifac |
---|
| 230 | END TYPE pbresonance |
---|
| 231 | !@3 ---------------------------------------------</br> |
---|
| 232 | TYPE vecresonance |
---|
| 233 | type (vecfield) cos,sin |
---|
| 234 | integer ifac |
---|
| 235 | END TYPE vecresonance |
---|
| 236 | !@3 ---------------------------------------------</br> |
---|
| 237 | TYPE taylorresonance |
---|
| 238 | type (taylor) cos,sin |
---|
| 239 | END TYPE taylorresonance |
---|
| 240 | |
---|
| 241 | TYPE beamenvelope !@2 A kind of Normal Form for Radiative Envelope |
---|
| 242 | ! radiation normalization |
---|
| 243 | type (damap) transpose ! Transpose of map which acts on polynomials |
---|
| 244 | type (taylor) bij ! Represents the stochastic kick at the end of the turn Env_f=M Env_f M^t + B |
---|
| 245 | TYPE (pbresonance) bijnr ! Equilibrium beam sizes in resonance basis |
---|
| 246 | real(dp) s_ij0(6,6) ! equilibrium beam sizes |
---|
| 247 | type (taylor) sij0 ! equilibrium beam sizes |
---|
| 248 | real(dp) emittance(3),tune(3),damping(3) |
---|
| 249 | logical(lp) AUTO,STOCHASTIC |
---|
| 250 | real(dp) KICK(3) ! fake kicks for tracking stochastically |
---|
| 251 | type (damap) STOCH ! Diagonalized of stochastic part of map for tracking stochastically |
---|
| 252 | END TYPE beamenvelope |
---|
| 253 | |
---|
| 254 | |
---|
| 255 | type tree_element !@1 USED FOR FAST TRACKING IN O_TREE_ELEMENT.F90 |
---|
| 256 | real(dp) , DIMENSION(:), POINTER :: CC |
---|
| 257 | real(dp) , DIMENSION(:), POINTER :: fix |
---|
| 258 | integer, DIMENSION(:), POINTER :: JL,JV |
---|
| 259 | INTEGER,POINTER :: N,ND2,no |
---|
| 260 | end type tree_element |
---|
| 261 | |
---|
| 262 | type spinmatrix |
---|
| 263 | type(real_8) s(3,3) |
---|
| 264 | end type spinmatrix |
---|
| 265 | |
---|
| 266 | type damapspin |
---|
| 267 | type(damap) M |
---|
| 268 | type(spinmatrix) s |
---|
| 269 | ! type(real_8) s(3,3) |
---|
| 270 | real(dp) e_ij(6,6) |
---|
| 271 | end type damapspin |
---|
| 272 | |
---|
| 273 | type normal_spin |
---|
| 274 | type(normalform) N ! regular orbital normal form |
---|
| 275 | type(damapspin) a1 ! brings to fixed point |
---|
| 276 | type(damapspin) ar ! normalises around the fixed point |
---|
| 277 | type(damapspin) as ! pure spin map |
---|
| 278 | type(damapspin) a_t ! !! (a_t%m,a_t%s) = (a1%m, I ) o (I ,as%s) o (ar%m,I) |
---|
| 279 | !!! extra spin info |
---|
| 280 | integer M(NDIM,NRESO),MS(NRESO),NRES ! orbital and spin resonances to be left in the map |
---|
| 281 | type(real_8) n0(3) ! n0 vector |
---|
| 282 | type(real_8) theta0 ! angle for the matrix around the orbit (analogous to linear tunes) |
---|
| 283 | real(dp) nu ! spin tune |
---|
| 284 | !!!Envelope radiation stuff |
---|
| 285 | real(dp) s_ij0(6,6) ! equilibrium beam sizes |
---|
| 286 | complex(dp) s_ijr(6,6) ! equilibrium beam sizes in resonance basis |
---|
| 287 | real(dp) emittance(3),tune(3),damping(3) ! equilibrium emittances (partially well defined only for infinitesimal damping) |
---|
| 288 | logical(lp) AUTO,STOCHASTIC |
---|
| 289 | real(dp) KICK(3) ! fake kicks for tracking stochastically |
---|
| 290 | real(dp) STOCH(6,6) ! Diagonalized of stochastic part of map for tracking stochastically |
---|
| 291 | real(dp) STOCH_inv(6,6) ! Diagonalized of stochastic part of map for tracking stochastically |
---|
| 292 | end type normal_spin |
---|
| 293 | |
---|
| 294 | include "a_def_frame_patch_chart.inc" |
---|
| 295 | include "a_def_all_kind.inc" |
---|
| 296 | include "a_def_sagan.inc" |
---|
| 297 | include "a_def_element_fibre_layout.inc" |
---|
| 298 | |
---|
| 299 | type(fibre), pointer :: lost_fibre |
---|
| 300 | type(integration_node), pointer :: lost_node |
---|
| 301 | |
---|
| 302 | type rf_phasor |
---|
| 303 | real(dp) x(2) |
---|
| 304 | real(dp) om |
---|
| 305 | real(dp) t |
---|
| 306 | end type rf_phasor |
---|
| 307 | |
---|
| 308 | type rf_phasor_8 |
---|
| 309 | type(real_8) x(2) |
---|
| 310 | type(real_8) om |
---|
| 311 | reaL(DP) t |
---|
| 312 | end type rf_phasor_8 |
---|
| 313 | |
---|
| 314 | type probe |
---|
| 315 | real(dp) x(6) |
---|
| 316 | type(spinor) s(ISPIN0R:ISPIN1R) |
---|
| 317 | type(rf_phasor) AC |
---|
| 318 | logical u |
---|
| 319 | type(integration_node),pointer :: lost_node |
---|
| 320 | end type probe |
---|
| 321 | |
---|
| 322 | type probe_8 |
---|
| 323 | type(real_8) x(6) ! Polymorphic orbital ray |
---|
| 324 | type(spinor_8) s(ISPIN0R:ISPIN1R) ! Polymorphic spin |
---|
| 325 | type(rf_phasor_8) AC ! Modulation |
---|
| 326 | real(dp) E_ij(6,6) ! Envelope |
---|
| 327 | ! stuff for exception |
---|
| 328 | logical u |
---|
| 329 | type(integration_node),pointer :: lost_node |
---|
| 330 | end type probe_8 |
---|
| 331 | |
---|
| 332 | type TEMPORAL_PROBE |
---|
| 333 | TYPE(probe) XS |
---|
| 334 | TYPE(INTEGRATION_NODE), POINTER :: NODE |
---|
| 335 | real(DP) DS,POS(6) |
---|
| 336 | END type TEMPORAL_PROBE |
---|
| 337 | |
---|
| 338 | type TEMPORAL_BEAM |
---|
| 339 | TYPE(TEMPORAL_PROBE), pointer :: TP(:) |
---|
| 340 | real(DP) a(3),ent(3,3),p0c,total_time |
---|
| 341 | integer n |
---|
| 342 | type(integration_node),pointer :: c ! pointer close to a(3) |
---|
| 343 | type(internal_state) state |
---|
| 344 | END type TEMPORAL_BEAM |
---|
| 345 | |
---|
| 346 | contains |
---|
| 347 | |
---|
| 348 | SUBROUTINE RESET_APERTURE_FLAG(complete) |
---|
| 349 | IMPLICIT NONE |
---|
| 350 | logical(lp), optional :: complete |
---|
| 351 | logical(lp) comp |
---|
| 352 | ! IF(c_%WATCH_USER) THEN |
---|
| 353 | ! IF(.NOT.ALLOW_TRACKING) THEN |
---|
| 354 | ! WRITE(6,*) " EXECUTION OF THE CODE MUST BE INTERRUPTED AT YOUR REQUEST" |
---|
| 355 | ! WRITE(6,*) " YOU DID NOT CHECK THE APERTURE STATUS" |
---|
| 356 | ! WRITE(6,*) " USING A CALL TO PRODUCE_APERTURE_FLAG" |
---|
| 357 | ! WRITE(6,*) " BEFORE CALLING A TRACKING FUNCTION" |
---|
| 358 | ! STOP 666 |
---|
| 359 | ! ENDIF |
---|
| 360 | ! IF(.NOT.c_%check_stable) THEN |
---|
| 361 | ! WRITE(6,*) " EXECUTION OF THE CODE MUST BE INTERRUPTED AT YOUR REQUEST" |
---|
| 362 | ! WRITE(6,*) " CODE MOST LIKELY DIED IN PURE DA/TPSA/LIE OPERATIONS " |
---|
| 363 | ! STOP 667 |
---|
| 364 | ! ENDIF |
---|
| 365 | ! ENDIF |
---|
| 366 | comp=.true. |
---|
| 367 | if(present(complete)) comp=complete |
---|
| 368 | c_%STABLE_DA =.TRUE. |
---|
| 369 | c_%CHECK_STABLE =.TRUE. |
---|
| 370 | c_%CHECK_MADX_APERTURE =.TRUE. |
---|
| 371 | c_%stable_da =.true. |
---|
| 372 | if(comp) then |
---|
| 373 | xlost=0.0_dp |
---|
| 374 | messagelost=" Aperture has been reset " |
---|
| 375 | nullify(lost_fibre) |
---|
| 376 | nullify(lost_node) |
---|
| 377 | |
---|
| 378 | endif |
---|
| 379 | END SUBROUTINE RESET_APERTURE_FLAG |
---|
| 380 | |
---|
| 381 | SUBROUTINE PRODUCE_APERTURE_FLAG(I) |
---|
| 382 | IMPLICIT NONE |
---|
| 383 | INTEGER I |
---|
| 384 | I=0 |
---|
| 385 | IF(.NOT.c_%CHECK_STABLE) THEN |
---|
| 386 | I=1 |
---|
| 387 | ENDIF |
---|
| 388 | |
---|
| 389 | |
---|
| 390 | ! ALLOW_TRACKING=.TRUE. |
---|
| 391 | |
---|
| 392 | END SUBROUTINE PRODUCE_APERTURE_FLAG |
---|
| 393 | |
---|
| 394 | ! moved here from sa_extend_poly.f90 |
---|
| 395 | |
---|
| 396 | REAL(DP) FUNCTION ROOT(X) ! REPLACES SQRT(X) |
---|
| 397 | IMPLICIT NONE |
---|
| 398 | REAL(DP),INTENT(IN)::X |
---|
| 399 | IF(.NOT.c_%CHECK_STABLE) then |
---|
| 400 | ROOT=1.0_dp |
---|
| 401 | return |
---|
| 402 | endif |
---|
| 403 | |
---|
| 404 | IF((X<0.0_dp).AND.c_%ROOT_CHECK) THEN |
---|
| 405 | ROOT=1.0_dp |
---|
| 406 | c_%CHECK_STABLE=.FALSE. |
---|
| 407 | messagelost="Root undefined " |
---|
| 408 | ELSEIF(X>=0.0_dp) THEN |
---|
| 409 | ROOT=SQRT(X) |
---|
| 410 | ELSE ! IF X IS NOT A NUMBER |
---|
| 411 | ROOT=1.0_dp |
---|
| 412 | c_%CHECK_STABLE=.FALSE. |
---|
| 413 | ENDIF |
---|
| 414 | |
---|
| 415 | END FUNCTION ROOT |
---|
| 416 | |
---|
| 417 | REAL(DP) FUNCTION ARCSIN(X) ! REPLACES ASIN(X) |
---|
| 418 | IMPLICIT NONE |
---|
| 419 | REAL(DP),INTENT(IN)::X |
---|
| 420 | IF(.NOT.c_%CHECK_STABLE) then |
---|
| 421 | ARCSIN=0.0_dp |
---|
| 422 | return |
---|
| 423 | endif |
---|
| 424 | IF((ABS(X)>1.0_dp).AND.c_%ROOT_CHECK) THEN |
---|
| 425 | ARCSIN=0.0_dp |
---|
| 426 | c_%CHECK_STABLE=.FALSE. |
---|
| 427 | messagelost="Arcsin undefined " |
---|
| 428 | ELSEIF(ABS(X)<=1.0_dp) THEN |
---|
| 429 | ARCSIN=ASIN(X) |
---|
| 430 | ELSE ! IF X IS NOT A NUMBER |
---|
| 431 | ARCSIN=0.0_dp |
---|
| 432 | c_%CHECK_STABLE=.FALSE. |
---|
| 433 | ENDIF |
---|
| 434 | |
---|
| 435 | END FUNCTION ARCSIN |
---|
| 436 | |
---|
| 437 | REAL(DP) FUNCTION ARCCOS(X) ! REPLACES ACOS(X) |
---|
| 438 | IMPLICIT NONE |
---|
| 439 | REAL(DP),INTENT(IN)::X |
---|
| 440 | IF(.NOT.c_%CHECK_STABLE) then |
---|
| 441 | ARCCOS=0.0_dp |
---|
| 442 | return |
---|
| 443 | endif |
---|
| 444 | IF((ABS(X)>1.0_dp).AND.c_%ROOT_CHECK) THEN |
---|
| 445 | ARCCOS=0.0_dp |
---|
| 446 | c_%CHECK_STABLE=.FALSE. |
---|
| 447 | messagelost="Arccos undefined " |
---|
| 448 | ELSEIF(ABS(X)<=1.0_dp) THEN |
---|
| 449 | ARCCOS=ACOS(X) |
---|
| 450 | ELSE ! IF X IS NOT A NUMBER |
---|
| 451 | ARCCOS=0.0_dp |
---|
| 452 | c_%CHECK_STABLE=.FALSE. |
---|
| 453 | ENDIF |
---|
| 454 | |
---|
| 455 | END FUNCTION ARCCOS |
---|
| 456 | |
---|
| 457 | REAL(DP) FUNCTION LOGE(X) ! REPLACES ACOS(X) |
---|
| 458 | IMPLICIT NONE |
---|
| 459 | REAL(DP),INTENT(IN)::X |
---|
| 460 | IF(.NOT.c_%CHECK_STABLE) then |
---|
| 461 | LOGE=0.0_dp |
---|
| 462 | return |
---|
| 463 | endif |
---|
| 464 | |
---|
| 465 | IF(X<=0.0_dp.AND.c_%ROOT_CHECK) THEN |
---|
| 466 | LOGE=0.0_dp |
---|
| 467 | c_%CHECK_STABLE=.FALSE. |
---|
| 468 | messagelost="Log undefined " |
---|
| 469 | ELSE |
---|
| 470 | LOGE=LOG(X) |
---|
| 471 | ENDIF |
---|
| 472 | |
---|
| 473 | END FUNCTION LOGE |
---|
| 474 | |
---|
| 475 | |
---|
| 476 | |
---|
| 477 | REAL(DP) FUNCTION COSEH(X) ! REPLACES COSH(X) |
---|
| 478 | IMPLICIT NONE |
---|
| 479 | REAL(DP),INTENT(IN)::X |
---|
| 480 | IF(.NOT.c_%CHECK_STABLE) then |
---|
| 481 | COSEH=1.0_dp |
---|
| 482 | return |
---|
| 483 | endif |
---|
| 484 | |
---|
| 485 | IF((ABS(X)>c_%hyperbolic_aperture).AND.c_%ROOT_CHECK) THEN |
---|
| 486 | COSEH=1.0_dp |
---|
| 487 | c_%CHECK_STABLE=.FALSE. |
---|
| 488 | messagelost="Coseh undefined " |
---|
| 489 | ELSEIF(ABS(X)<=c_%hyperbolic_aperture) THEN |
---|
| 490 | COSEH=COSH(X) |
---|
| 491 | ELSE ! IF X IS NOT A NUMBER |
---|
| 492 | COSEH=1.0_dp |
---|
| 493 | c_%CHECK_STABLE=.FALSE. |
---|
| 494 | ENDIF |
---|
| 495 | |
---|
| 496 | END FUNCTION COSEH |
---|
| 497 | |
---|
| 498 | REAL(DP) FUNCTION SINEH(X) ! REPLACES SINH(X) |
---|
| 499 | IMPLICIT NONE |
---|
| 500 | REAL(DP),INTENT(IN)::X |
---|
| 501 | IF(.NOT.c_%CHECK_STABLE) then |
---|
| 502 | SINEH=0.0_dp |
---|
| 503 | return |
---|
| 504 | endif |
---|
| 505 | |
---|
| 506 | IF((ABS(X)>c_%hyperbolic_aperture).AND.c_%ROOT_CHECK) THEN |
---|
| 507 | SINEH=0.0_dp |
---|
| 508 | c_%CHECK_STABLE=.FALSE. |
---|
| 509 | messagelost="Sineh undefined " |
---|
| 510 | ELSEIF(ABS(X)<=c_%hyperbolic_aperture) THEN |
---|
| 511 | SINEH=SINH(X) |
---|
| 512 | ELSE ! IF X IS NOT A NUMBER |
---|
| 513 | SINEH=0.0_dp |
---|
| 514 | c_%CHECK_STABLE=.FALSE. |
---|
| 515 | ENDIF |
---|
| 516 | |
---|
| 517 | END FUNCTION SINEH |
---|
| 518 | |
---|
| 519 | REAL(DP) FUNCTION arctan(X) ! REPLACES SINH(X) |
---|
| 520 | IMPLICIT NONE |
---|
| 521 | REAL(DP),INTENT(IN)::X |
---|
| 522 | IF(.NOT.c_%CHECK_STABLE) then |
---|
| 523 | arctan=0.0_dp |
---|
| 524 | return |
---|
| 525 | endif |
---|
| 526 | |
---|
| 527 | IF((ABS(X)>c_%hyperbolic_aperture).AND.c_%ROOT_CHECK) THEN |
---|
| 528 | arctan=0.0_dp |
---|
| 529 | c_%CHECK_STABLE=.FALSE. |
---|
| 530 | messagelost="Arctan undefined " |
---|
| 531 | ELSEIF(ABS(X)<=c_%hyperbolic_aperture) THEN |
---|
| 532 | arctan=atan(X) |
---|
| 533 | ELSE ! IF X IS NOT A NUMBER |
---|
| 534 | arctan=0.0_dp |
---|
| 535 | c_%CHECK_STABLE=.FALSE. |
---|
| 536 | ENDIF |
---|
| 537 | |
---|
| 538 | END FUNCTION arctan |
---|
| 539 | |
---|
| 540 | |
---|
| 541 | |
---|
| 542 | end module definition |
---|