source: PSPA/madxPSPA/libs/ptc/src/a_scratch_size.f90 @ 430

Last change on this file since 430 was 430, checked in by touze, 11 years ago

import madx-5.01.00
File size: 39.3 KB
Line 
1!The Full Polymorphic Package
2!Copyright (C) Etienne Forest and Frank Schmidt
3
4! This library is free software; you can redistribute it and/or
5! modify it under the terms of the The Clarified Artistic License.
6! However the source files C_dabnew.f90 and D_Lielib.f90 are derived from
7! the LBNL versions of Berz's DA-package and Forest's analysis package.
8! Their distribution and commercial usage are governed by the laws of
9! the United States of America and more specifically by the USA Department
10! of Energy. The above license may be partly or totally void with respect to
11! these two files.
12
13! All the other files were created as part of Forest's work as an employee
14! of the Japanese Ministry of Culture and Education and are, to best of our
15! knowledge, his intellectual property under Japanese Law.
16
17! THIS PACKAGE IS PROVIDED "AS IS" AND WITHOUT ANY EXPRESS OR IMPLIED WARRANTIES,
18! INCLUDING, WITHOUT LIMITATION, THE IMPLIED WARRANTIES OF MERCHANTIBILITY AND
19! FITNESS FOR A PARTICULAR PURPOSE.
20
21
22module precision_constants
23  implicit none
24  public
25  integer,parameter  :: newscheme_max =200
26  integer,private,parameter::n_read_max=20,NCAR=120
27  private EQUAL_i,EQUAL_Si,EQUAL_r,EQUAL_c !,WRITE_G
28  private read_d,read_int,read_int_a,read_d_a
29  !Double precision
30  real(kind(1d0)) :: doublenum = 0d0
31  !Precision
32  integer,parameter::nlp=24
33  integer,parameter::vp=nlp
34  integer,parameter::lp=4
35  ! double precision
36  integer,parameter::sp=kind(1.e0)
37  integer,parameter::dp=selected_real_kind(2*precision(1.e0))
38  ! quadrupole precision
39  !  integer,parameter::sp=selected_real_kind(2*precision(1.e0))
40  !  integer,parameter::dp=selected_real_kind(4*precision(1.e0))
41  !Logicals
42  logical(lp),parameter:: my_true=.true.
43  logical(lp),parameter:: my_false=.false.
44  logical(lp),target :: global_verbose=.false.
45  logical(lp),target :: no_hyperbolic_in_normal_form=.true.
46  !Numbers double
47  real(dp),parameter::zero=0e0_dp,one=1e0_dp,two=2e0_dp,three=3e0_dp,four=4e0_dp,five=5e0_dp
48  real(dp),parameter::six=6e0_dp,seven=7e0_dp,eight=8e0_dp,nine=9e0_dp,ten=10e0_dp
49  real(dp),parameter::eleven=11e0_dp,twelve=12e0_dp,c_14=14e0_dp,c_15=15e0_dp
50  real(dp),parameter::c_16=16e0_dp,c_20=20e0_dp,c_24=24e0_dp,c_27=27e0_dp,c_30=30e0_dp
51  real(dp),parameter::c_32=32e0_dp,c_40=40e0_dp,c_48=48e0_dp,c_50=50e0_dp
52  real(dp),parameter::c_55=55e0_dp,c_72=72e0_dp,c_80=80e0_dp,c_85=85e0_dp
53  real(dp),parameter::c_90=90e0_dp,c_100=100e0_dp,c_120=120e0_dp,c_160=160e0_dp
54  real(dp),parameter::c_180=180e0_dp,c_1024=1024e0_dp
55  real(dp),parameter::half=0.5e0_dp,c_0_8=0.8e0_dp,c_0_1=0.1e0_dp,c_0_75=0.75e0_dp
56  real(dp),parameter::c_0_4375=0.4375e0_dp,c_4d_1=0.4e0_dp,c_0_125=0.125e0_dp
57  real(dp),parameter::c_1d_2=1e-2_dp,c_1d_3=1e-3_dp,c_0_0001=1e-4_dp
58  real(dp),parameter::c_0_25d_3=0.25e-3_dp,c_0_25=0.25e0_dp,c_0_5d_3=0.5e-3_dp
59  real(dp),parameter::c_1d_8=1e-8_dp,c_1_2=1.2e0_dp,c_0_2=0.2e0_dp
60  real(dp),parameter::c_1d3=1e3_dp,c_1d5=1e5_dp,c_1d6=1e6_dp
61  real(dp),parameter::c_1d9=1e9_dp,c_111110=1.1111e5_dp,c_2d5=2e5_dp
62  real(dp),parameter::c_1d8=1e8_dp,c_1d10=1e10_dp,c_1d4=1e4_dp
63  real(dp),parameter::c_1d30=1e30_dp,c_1d36=1e36_dp
64  real(dp),parameter::c_221=221e0_dp,c_981=981e0_dp,c_867=867e0_dp,c_102=102e0_dp
65  real(dp),parameter::c_183=183e0_dp,c_678=678e0_dp,c_472=472e0_dp,c_66=66e0_dp
66  real(dp),parameter::c_716=716e0_dp,c_2079=2079e0_dp,c_1002=1002e0_dp,c_834=834e0_dp
67  real(dp),parameter::c_454=454e0_dp,c_82=82e0_dp,c_41=41e0_dp,c_216=216e0_dp,c_272=272e0_dp
68  real(dp),parameter::c_840=840e0_dp,c_360=360e0_dp,c_300=300e0_dp,c_137=137e0_dp
69  real(dp),parameter::c_0_28=0.28e0_dp,c_0_31=0.31e0_dp,c_1_8=1.8e0_dp
70  real(dp),parameter::c_1d_5=1e-5_dp, c_0_3079=0.3079e0_dp
71  !Mathematical Constants
72  real(dp),TARGET :: hyperbolic_aperture=10.0_dp
73  real(dp),parameter::pi=3.141592653589793238462643383279502e0_dp,twopi=2.0_dp*pi,pih=pi*0.5_dp
74  real(dp),parameter::twopii=1.0_dp/twopi,pil=pih-0.4_dp,pim=pih+0.4_dp
75  !  real(dp),parameter::rpi4=1.772453850905516027298167e0_dp/two ! rpi4=sqrt(pi)/2
76  real(dp)::rpi4
77  real(dp),parameter::RAD_TO_DEG_=180.0_dp/pi,DEG_TO_RAD_=pi/180.0_dp
78  !Physical Constants
79  real(dp),parameter::A_ELECTRON=1.15965218111e-3_dp  !frs NIST CODATA 2006
80  real(dp),parameter::A_MUON=1.16592069e-3_dp         !frs NIST CODATA 2006
81  real(dp),parameter::A_PROTON=1.79284735e-0_dp       !frs (approx) NIST CODATA 2006
82  real(dp),parameter:: pmaMUON = 105.6583668E-3_DP    !frs NIST CODATA 2006
83  real(dp) :: e_muon = 0.d0
84 !  real(dp),parameter:: pmadt = 1.875612793e0_dp    ! sateesh
85  !  real(dp),parameter:: pmah3 = 2.808391e0_dp    ! sateesh
86  !  real(dp),parameter:: A_dt = -0.142987272e0_dp    ! sateesh
87  !  real(dp),parameter:: a_h3 =-4.183963e0_dp    ! sateesh
88
89
90  real(dp) :: A_particle = A_ELECTRON
91  real(dp),parameter::pmae=5.10998910e-4_dp           !frs NIST CODATA 2006
92  real(dp),parameter::pmae_amu=5.4461702177e-4_dp     !frs NIST CODATA 2006
93  ![GeV]
94  real(dp),parameter::pmap=0.938272013e0_dp           !frs NIST CODATA 2006
95  ![GeV]
96  real(dp),parameter::CLIGHT=2.99792458e8_dp          ! exact
97  ![m/s]
98  real(dp),parameter::hbar=6.58211889e-25_dp          !frs NIST CODATA 2006
99  ![GeV*s]
100  real(dp),parameter::dhbar=1.054571628e-34_dp        !frs NIST CODATA 2006
101  ![J*s]
102  real(dp),parameter::qelect=1.602176487e-19_dp       !frs NIST CODATA 2006
103  ![A*s]
104  real(dp),parameter::eps_0=8.854187817e-12_dp        ! exact
105  ![A*S/V*m]
106  real(dp),parameter::class_e_radius=qelect/4.0_dp/pi/eps_0/pmae/1e9_dp
107  ![m]
108  real(dp),parameter::CGAM=4.0_dp*pi/3.0_dp*class_e_radius/pmae**3
109  ![m/Gev^3] old: 8.846056192e-5_dp
110  real(dp),parameter::HBC=hbar*CLIGHT
111  ![GeV*m] old: HBC=1.9732858e-16_dp
112  !Smallness Parameters
113   !Smallness Parameters
114  real(dp),parameter::epsmac=1e-7_dp,c_1d_20=1e-20_dp,c_1d_37=1e-37_dp
115  real(dp) :: epsflo=1.e-10_dp
116  real(dp),parameter::mybig=1e38_dp
117  real(dp),parameter::machep=1e-17_dp
118  real(dp),parameter::PUNY=1e-38_dp
119  real(dp),parameter::EPSdolmac=1e-7_dp
120  real(dp),parameter::eps_tpsalie=1e-9_dp
121  real(dp),parameter::eps_real_poly=1e-6_dp
122  real(dp),parameter::eps_rot_mis1=1e-6_dp,eps_rot_mis2=1e-9_dp,epsdif=1e-6_dp
123  real(dp),parameter::eps_extend_poly=1e-10_dp
124  real(dp),parameter::eps_fitted=1e-11_dp
125  real(dp),parameter::c_1d_11=1e-11_dp
126  real(dp),parameter::c_1d_16=1e-16_dp
127  real(dp),parameter::eps_def_kind=1e-9_dp
128  real(dp),parameter::deps_tracking=1e-6_dp
129  real(dp),parameter::c_1d_7=1e-7_dp
130  real(dp),parameter::c_1d_38=1e-38_dp
131  real(dp),parameter::c_1d_40=1e-40_dp
132  real(dp),parameter::c_1d_15=1e-15_dp
133  real(dp),parameter::c_1d_6=1e-6_dp
134  real(dp),parameter::c_1d_9=1e-9_dp
135  real(dp),parameter::c_1d_10=1e-10_dp
136  real(dp),parameter::c_2_2d_7=2.2e-7_dp,c_1_35d_8=1.35e-8_dp,c_2_2d_8=2.2e-8_dp
137  real(dp),parameter::c_2_7d_8=2.7e-8_dp
138  !Magic Numbers
139  real(dp),parameter::c_9999_12345=9999.12345e0_dp
140  real(dp),parameter::c_0_284496736=0.284496736e0_dp
141  real(dp),parameter::c_0_254829592=0.254829592e0_dp
142  real(dp),parameter::c_0_3275911=0.3275911e0_dp
143  real(dp),parameter::c_1_061405429=1.061405429e0_dp
144  real(dp),parameter::c_1_421413741=1.421413741e0_dp
145  real(dp),parameter::c_1_453152027=1.453152027e0_dp
146  real(dp),parameter::c_720=720e0_dp,c_30240=30240e0_dp
147  real(dp),parameter::c_1209600=1209600e0_dp,c_21772800=21772800e0_dp
148  real(dp),parameter::c_1_17767998417887=1.17767998417887e0_dp
149  real(dp),parameter::c_0_235573213359357=0.235573213359357e0_dp
150  real(dp),parameter::c_0_78451361047756=0.78451361047756e0_dp
151  real(dp),parameter::c_0_999=0.999e0_dp,c_4_2d_2=4.2e-2_dp
152  real(dp),parameter::c_4_2d_3=4.2e-3_dp,c_6_6d_8=6.6e-8_dp
153  real(dp),parameter::c_2_2d_3=2.2e-3_dp,c_2_2d_4=2.2e-4_dp
154  real(dp),parameter::c_0_9=0.9e0_dp,c_2_5=2.5e0_dp,c_0_148=0.148e0_dp
155  real(dp),parameter::c_0_005=5e-3_dp,c_0_012=1.2e-2_dp,c_1_5=1.5e0_dp
156  real(dp),parameter::c_0_002=2e-3_dp,c_0_05=5e-2_dp,c_0_216=0.216e0_dp
157  real(dp),parameter::c_0_7=0.7e0_dp,c_1_2d_5=1.2e-5_dp,c_1d7=1e7_dp
158  ! Constant Symplectic integrator schemes
159  real(dp) YOSK(0:4), YOSD(4)    ! FIRST 6TH ORDER OF YOSHIDA
160  real(dp),parameter::AAA=-0.25992104989487316476721060727823e0_dp  ! fourth order integrator
161  real(dp),parameter::FD1=0.5_dp/(1.0_dp+AAA),FD2=AAA*FD1,FK1=1.0_dp/(1.0_dp+AAA),FK2=(AAA-1.0_dp)*FK1
162  ! end of symplectic integrator coefficients
163  !Initialized numbers
164  real(dp)::eps=1e-38_dp
165  real(dp)::EPSdol=1e-37_dp
166  LOGICAL(lp),target  :: s_aperture_CHECK=.TRUE.
167  LOGICAL(lp),TARGET  :: ROOT_CHECK=.TRUE.
168  LOGICAL(lp),TARGET  :: CHECK_STABLE=.TRUE.
169  LOGICAL(lp),TARGET  :: WATCH_USER=.FALSE.
170  !  LOGICAL(lp) :: ALLOW_TRACKING=.true.
171  LOGICAL(lp),TARGET  :: CHECK_MADX_APERTURE=.TRUE.
172  LOGICAL(lp),TARGET  :: APERTURE_FLAG=.true.
173
174  REAL(dp),TARGET   :: absolute_aperture=1.0_dp
175  integer,TARGET :: wherelost=0
176  logical(lp),TARGET :: stable_da =.true.
177  logical(lp),TARGET :: check_da =.true.
178  logical(lp),TARGET :: print_frame =.true.
179  logical(lp),TARGET :: sixtrack_compatible =.false.
180  integer ,target ::  spin_normal_position=2
181  real(dp),target ::  da_absolute_aperture=1e6_dp
182  real(dp),pointer :: crash
183  INTEGER,  TARGET :: NPARA_original
184  logical  :: default_tpsa=.false.
185  logical, target :: lingyun_yang=.false.
186  integer, target :: last_tpsa=0
187  integer :: mf_herd=0
188  character*255 :: print_herd="PRINT_HERD.TXT"
189  character*255 :: initial_setting="FINAL_SETTINGS.TXT"
190  character*255 :: final_setting="FINAL_SETTINGS.TXT"
191  character*255 :: def_orbit_node="def_orbit_node.txt"
192  character*255 :: file_block_name="noprint"
193  real(dp) :: lmax=1.e38_dp
194  logical(lp) :: printdainfo=my_false
195  integer   lielib_print(12)
196  DATA lielib_print /0,0,0,0,0,0,0,0,0,0,0,1/
197  INTEGER,TARGET :: SECTOR_NMUL_MAX=10
198  INTEGER, target :: SECTOR_NMUL = 10
199  integer, parameter :: no_e=5  !  electric
200
201  logical(lp) :: change_sector=my_true
202  real(dp) :: xlost(6)=0.0_dp
203  character(255) :: messagelost
204  !  logical(lp) :: fixed_found
205  !  lielib_print(1)=1   lieinit prints info
206  !  lielib_print(2)=1   expflo warning if no convergence
207  !  lielib_print(3)=1   Shows details in flofacg
208  !  lielib_print(4)=1   tunes and damping
209  !  lielib_print(5)=1  order in orbital normal form
210  !  lielib_print(6)=1  symplectic condition
211  !  lielib_print(7)=-1  go manual in normal form  (use auto command in fpp)
212  !  lielib_print(8)=-1  To use nplane from FPP normalform%plane
213  !  lielib_print(9)=1  print in checksymp(s1,norm) in j_tpsalie.f90
214  !  lielib_print(10)=1  print lingyun's checks
215  !  lielib_print(11)=1  print warning about Teng-Edwards
216  !  lielib_print(12)=1  print info in make_node_layout
217
218
219  type info_window
220     character(3) adv
221     integer nc,nr,ni
222     character(NCAR) c(n_read_max)
223     character(NCAR) Fc
224     real(dp) r(n_read_max)
225     character(NCAR) FR
226     integer i(n_read_max)
227     character(NCAR) FI
228  end type info_window
229
230  type CONTROL
231     ! Da stuff
232     real(dp),pointer :: total_da_size !  in megabytes
233     integer,pointer :: lda_used       !  maximum number of da variables in Berz's
234     logical(lp),pointer  :: OLD    ! = true  = bERZ
235     logical(lp),pointer  :: real_warning  ! = true
236     integer,pointer :: no       ! order of da
237     integer,pointer :: nv       ! number of variables
238     integer,pointer :: nd       ! degrees of freedom
239     integer,pointer :: nd2      ! phase space dimension
240     integer,pointer :: np       ! number of parameters in fpp
241     integer,pointer :: nspin       ! number of spin variables (0 or 3)
242     integer,pointer :: SPIN_pos       ! position of spin variables (0 or 3)
243     integer,pointer :: ndpt     ! constant energy variable position is different from zero
244     integer,pointer :: NPARA     ! PARAMETER LOCATION IN PTC in fpp
245     integer,pointer :: npara_fpp     ! PARAMETER LOCATION IN FPP or PTC
246     integer,pointer :: np_pol     ! parameters produced through pol_block
247     logical(lp),pointer :: knob
248     logical(lp),pointer :: valishev
249     !     integer, pointer :: NDPT_OTHER
250     logical(lp), pointer :: setknob
251     REAL(dp),pointer     :: da_absolute_aperture  ! in case one tracks with da.
252     !
253
254     integer,pointer :: wherelost     ! counting lost particles in integration nodes
255     logical(lp),pointer  :: ROOT_CHECK   !=.TRUE. performs check in roots and hyperbolic if true
256     logical(lp),pointer  :: CHECK_STABLE !=.TRUE. particle status
257     logical(lp),pointer  :: CHECK_MADX_APERTURE  !=.TRUE. false means particle lost in aperture
258     logical(lp),pointer  :: APERTURE_FLAG       !=.TRUE. aperture checks globally done (default)
259     logical(lp),pointer  :: s_aperture_CHECK       !=.TRUE. aperture checks globally done (default)
260
261
262     logical(lp),pointer  :: WATCH_USER     ! FALSE NORMALLY : WATCHES USER FOR FAILING TO CHECK APERTURES
263
264     REAL(dp),pointer     :: absolute_aperture     !=1e3_dp generic aperture check
265     real(dp),pointer :: hyperbolic_aperture  ! controls crashes in exponentials
266
267
268     ! influence fibre creation
269
270     integer, pointer :: MADTHICK        !
271     integer, pointer :: MADTHIN_NORMAL
272     integer, pointer :: MADTHIN_SKEW
273     integer, pointer::  NSTD,METD       ! number of steps and integration method
274     logical(lp), pointer :: MADLENGTH   ! =.false. rbend crazy length in mad8 as input
275     logical(lp), pointer :: MAD         !=.false. mad definition of multipole for input only
276     logical(lp), pointer :: EXACT_MODEL != .false. exact model used
277     logical(lp), pointer :: ALWAYS_EXACTMIS  !=.TRUE. exact formula in tracking used for that element
278     logical(lp),pointer :: ALWAYS_knobs  !=.false. ptc knob default status
279     logical(lp),pointer :: recirculator_cheat  ! =.false.  if true energy patches use the time formula always
280     logical(lp),pointer :: sixtrack_compatible !  to insure some sixtrack compatibility default=false
281     integer, pointer:: CAVITY_TOTALPATH ! REAL PILL B0X =1 , FAKE =0  default
282     integer,pointer :: HIGHEST_FRINGE !=2  quadrupole fringe ON IF FRINGE PRESENT
283     logical(lp),pointer :: do_beam_beam   ! obvious meaning: false normally
284     ! creates a reverse propagator
285     integer,pointer ::FIBRE_DIR         !=1 or -1 for reversed
286     real(dp),pointer ::INITIAL_CHARGE         ! =1 or -1 AND  ADJUST THE MASS IS THE PREFERED MODE
287     ! creates a reverse propagator and a reversed ring in combination with above
288     logical(lp),pointer ::FIBRE_flip    !=.true.
289     !  x_prime true means noncanonical outside magnets. x(5) variables stays the same.
290     real(dp), pointer :: eps_pos
291     ! fill once and never touch again
292
293     integer, pointer :: SECTOR_NMUL_MAX     != 10 maxwell equations is solved to order 10 in exact sectors
294     integer, pointer :: SECTOR_NMUL     != 4  MULTIPOLES IN TEAPOT BEND ALLOWED BY DEFAULT
295     real(dp), pointer :: wedge_coeff(:)     ! QUAD_KICK IN WEDGE
296     logical(lp), pointer :: MAD8_WEDGE      ! QUAD_KICK + FRINGE IF FRINGE IS OUT.
297
298     logical(lp), pointer:: electron     !  electron if true otherwise proton
299     real(dp), pointer :: massfactor     !=one  sets variable muon and electron must be true
300     ! global on the fly
301     logical(lp), pointer :: compute_stoch_kick != .false. store stochastic kick for stochastic tracking
302     logical(lp),pointer :: FEED_P0C   !=.FALSE.  work takes p0c instead of energy
303     logical(lp),pointer :: ALWAYS_EXACT_PATCHING  !=.TRUE. patching done correctly
304     ! used to output horror messages
305     logical(lp),pointer :: stable_da !=.true.  interrupts DA if check_da is true
306     logical(lp),pointer :: check_da  !=.true.
307     logical(lp),pointer :: OLD_IMPLEMENTATION_OF_SIXTRACK  !=.true.
308     real(dp),pointer :: phase0 ! default phase in cavity
309     logical(lp), pointer :: global_verbose
310     logical(lp), pointer :: no_hyperbolic_in_normal_form ! unstable produces exception
311  end type CONTROL
312
313  type(control) c_
314
315  type(info_window),TARGET:: w_i
316  type(info_window),TARGET:: w_ii
317  type(info_window),TARGET::  r_i
318  type(info_window),pointer:: W_P
319  type(info_window),pointer:: R_P
320
321  INTERFACE assignment (=)
322     MODULE PROCEDURE EQUAL_Si
323     MODULE PROCEDURE EQUAL_i
324     MODULE PROCEDURE EQUAL_r
325     MODULE PROCEDURE EQUAL_c
326  end  INTERFACE
327
328  !  INTERFACE ! WRITE_I
329  !     MODULE PROCEDURE WRITE_G  ! not private
330  !  END INTERFACE
331  !  INTERFACE WRITE_E
332  !     MODULE PROCEDURE WRITE_G  ! not private
333  !  END INTERFACE
334
335  INTERFACE read
336     MODULE PROCEDURE read_d
337     MODULE PROCEDURE read_int
338     MODULE PROCEDURE read_int_a
339     MODULE PROCEDURE read_d_a
340  END INTERFACE
341
342contains
343
344  function mat_norm(m)
345    implicit none
346    real(dp) mat_norm
347    real(dp) m(:,:)
348    integer i,j
349
350    mat_norm=0.0_dp
351    do i=1,size(m,dim=1)
352       do j=1,size(m,dim=2)
353          mat_norm=mat_norm+abs(m(i,j))
354       enddo
355    enddo
356
357  end function mat_norm
358
359  SUBROUTINE  check_stability(S1)
360    implicit none
361    REAL(DP),INTENT(INOUT)::S1(6)
362    INTEGER I
363
364    DO I=1,5
365       IF(ABS(S1(I))>C_%ABSOLUTE_APERTURE) THEN
366          S1=PUNY
367          C_%CHECK_STABLE=.FALSE.
368          EXIT
369       ENDIF
370    ENDDO
371
372  end   SUBROUTINE check_stability
373
374  ! Symplectic integrator routines setting coefficients
375
376  SUBROUTINE MAKE_YOSHIDA
377    IMPLICIT NONE
378    integer i
379
380    YOSK(4)=0.0_dp
381    YOSK(3)=0.78451361047756e0_dp
382    YOSK(2)=0.235573213359357e0_dp
383    YOSK(1)=-1.17767998417887e0_dp
384    YOSK(0)=1.0_dp-2.0_dp*(YOSK(1)+YOSK(2)+YOSK(3))
385
386    do i=4,1,-1
387       YOSD(i)=(YOSK(i)+YOSK(i-1))/2.0_dp
388    enddo
389
390    do i=3,0,-1
391       YOSK(i+1)=YOSK(I)
392    enddo
393
394
395  END SUBROUTINE MAKE_YOSHIDA
396
397  SUBROUTINE input_sector(se2,se1)
398    implicit none
399    logical ttt,t1,t2
400    integer se1,se2
401
402    !    if(.not.change_sector) return
403
404    t1=(SECTOR_NMUL_MAX/=se1)
405    t2=(SECTOR_NMUL/=se2)
406
407    ttt=t1.or.t2
408
409    if(ttt) then
410       if(change_sector) then
411          write(6,*) " SECTOR_NMUL_MAX is changed from ",SECTOR_NMUL_MAX," to ",se1
412          write(6,*) " SECTOR_NMUL is changed from ",SECTOR_NMUL," to ",se2
413          write(6,*) " GLOBAL VARIABLES that can no longer be changed"
414          SECTOR_NMUL_MAX=se1
415          SECTOR_NMUL=se2
416       else
417          if(t1) write(6,*) " sector_nmul_max CANNOT be changed from ",SECTOR_NMUL_MAX," to ",se1
418          if(t2) write(6,*) " sector_nmul CANNOT be changed from ",SECTOR_NMUL," to ",se2
419          write(6,*) " Watch out : The are GLOBAL VARIABLES "
420       endif
421    endif
422
423  end subroutine input_sector
424
425  !  Printing routines for screen
426
427  subroutine  get_ncar(n)
428    implicit none
429    integer n
430    n=ncar
431  end subroutine  get_ncar
432
433  SUBROUTINE  EQUAL_Si(S2,S1)
434    implicit none
435    type (info_window),INTENT(inOUT)::S2
436    integer,INTENT(IN)::S1
437    integer i
438    if(s1==1)then
439       s2%adv='NO'
440    else
441       s2%adv='YES'
442    endif
443    s2%ni=0
444    s2%nC=0
445    s2%nR=0
446    s2%Fi=' '
447    s2%FC=' '
448    s2%FR=' '
449    do i=1,n_read_max
450       s2%i(i)=0
451       s2%R(i)=0.0_dp
452       s2%C(i)=' '
453    enddo
454  END SUBROUTINE EQUAL_Si
455
456  SUBROUTINE  EQUAL_i(S2,S1)
457    implicit none
458    type (info_window),INTENT(inOUT)::S2
459    integer,INTENT(IN)::S1(:)
460    integer n,i
461    n=size(s1)
462    s2%ni=n
463    do i=1,n
464       s2%i(i)=s1(i)
465    enddo
466  END SUBROUTINE EQUAL_i
467
468  SUBROUTINE  EQUAL_r(S2,S1)
469    implicit none
470    type (info_window),INTENT(inOUT)::S2
471    real(dp),INTENT(IN)::S1(:)
472    integer n,i
473    n=size(s1)
474    s2%nr=n
475    do i=1,n
476       s2%r(i)=s1(i)
477    enddo
478  END SUBROUTINE EQUAL_r
479
480  SUBROUTINE  EQUAL_c(S2,S1)
481    implicit none
482    type (info_window),INTENT(inOUT)::S2
483    character(*),INTENT(IN)::S1(*)
484    integer i
485    do i=1,s2%nc
486       s2%c(i)=' '
487       s2%c(i)=s1(i)
488    enddo
489  END SUBROUTINE EQUAL_c
490
491  !  SUBROUTINE WRITE_G(IEX)
492  !    IMPLICIT NONE
493  !    integer, OPTIONAL :: IEX
494  !    integer I,MYPAUSE,IPAUSE
495  !    if(.not.global_verbose) return
496  !    IF(W_P%NC/=0) THEN
497  !       if(W_P%FC/=' ') then
498  !          WRITE(6,W_P%FC,advance=W_P%ADV) (W_P%C(I), I=1,W_P%NC)
499  !       else
500  !          do i=1,W_P%NC
501  !             WRITE(6,*) W_P%C(I)
502  !          enddo
503  !       endif
504  !    ENDIF
505  !    IF(W_P%NI/=0) THEN
506  !       if(W_P%FI/=' ') then
507  !          WRITE(6,W_P%FI,advance=W_P%ADV) (W_P%I(I), I=1,W_P%NI )
508  !       else
509  !          do i=1,W_P%NI
510  !             WRITE(6,*) W_P%I(I)
511  !          enddo
512  !       endif
513  !    ENDIF
514  !    IF(W_P%NR/=0) THEN
515  !       if(W_P%FR/=' ') then
516  !          WRITE(6,W_P%FR,advance=W_P%ADV) (W_P%R(I), I=1,W_P%NR)
517  !       else
518  !          do i=1,W_P%NR
519  !             WRITE(6,*) W_P%R(I)
520  !          enddo
521  !       endif
522  !    ENDIF
523  !    if(W_P%ADV=='NO') then
524  !       WRITE(6,*) " "
525  !    endif
526  !    IF(PRESENT(IEX)) THEN
527  !       if(iex==-1) stop
528  !       IPAUSE=MYPAUSE(IEX)
529  !    ENDIF
530  !  END SUBROUTINE WRITE_G
531  !
532  !  SUBROUTINE WRITE_a(IEX)
533  !    IMPLICIT NONE
534  !    integer, OPTIONAL :: IEX
535  !    logical(lp) temp
536  !    temp=global_verbose
537  !    global_verbose=.true.
538  !    ! call ! WRITE_I(IEX)
539  !    global_verbose=temp
540  !
541  !  END SUBROUTINE WRITE_a
542  !
543  SUBROUTINE read_int(IEX)
544    IMPLICIT NONE
545    integer, intent(inout):: iex
546    read(5,*)  iex
547  END SUBROUTINE read_int
548
549  SUBROUTINE read_int_a(IEX,n)
550    IMPLICIT NONE
551    integer, intent(inout):: iex(:)
552    integer, intent(in):: n
553    integer i
554    read(5,*) (iex(i),i=1,n)
555  END SUBROUTINE read_int_a
556
557  SUBROUTINE read_d(IEX)
558    IMPLICIT NONE
559    real(dp), intent(inout)::   iex
560    read(5,*) iex
561  END SUBROUTINE read_d
562
563  SUBROUTINE read_d_a(IEX,n)
564    IMPLICIT NONE
565    real(dp), intent(inout)::   iex(:)
566    integer, intent(in):: n
567    integer i
568    read(5,*) (iex(i),i=1,n)
569  END SUBROUTINE read_d_a
570
571!!!!!!!!!!!!!!!!!! old   special for lielib: others in sa_extend_poly
572  REAL(DP) FUNCTION  ARCCOS_lielib(X)  ! REPLACES ACOS(X)
573    IMPLICIT NONE
574    REAL(DP),INTENT(IN)::X
575    IF(.NOT.c_%CHECK_STABLE) then
576       ARCCOS_lielib=0.0_dp
577       return
578    endif
579    IF((ABS(X)>1.0_dp).AND.c_%ROOT_CHECK) THEN
580       ARCCOS_lielib=0.0_dp
581       c_%CHECK_STABLE=.FALSE.
582    ELSEIF(ABS(X)<=1.0_dp) THEN
583       ARCCOS_lielib=ACOS(X)
584    ELSE      !  IF X IS NOT A NUMBER
585       ARCCOS_lielib=0.0_dp
586       c_%CHECK_STABLE=.FALSE.
587    ENDIF
588
589  END FUNCTION ARCCOS_lielib
590
591  REAL(DP) FUNCTION  LOGE_lielib(X)  ! REPLACES ACOS(X)
592    IMPLICIT NONE
593    REAL(DP),INTENT(IN)::X
594    IF(.NOT.c_%CHECK_STABLE) then
595       LOGE_lielib=0.0_dp
596       return
597    endif
598
599    IF(X<=0.0_dp.AND.c_%ROOT_CHECK) THEN
600       LOGE_lielib=0.0_dp
601       c_%CHECK_STABLE=.FALSE.
602    ELSE
603       LOGE_lielib=LOG(X)
604    ENDIF
605
606  END FUNCTION LOGE_lielib
607
608end module precision_constants
609
610
611module scratch_size
612  implicit none
613  public
614  integer,parameter::ndumt=10                   !  Number of scratch level
615end module scratch_size
616
617module file_handler
618  use precision_constants
619  implicit none
620  public
621  private INTFILE,INTFILE_K
622  integer,private,parameter::nf=3,MFI=20,MFO=99
623  integer,private :: winterfile(nf) =(/40,41,42/)
624  logical(lp) :: myfile(MFI:MFO)
625  logical(lp),private :: doit=.true.,ginofile=.false.
626
627  type file_
628     logical(lp) MF
629     ! AIMIN CHANGES FOR MS4.0
630     !   logical(lp) :: mf=.false.
631  end type file_
632
633  type file_K
634     logical(lp) MF
635     ! AIMIN CHANGES FOR MS4.0
636     !  logical(lp) :: mf=.false.
637  end type file_K
638
639  type(FILE_)  NEWFILE
640  type(FILE_K)  closeFILE
641
642
643  INTERFACE ASSIGNMENT (=)
644     MODULE PROCEDURE INTFILE
645     MODULE PROCEDURE INTFILE_K
646  END  INTERFACE
647
648CONTAINS
649
650
651  SUBROUTINE  INTFILE(S2,S1)
652    ! gino module here
653    implicit none
654    integer,INTENT(inout)::S2
655    type (FILE_),INTENT(in)::S1
656    integer i
657
658    if(ginofile) then
659       ! put gino compatible code here
660
661    else
662       if(doit) then
663          call zerofile
664          doit=.false.
665       endif
666       I=MFI
667       DO WHILE(myfile(i).AND.I<=MFO)
668          I=I+1
669       ENDDO
670       IF(I==MFO) THEN
671          !       w_p=0
672          !       w_p%nc=1
673          !       w_p%fc='(1(1X,A120))'
674          !       w_p%c(1)=  " NO MORE UNITS AVAILABLE: INTFILE "
675          !       ! call !write_e(MFO)
676       ENDIF
677       S2=I
678       myfile(I)=.true.
679       if(s1%mf) then
680          !      w_p=0
681          !      w_p%nc=1
682          !      w_p%fc='(1(1X,A120))'
683          !      write(w_p%c(1),'(1x,L1,1x)')   s1%mf
684          !      ! call ! WRITE_I
685       endif
686    endif
687  END SUBROUTINE INTFILE
688
689
690  SUBROUTINE  INTFILE_K(S2,S1)
691    implicit none
692    integer,INTENT(inOUT)::S2
693    type (FILE_K),INTENT(IN)::S1
694    integer IPAUSE,MYPAUSES
695    character(120) line
696
697
698    if(ginofile) then
699       ! put gino compatible code here
700
701    else
702       IF(S2>=MFI.AND.S2<MFO) THEN
703          IF(myfile(S2)) THEN
704             myfile(S2)=.false.
705          ELSE
706             WRITE(line,'(a30)') "PROBLEMS WITH UNITS: INTFILE_K"
707             WRITE(line(31:120),'(1x,i4,1x,L1)') S2, myfile(S2)
708             IPAUSE=MYPAUSES(1,line)
709          ENDIF
710       ELSE
711          WRITE(line(1:30),'(a30)') "PROBLEMS WITH UNITS: INTFILE_K"
712          WRITE(line(31:120),'(1x,i4,1x,L1)') S2, myfile(S2)
713          IPAUSE=MYPAUSES(2,line)
714       ENDIF
715       CLOSE(S2)
716       S2=-S2
717       if(s1%mf) THEN
718          WRITE(line,'(a9,L1)') " s1%mf = ",s1%mf
719          IPAUSE=MYPAUSES(3,line)
720       ENDIF
721    endif
722  END SUBROUTINE INTFILE_K
723
724
725  SUBROUTINE  ZEROFILE
726    implicit none
727    integer i
728
729    DO I=MFI,MFO
730       myfile(i)=.FALSE.
731    ENDDO
732    do i=1,nf
733       myfile(winterfile(i))=.true.
734    enddo
735  END SUBROUTINE ZEROFILE
736
737  SUBROUTINE KanalNummer(iff,file)
738    implicit none
739    integer, INTENT(OUT) :: iff
740    character(*),optional :: file
741    logical :: opened, exists
742    integer :: i
743
744    DO i= 9999, 7, -1
745       INQUIRE(UNIT= i, EXIST= exists, OPENED= opened)
746       IF (exists .AND. (.NOT. opened)) GOTO 20
747    ENDDO
748    WRITE (UNIT= *, FMT= *) ' cannot find free unit within the range 7-9999..'
749    CALL ReportOpenFiles
750    STOP
75120  CONTINUE
752    iff= I
753
754    if(present(file)) then
755       open(unit=iff,file=file)
756    endif
757  END SUBROUTINE KanalNummer
758
759  SUBROUTINE ReportOpenFiles
760    implicit none
761
762    logical :: opened, exists, named
763    CHARACTER(LEN= 400) :: name
764    integer :: i
765
766    DO i= 9999, 7, -1
767       INQUIRE(UNIT= i, EXIST= exists, OPENED= opened)
768       IF (exists .AND. opened) THEN
769          INQUIRE(UNIT= i, NAMED= named, NAME= name)
770          write (*,7010) i, name(:LEN_TRIM(name))
771       ENDIF
772    ENDDO
7737010 FORMAT(' iUnit:',I3,', name: "',A,'"')
774
775  END SUBROUTINE ReportOpenFiles
776
777
778  SUBROUTINE CONTEXT( STRING, nb )
779    IMPLICIT NONE
780    CHARACTER(*) STRING
781    CHARACTER(1) C1
782    integer, optional :: nb
783    integer I,J,K,nb0,count
784    nb0=0
785    if(present(nb)) nb0=1
786    J = 0
787    count=0
788    DO I = 1, LEN (STRING)
789       C1 = STRING(I:I)
790       STRING(I:I) = ' '
791       IF( C1 .NE. ' ' ) THEN
792          if(count/=0.and.nb0==1) then
793             J = J + 1
794             STRING(J:J) = ' '
795             count=0
796          endif
797          J = J + 1
798          K = ICHAR( C1 )
799          IF( K .GE. ICHAR('a') .AND. K .LE. ICHAR('z') ) THEN
800             C1 = CHAR( K - ICHAR('a') + ICHAR('A') )
801          ENDIF
802          STRING(J:J) = C1
803       else
804          count=count+1
805       ENDIF
806    ENDDO
807    string=string(1:len_trim(string))
808    RETURN
809  END  SUBROUTINE CONTEXT
810
811  SUBROUTINE create_name(name_root,ind,suffix,filename)
812    implicit none
813    CHARACTER(*) name_root,suffix,filename
814    integer :: ind
815    character(20) temp
816    call context(name_root)
817    call context(suffix)
818
819    write(temp,*) ind
820    call context(temp)
821
822    filename=' '
823
824    filename=name_root(1:len_trim(name_root))//temp(1:len_trim(temp))//'.'//suffix(1:len_trim(suffix))
825
826  END SUBROUTINE create_name
827
828end module file_handler
829
830
831module my_own_1D_TPSA
832  USE precision_constants
833  implicit none
834  !  public
835  private input_real_in_my_1D_taylor
836  integer :: n_tpsa_exp = 10
837  integer, parameter :: N_my_1D_taylor=31  ! SHOULD BE AS ENGE_N
838  integer, private :: No_my_1D_taylor=N_my_1D_taylor  ! SHOULD BE AS ENGE_N
839
840  private mul,dmulsc,dscmul,INV
841  private div,ddivsc,dscdiv,Idivsc
842  private add,unaryADD,daddsc,dscadd
843  private subs,unarySUB,dsubsc,dscsub,POW
844  private  DEXPT,DLOGT,DSQRTT
845  private  DCOST,DSINT
846
847  ! Definitions
848  type my_1D_taylor
849     real(dp) a(0:N_my_1D_taylor)
850
851  END type my_1D_taylor
852
853  INTERFACE assignment (=)
854     MODULE PROCEDURE input_my_1D_taylor_in_real   !@1 &nbsp; Taylor=real
855     MODULE PROCEDURE input_real_in_my_1D_taylor  !@1 &nbsp; real=Taylor
856  end  INTERFACE
857
858  INTERFACE OPERATOR (*)
859     MODULE PROCEDURE mul   !@1 &nbsp; Taylor * Taylor
860     MODULE PROCEDURE dmulsc    !@1 &nbsp; Taylor * Real(dp)
861     MODULE PROCEDURE dscmul     !@1 &nbsp;  Real(dp) * Taylor
862  END INTERFACE
863
864  INTERFACE OPERATOR (/)
865     MODULE PROCEDURE div    !@1 &nbsp; Taylor / Taylor
866     MODULE PROCEDURE ddivsc  !@1 &nbsp; Taylor / Real(dp)
867     MODULE PROCEDURE dscdiv  !@1 &nbsp; Real(dp) / Taylor
868     MODULE PROCEDURE Idivsc  !@1 &nbsp; Taylor / Integer  <font color="#FF0000">&nbsp; &#8594; &nbsp; added because useful in example code</font>
869  END INTERFACE
870
871  INTERFACE OPERATOR (+)
872     MODULE PROCEDURE add       !@1 &nbsp; Taylor + Taylor
873     MODULE PROCEDURE unaryADD  !@1 &nbsp;  +Taylor
874     MODULE PROCEDURE daddsc    !@1 &nbsp;  Taylor + Real(dp)
875     MODULE PROCEDURE dscadd    !@1 &nbsp;  Real(dp) + Taylor
876  END INTERFACE
877
878  INTERFACE OPERATOR (-)
879     MODULE PROCEDURE subs      !@1 &nbsp; Taylor - Taylor
880     MODULE PROCEDURE unarySUB     !@1 &nbsp;  -Taylor
881     MODULE PROCEDURE dsubsc   !@1 &nbsp;  Taylor - Real(dp)
882     MODULE PROCEDURE dscsub    !@1 &nbsp;  Real(dp) - Taylor
883  END INTERFACE
884
885
886  INTERFACE OPERATOR (**)
887     MODULE PROCEDURE POW    !@1 &nbsp; Taylor ** Integer
888  END INTERFACE
889
890
891  !@3  <p><i><font size="5">&nbsp;Overloading standard procedures </font></i></p>
892
893  INTERFACE exp
894     MODULE PROCEDURE DEXPT   !@1 &nbsp; exp(Taylor)
895  END INTERFACE
896
897  INTERFACE LOG
898     MODULE PROCEDURE DLOGT   !@1 &nbsp; log(Taylor)
899  END INTERFACE
900
901  INTERFACE SQRT
902     MODULE PROCEDURE DSQRTT    !@1 &nbsp; sqrt(Taylor)
903  END INTERFACE
904
905  INTERFACE COS
906     MODULE PROCEDURE DCOST    !@1 &nbsp; cos(Taylor)
907  END INTERFACE
908
909  INTERFACE SIN
910     MODULE PROCEDURE DSINT   !@1 &nbsp; sin(Taylor)
911  END INTERFACE
912
913
914
915  !@3  <p><i><font size="5">User defined operator</font></i></p>
916
917
918  ! Destructors and Constructors for my_1D_taylor
919
920
921contains
922
923  subroutine set_my_taylor_no(no)
924    implicit none
925    integer no
926
927    if(no>N_my_1D_taylor) then
928       No_my_1D_taylor=N_my_1D_taylor
929       write(6,*) " warning NO too big in set_my_taylor_no: recompile FPP if needed "
930    else
931       No_my_1D_taylor=no
932    endif
933
934  end subroutine set_my_taylor_no
935
936  FUNCTION add( S1, S2 )
937    implicit none
938    type (my_1D_taylor) add
939    type (my_1D_taylor), INTENT (IN) :: S1, S2
940
941    add%a=S1%a + S2%a
942
943  END FUNCTION add
944
945  FUNCTION daddsc( S1, sc )
946    implicit none
947    type (my_1D_taylor) daddsc
948    type (my_1D_taylor), INTENT (IN) :: S1
949    real(dp), INTENT (IN) :: sc
950
951    daddsc=s1
952    daddsc%a(0)= s1%a(0) + sc
953
954  END FUNCTION daddsc
955
956  FUNCTION dscadd( sc ,  S1)
957    implicit none
958    type (my_1D_taylor) dscadd
959    type (my_1D_taylor), INTENT (IN) :: S1
960    real(dp), INTENT (IN) :: sc
961
962    dscadd=s1
963    dscadd%a(0)= s1%a(0) + sc
964
965  END FUNCTION dscadd
966
967  FUNCTION unaryadd( S1 )
968    implicit none
969    type (my_1D_taylor) unaryadd
970    type (my_1D_taylor), INTENT (IN) :: S1
971
972    unaryadd=s1
973
974  END FUNCTION unaryadd
975
976
977  FUNCTION subs( S1, S2 )
978    implicit none
979    type (my_1D_taylor) subs
980    type (my_1D_taylor), INTENT (IN) :: S1, S2
981
982    subs%a=S1%a - S2%a
983
984  END FUNCTION subs
985
986  FUNCTION unarySUB( S1 )
987    implicit none
988    type (my_1D_taylor) unarySUB
989    type (my_1D_taylor), INTENT (IN) :: S1
990
991    unarySUB%a=-s1%a
992
993  END FUNCTION unarySUB
994
995  FUNCTION dsubsc( S1, sc )
996    implicit none
997    type (my_1D_taylor) dsubsc
998    type (my_1D_taylor), INTENT (IN) :: S1
999    real(dp), INTENT (IN) :: sc
1000
1001    dsubsc=s1
1002    dsubsc%a(0)= s1%a(0) - sc
1003
1004  END FUNCTION dsubsc
1005
1006
1007  FUNCTION dscsub( sc , S1 )
1008    implicit none
1009    type (my_1D_taylor) dscsub
1010    type (my_1D_taylor), INTENT (IN) :: S1
1011    real(dp), INTENT (IN) :: sc
1012
1013    dscsub=-s1   ! uses unary sub
1014    dscsub=sc + dscsub   ! uses add
1015
1016  END FUNCTION dscsub
1017
1018
1019  FUNCTION mul( S1, S2 )
1020    implicit none
1021    type (my_1D_taylor) mul
1022    type (my_1D_taylor), INTENT (IN) :: S1, S2
1023    integer I,J
1024    mul%a=0.0_dp
1025    DO I=0,No_my_1D_taylor
1026       DO J=0,No_my_1D_taylor
1027          IF(I+J>No_my_1D_taylor) CYCLE
1028          mul%a(I+J)=S1%a(I)*S2%a(J)+ mul%a(I+J)
1029       ENDDO
1030    ENDDO
1031
1032  END FUNCTION mul
1033
1034  FUNCTION dmulsc( S1, sc )
1035    implicit none
1036    type (my_1D_taylor) dmulsc
1037    type (my_1D_taylor), INTENT (IN) :: S1
1038    real(dp), INTENT (IN) :: sc
1039
1040    dmulsc%a= s1%a*sc
1041
1042  END FUNCTION dmulsc
1043
1044  FUNCTION dscmul( sc ,S1 )
1045    implicit none
1046    type (my_1D_taylor) dscmul
1047    type (my_1D_taylor), INTENT (IN) :: S1
1048    real(dp), INTENT (IN) :: sc
1049
1050    dscmul%a= s1%a*sc
1051
1052  END FUNCTION dscmul
1053
1054  FUNCTION ddivsc( S1, sc )
1055    implicit none
1056    type (my_1D_taylor) ddivsc
1057    type (my_1D_taylor), INTENT (IN) :: S1
1058    real(dp), INTENT (IN) :: sc
1059
1060    ddivsc%a= s1%a/sc
1061
1062  END FUNCTION ddivsc
1063
1064  FUNCTION Idivsc( S1, sc )
1065    implicit none
1066    type (my_1D_taylor) Idivsc
1067    type (my_1D_taylor), INTENT (IN) :: S1
1068    integer, INTENT (IN) :: sc
1069
1070    Idivsc%a= s1%a/sc
1071
1072  END FUNCTION Idivsc
1073
1074  FUNCTION POW( S1,N )
1075    implicit none
1076    type (my_1D_taylor) POW , T
1077    type (my_1D_taylor), INTENT (IN) :: S1
1078    integer, INTENT (IN) :: N
1079    integer I
1080
1081    POW=1.0_dp
1082
1083    IF(N<=0) THEN
1084       T=1.0_dp/S1
1085       DO I=1,-N
1086          POW=POW*T
1087       ENDDO
1088    ELSE
1089       DO I=1,N
1090          POW=POW*S1
1091       ENDDO
1092
1093    ENDIF
1094
1095  END FUNCTION POW
1096
1097  FUNCTION inv( S1 )
1098    implicit none
1099    type (my_1D_taylor) inv,t,TT
1100    type (my_1D_taylor), INTENT (IN) :: S1
1101    integer I
1102    T=S1/S1%A(0)
1103    T%A(0)=0.0_dp
1104
1105    TT=T
1106    inv=1.0_dp
1107
1108    DO I=1,No_my_1D_taylor
1109       INV=INV-TT
1110       TT=-TT*T
1111    ENDDO
1112
1113    INV=INV/S1%A(0)
1114
1115  END FUNCTION inv
1116
1117  FUNCTION DIV( S1, S2 )
1118    implicit none
1119    type (my_1D_taylor) DIV
1120    type (my_1D_taylor), INTENT (IN) :: S1, S2
1121
1122    DIV=INV(S2)
1123    DIV=S1*DIV
1124  END FUNCTION DIV
1125
1126  FUNCTION dscdiv(  sc , S1)
1127    implicit none
1128    type (my_1D_taylor) dscdiv
1129    type (my_1D_taylor), INTENT (IN) :: S1
1130    real(dp), INTENT (IN) :: sc
1131
1132    dscdiv=INV(S1)
1133    dscdiv= SC *  dscdiv
1134
1135  END FUNCTION dscdiv
1136
1137  !  DEFININING my_1D_taylor=CONSTANT
1138  subroutine input_real_in_my_1D_taylor( S2, S1 )
1139    implicit none
1140    real(dp), INTENT (IN) :: S1
1141    type (my_1D_taylor), INTENT (inout) :: S2
1142
1143    S2%a=0.0_dp
1144    S2%a(0)=s1
1145
1146  END subroutine input_real_in_my_1D_taylor
1147
1148  subroutine input_my_1D_taylor_in_real( S2, S1 )
1149    implicit none
1150    type (my_1D_taylor), INTENT (IN) :: S1
1151    real(dp), INTENT (inout) :: S2
1152
1153    S2=S1%a(0)
1154
1155  END subroutine input_my_1D_taylor_in_real
1156
1157  !  DEFININING EXP
1158
1159  FUNCTION DEXPT( S1 )
1160    implicit none
1161    type (my_1D_taylor) DEXPT,t,tt
1162    type (my_1D_taylor), INTENT (IN) :: S1
1163    integer I
1164    T=S1
1165    T%A(0)=0.0_dp
1166
1167
1168    DEXPT=1.0_dp
1169    tt=1.0_dp
1170
1171    do i=1,No_my_1D_taylor
1172       tt=tt*t/i
1173       DEXPT=DEXPT + tt
1174    enddo
1175
1176    DEXPT=DEXPT*EXP(S1%A(0))
1177
1178  END FUNCTION DEXPT
1179
1180
1181  FUNCTION DLOGT( S1 )
1182    implicit none
1183    type (my_1D_taylor) DLOGT,t,TT
1184    type (my_1D_taylor), INTENT (IN) :: S1
1185    integer I
1186
1187    T=S1/S1%A(0)
1188    T%A(0)=0.0_dp
1189    DLOGT=0.0_dp
1190    TT=T
1191    do i=1,No_my_1D_taylor
1192       DLOGT=TT/I+DLOGT
1193       TT=-TT*T
1194    ENDDO
1195
1196    DLOGT=DLOGT+ LOG(S1%A(0))
1197
1198  END FUNCTION  DLOGT
1199
1200  FUNCTION DSQRTT( S1 )
1201    implicit none
1202    type (my_1D_taylor) DSQRTT,t
1203    type (my_1D_taylor), INTENT (IN) :: S1
1204    ! WRITE(6,*) " MARDE "
1205    ! STOP 666
1206    ! T=S1/S1%A(0)
1207    ! T%A(0)=zero
1208
1209    ! DSQRTT=one + T/two - T**2/eight+T**3/c_16
1210    ! DSQRTT=DSQRTT* SQRT(S1%A(0))
1211    DSQRTT=log(s1)/2.0_dp
1212    DSQRTT=exp(DSQRTT)
1213  END FUNCTION  DSQRTT
1214
1215  FUNCTION DCOST(S1)
1216    implicit none
1217    type (my_1D_taylor) DCOST,t
1218    type (my_1D_taylor), INTENT (IN) :: S1
1219    WRITE(6,*) " MARDE "
1220    STOP 666
1221
1222    T=S1
1223    T%A(0)=0.0_dp
1224
1225    DCOST=COS(S1%A(0))*(1.0_dp-T**2/2.0_dp)-SIN(S1%A(0))*(T-T**3/6.0_dp)
1226
1227  END FUNCTION  DCOST
1228
1229  FUNCTION DSINT(S1)
1230    implicit none
1231    type (my_1D_taylor) DSINT,t
1232    type (my_1D_taylor), INTENT (IN) :: S1
1233
1234    T=S1
1235    T%A(0)=0.0_dp
1236
1237    DSINT=SIN(S1%A(0))*(1.0_dp-T**2/2.0_dp)+COS(S1%A(0))*(T-T**3/6.0_dp)
1238
1239  END FUNCTION  DSINT
1240
1241end module my_own_1D_TPSA
1242
1243module gauss_dis
1244
1245  use precision_constants
1246  public
1247  private alex_ISEED
1248  integer :: alex_ISEED=1000
1249  !cccccccccccccccccccccccc     Program 2.13     ccccccccccccccccccccccccc
1250  !
1251  !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1252  !                                                                      c
1253  ! Please Note:                                                         c
1254  !                                                                      c
1255  ! (1) This computer program is part of the book, "An Introduction to   c
1256  !     Computational Physics," written by Tao Pang and published and    c
1257  !     copyrighted by Cambridge University Press in 1997.               c
1258  !                                                                      c
1259  ! (2) No warranties, express or implied, are made for this program.    c
1260  !                                                                      c
1261  !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1262  !
1263contains
1264
1265  SUBROUTINE gaussian_seed(i)
1266    implicit none
1267    integer i
1268    alex_ISEED=i
1269  end SUBROUTINE gaussian_seed
1270
1271  SUBROUTINE GRNF(X,cut)
1272    implicit none
1273    real(dp) r1,r2,x,cut
1274
1275
12761   R1 = -LOG(1.0_dp-RANF())
1277    R2 = 2.0_dp*PI*RANF()
1278    R1 = SQRT(2.0_dp*R1)
1279    X  = R1*COS(R2)
1280    if(abs(x)>cut) goto 1
1281    ! Y  = R1*SIN(R2)
1282    RETURN
1283  END SUBROUTINE GRNF
1284  !
1285  real(dp) FUNCTION RANF()
1286    implicit none
1287    integer ia,ic,iq,ir,ih,il,it
1288    DATA IA/16807/,IC/2147483647/,IQ/127773/,IR/2836/
1289    IH = alex_ISEED/IQ
1290    IL = MOD(alex_ISEED,IQ)
1291    IT = IA*IL-IR*IH
1292    IF(IT.GT.0) THEN
1293       alex_ISEED = IT
1294    ELSE
1295       alex_ISEED = IC+IT
1296    END IF
1297    RANF = alex_ISEED/FLOAT(IC)
1298    RETURN
1299  END FUNCTION RANF
1300
1301end module gauss_dis
1302
1303integer function mypause(i)
1304  use precision_constants
1305  implicit none
1306  !
1307  ! Replaces obsolescent feature pause
1308  !
1309  integer i
1310  !
1311  !  write (*,'(A,i6)') ' PAUSE: ',i
1312  w_p=1
1313  w_p=(/i/); w_p%fi='(1x,i4)'
1314  w_p%nc=1
1315  w_p%c(1)=' ipause=mypause(0)  ';w_p%fc='((A8,1x))'
1316
1317  ! call ! WRITE_I
1318  ! read(*,*) I
1319  mypause=i
1320  ! mypause=sqrt(dble(-i))
1321end function mypause
1322
1323integer function mypauses(i,string)
1324  use precision_constants
1325  implicit none
1326  !
1327  ! Replaces obsolescent feature pause
1328  !
1329  integer i,l,n
1330  character(*) string
1331  !
1332  !  write (*,'(A,i6)') ' PAUSE: ',i
1333  w_p=1
1334  w_p=(/i/); w_p%fi='(1x,i4)'
1335  w_p%nc=2
1336  l=len(string)
1337  call get_ncar(n)
1338  if(l>n) l=120
1339  w_p%c(1)=string(1:l)
1340  w_p%c(2)=' ipause=mypause(0)  ';w_p%fc='((A120,1x,/,a8,1x,))'
1341
1342  ! call ! WRITE_I
1343  !  read(*,*) I
1344  mypauses=i
1345  !  mypauses=sqrt(dble(-i))
1346end function mypauses
Note: See TracBrowser for help on using the repository browser.