IMPLICIT REAL*8(a-h,o-z)
REAL*8 Ki(500,100),molefrac(0:100)
PARAMETER(Npress= 16.Ncomp= l0)
DIMENSION alpha(500), beta(500), coefficien(0:100),
@ Pi(500), Ti(500), tarray(2), xalpha(500)
*
* Data Input
*
* The number of components (Ncomp) and the number of data sets
* to be run (Npress) are specified as PARAMETERs'
*
* Open and Rewind Input and Output Files
*
OPEN(unit= 1 ,file='indata'status= 'old')
OPEN(unit=-7,file='table',status='unknown')
OPEN(unit=8,file= 'plot',status='unknown')
REWIND(unit= 1)
REWIND(unit=-7)
REWIND(unit=-8)
read(1,*) (molefrac(i), i = 1, Ncomp)
do 1000 j = 1, Npress
read(1,*) Pi(j), Ti(j), beta(j)
read(l,*) (Kij,i), i = 1, Ncomp)
xalpha(j) = 1.d0 - beta(j)
1000 continue
*
* Choose between single or multiple runs
*
write(6,*) 'Evaluate one data set? enter 1'
write(6,*) 'Evaluate one data set? enter 2'
read(5,*') numsets
if(numsets .EQ. 1) then
write(6,*) 'Enter number of data set for this run'
read(5,*) j
go to 2 100
end if
do 2000 j = 1, Npress
2100 write(7,*) ' '
write(7,*) '
write(7,*) ' RUN 'j
write(6,*) 'J =',j
write(7,2500) Pi(j),Ti(j),beta(j)
2500 format('PIrssure = ',f6.1,' psia Temperature = ,f6.1,' F
@Liquid Mole Fraction = ',f6.4)
*
* Call subroutines
*
* Calculate coefficients of polynomial
*
call APHACOEFF(Ncomp,Npressj,molefrac,Ki,coefficient)
*
* Predict the number of roots on [0,1] by Founier-Budan
theorem
*
call BUDAN(jNcomp,coefficientnumroot)
*
* Solve for the roots by Newton-Raphson method
*
call ALPHAROOT(j,Ncomp,coefficient,xalpha,numroot,alpha)
*
* Generate various plots (EDIT the file to remove comments
for specific
* options)
*
call ALPHAPLOT(Ncompj,molefrac,alpha,coefficient,Ki)
* ALPHAROOT has internal output section to compile a table
* listing statistics on the determination of alpha
2000 continue
*****************************************************************************
* Produce this format to plot data points as dots:
* (PLOTFAT=20)
*
* 2
* x(1) y(1)
* x(1) y(l)
* 2
* x(2) y(2)
* x(2) y(2)
* etc.
do 3000 j = 1, Npress
write(8,3500) alpha(j),xalphaj),alpha(j),xalpha(j)
3500 format("2 'j,el 6.9,10x,e 16.9j,e 16.9,10x,e 16.9)
3000 continue
CLOSE(unit=1)
CLOSE(unit=7)
CLOSE(unit=8)
stop
end
*****************************************************************************
You don't say why you are posting this.
It is not F77 code but F90 code (e.g. the use of KIND and
status='unknown' and varibe names longer than 6 characters and
comments defiened as an asterisk in column 1 instead of "C").
And no compiler I have used recently will accept a negative integer as
a unit number (perhaps this is permitted in some compilers of this
century). It may very well be be a compiler-specific trick to pass a
parameter to the open statement. I see the same positive unit numbers
are used later.
What's the problem? (Unless you cannot compile it with and F77
compiler; in which case try F90.)
Therence,
There is no question in this post.
It's a regular spam from Mr. Musatov, I assume.
Thanks, how were you able to distinguish?
I'm puzzled by that one also that either.
> It is not F77 code but F90 code (e.g. the use of KIND and
> status='unknown' and varibe names longer than 6 characters and
> comments defiened as an asterisk in column 1 instead of "C").
I see no usage of kind here. The real*8 syntax is not a kind parameter.
It is neither Fortran 77 nor Fortran 90, but instead is a common
nonstandard feature. Although the feature is common in both f77 and f90
compilers, it is a feature associated more with f66 and f77 than with
f90 as it predates kind parameters.
Status='unknown' and comments using an asterisk in column 1 are both
perfectly standard f77. If anything, the asterisk comment style is more
an f77 one than an f90 one, insomuch as it applies only to fixed source
form, which is even obsolescent as of f95.
Variable names longer than 6 characters are indeed a feature new to the
standard as of f90, but they are one of the most common extensions out
there in f77 compilers. If one is going to be quite that strict about
the matter, one might also note the use of lower case as being a
similarly nonstandard feature in f77, along with the use of the @
character, which isn't in the f77 character set.
If I were going to get real picky about "legalisms", I'd probably note
that the "in Perpetuity" part of "All Rights Reserved in Perpetuity"
doesn't follow standards either, but that's a very different set of
standards from the Fortran ones. :-)
> And no compiler I have used recently will accept a negative integer as
> a unit number
There are several things here that no compiler of any vintage would
accept. I suppose it is vaguely possible that the consequent compilation
errors are supposed to be the question. Most of them are easy to miss in
a quick skim. I don't know whether this was posted here before throwing
it at a compiler or perhaps this was manually transcribed here or done
with OCR. For example, I see a missing parens in the statement before
3500. And the format in 3500 just looks garbled in a wy I can't figure
out. I suspect the -7 unit number to be some kind of simillar
transcription error or typo, as positive 7 seems to be used later.
--
Richard Maine | Good judgment comes from experience;
email: last name at domain . net | experience comes from bad judgment.
domain: summertriangle | -- Mark Twain
>You don't say why you are posting this.
>It is not F77 code but F90 code (e.g. the use of KIND and
>status='unknown' and varibe names longer than 6 characters and
>comments defiened as an asterisk in column 1 instead of "C").
Looks like F77 to me.
* was common in F77.
>And no compiler I have used recently will accept a negative integer as
>a unit number (perhaps this is permitted in some compilers of this
>century). It may very well be be a compiler-specific trick to pass a
>parameter to the open statement. I see the same positive unit numbers
>are used later.
Looks like a typo.
I have three F77 compilers (Burroughs, IBM and Microsoft).
Fisrt problem was the cut-and paste operation via notepad inserted
newline symbols in position 73. When these were cleared up the
program still failed for what I also identified as did Richard, the
incomplete copying, possible by OCR, from the origin of the source
posted. Then there was the use of @ as a continuation symbol, (and in
the wrong column) problems with zero indices to start a matrix and so
on, mis-matched apostophy types, missing index brackets and more.
I missed rhe valid asterisk use, as they all got translated to a
windows-coded special symbol and not the hex 2A ascii symbol during
the cut-and-paste and subsequent editing, and I was ignoring these
errors, and globally fixed them, but noting yet another "error" type.
Still, there were other problems that were just plain rejected.
But by that point it was clear the code was no way near being of F77
origin.
There may well be extended versions by some compiler vendors that
allow some non-standard features to get by, but I doubt any accepts
all of the usage claimed to be F77 code, even after removing the
obvious errors. I admit I missed the zero index clincher because the
error message assumed some other problem ("upper bound lower that
lower bound"); and I didn't go back to read the original.
> >And no compiler I have used recently will accept a negative integer as
> >a unit number (perhaps this is permitted in some compilers of this
> >century).
I have used several f77 compilers that preconnected specific
negative unit numbers to stdin, stdout, and stderr. This was fairly
common starting in the 1980's. I don't think you could use negative
unit numbers in general (e.g. in open statements), but these
preconnected units were treated differently. Also, it was not very
portable, one compiler might use 0, -1, and -2 while another might
use -1, -2, and -3.
I'm assuming here that the 1980's was last century. If by the last
century you mean before 1909, then that is a different matter. :-)
> It may very well be be a compiler-specific trick to pass a
> >parameter to the open statement. I see the same positive unit numbers
> >are used later.
>
> Looks like a typo.
In this case that is probably right, a typo or an OCR error.
$.02 -Ron Shepard
The common numbers were 1, 2, and 3 and 5, and 6.
| I'm assuming here that the 1980's was last century. If by the last
| century you mean before 1909, then that is a different matter. :-)
In omitting the poster, you have confused who sent what.
I did not say that; it was Terence (whose POST ID I have restored).
Are you still with me?
Sorry, but I am stuck with Google access and the quirks of waht
happens if you cut anything, or answer before a opposed to after the
original positied lines.
Hey-ho.
Anyway, to me, the "last century" is 1900 plus any two digits between
00 and 99.
If you need my rememberances of pre 1900, I doh't recall many (>=0).
musatov claims his motive is to help children. Yet he
keeps reminding everyone he is trademarking and copyrighting everything he does -- while showing no
qualms in repeatedly pasting here without attributing sources .
It is dificult to believe, to say the least, that someone who claims to be interested only in doing charitable work shows such zeal in copyrighting and trademarking everything he does. Arguably, this character is also using others' posts --without permission --for his personal research.
While musatov claims to be interested uniquely in
helping children, he explicitly describes himself as
a 'web enterpreneur' . This seems troublesome.
This character musatov is so far gone that he has accused those he claims are lying, of being Satan -- I am not kidding here. Yet, musatov has had no trouble whatsoever in falsely claiming that his search engine had surpassed google in several countries.
This self-appointed champion of honesty also appended 'AP' at the beginning of his article , to make it seem as if the associated press had written said article. Only when several called him up on it did he
admit this was false. Maybe 'Satan' possesed him while he was lying about it.
I suspect any content of your post will be used to produce profit for musatov. musatov is zealous about defending _his_ intellectual property. But his zeal does not extend to respecting others' property.
<snip>
> Anyway, to me, the "last century" is 1900 plus any two digits
> between 00 and 99.
To you, does 2 + 2 equal 5? This century is the 21st Century, so last
century was the 20th Century. The first Century ran from 0001AD to
0100AD (there was no year 0), so the 20th Century ran from 1901AD to
2000AD.
<snip>
--
Richard Heathfield <http://www.cpax.org.uk>
Email: -http://www. +rjh@
"Usenet is a strange place" - dmr 29 July 1999
Sig line vacant - apply within
Sounds fine to me.
What decade is it currently.
By what year, according to Kennedy's speech, was the moon landing
supposed to occur?
-- glen
> In comp.lang.fortran Richard Heathfield <r...@see.sig.invalid> wrote:
>
>> To you, does 2 + 2 equal 5? This century is the 21st Century, so
>> last century was the 20th Century. The first Century ran from
>> 0001AD to 0100AD (there was no year 0), so the 20th Century ran
>> from 1901AD to 2000AD.
>
> Sounds fine to me.
>
> What decade is it currently.
The 201st AD.
> By what year, according to Kennedy's speech, was the moon landing
> supposed to occur?
The decade in question was the 197th, 1961-1970.
Why did NASA decide that the moon landing had to be in 1969?
It seems that they didn't believe that they had one more year.
-- glen
JFK set 1970 as a limit, not a target. As you rightly asked, "by what
year", not "in what year".
You only have to go back 9 posts to see that you were in fact the poster of that sentence.
Quote:
________________________________
From: "Terence" <tbwr...@cantv.net>
Newsgroups: comp.lang.fortran,comp.programming,sci.math,comp.infosystems,rec.org.mensa
Sent: Saturday, 10 October 2009 6:09 PM
Subject: Re: * Program ALPHATEST (FORTRAN 77)
You don't say why you are posting this.
It is not F77 code but F90 code (e.g. the use of KIND and
status='unknown' and varibe names longer than 6 characters and
comments defiened as an asterisk in column 1 instead of "C")..
> JFK set 1970 as a limit, not a target. As you rightly asked, "by what
> year", not "in what year".
Maybe, but as I understand it NASA considered the goal to be 1969.
July gave them some more chances before December 1969.
Note that unlike centuries (20th, 21st) decades are rarely
if ever described that way, but usually the 60's, or the 70's.
Technically NASA had until the end of December 1970, but the
public wasn't going to give them that long, and so they couldn't
schedule it that way.
-- glen
Ron Shepard poted the message responded to by Robin.
...
>| century you mean before 1909, then that is a different matter. :-)
Then Robin responded to this, WHETHER OR NOT HE MEANT TO, BY :-
>In omitting the poster, you have confused who sent what.
>I did not say that; it was Terence (whose POST ID I have restored).
But Terence's POST ID was in ever message he (I) posted!!
Which is why I responded to Robin that it WASn't my posting but Ron's
response to my posting.
And OK, a centry runs xx*100 +01 to (XX+1)*100, and not xx00-xx99.
I still think I get 99 out of 100.
The reference was STILL "last century" as I said, which I lived over
2/3 of, plus all I clock up in this one, still (of course) programming
in F77. (The WORD as given). :o)>
Even the concept "AD" is pretty well considered to be off by about 4
years anyway, just as the census in Bethleham was in June -3 as per
Roman records, and not December 25 0001.
abcdefghijklm o
r
t
a
b
l
e
c
o
d
e
c o
m p r
e s
s
abcdefghi o
abcdefghijklmn
abcdefghijklmnopq
abcdefghijklmnopqrstu
abcdefghijklmno
abcd
abcde
abcdefghijklmnopqr
a
abcdefghijklmnopqrst
abcd
abcde
abcdefghijklm
abcdefghijklmno
abcdefghijklmn
abcdefghijklmnonpqrs
abcdefghijklmnopqrst
abcdefghijklmnopqr
a
abcdefghijklmn
abcd
abcdefghijklmnopqrstu
abcdefghijklmnop====
1. Q.E.D. - Wikipedia, the free encyclopedia Q.E.D. is an abbreviation
of the Latin phrase quod erat demonstrandum, which literally means
"which was to be demonstrated". The phrase is written in
its ...Etymology and early use - Modern philosophy - QEF
en.wikipedia.org/wiki/Q.E.D. - 2. Urban Dictionary: quod erat
demonstrandum Originally Latin meaning "quod erat demonstrandum" or
"which was to be shown or proven", now used mainly by physics students
to insult someone when ...www.urbandictionary.com/define.php?...quod
+erat+demonstrandum
3. Quod Erat Demonstrandum... computer graphics - gallery 10 Apr
2009 ... Mysterious story teller with visual imagery (2d and 3d).
Inviting, exciting, and visually imaginative works: Surreal, Fantasy
and SciFi ...qed.inet24.pl/indexa.htm
4. Quod erat demonstrandum - Idioms - by the Free
Dictionary ...Definition of Quod erat demonstrandum in the Idioms
Dictionary. Quod erat demonstrandum phrase. What does Quod erat
demonstrandum expression mean?idioms.thefreedictionary.com/Quod+erat
+demonstrandum
5. quod erat demonstrandum - What is a(n) quod erat
demonstrandum ...quod erat demonstrandum - What is a(n) quod erat
demonstrandum from The Oxford Dictionary of Phrase and Fable at
Encyclopedia.com.www.encyclopedia.com/.../1O214-
quoderatdemonstrandum.html
6. LOREM IPSUM - Quod erat demonstrandum Aliquam erat volutpat. In at
mi ac nisi mollis posuere. ... Praesent venenatis, erat nec vulputate
molestie, elit felis dictum leo, non faucibus risus
erat ...www.loremipsum.ro/
7. Quod Erat Demonstrandum (Q.E.D.) - Definition of Q.E.D. Definition:
Quod erat demonstrandum (Q.E.D.) is the Latin for that which was to be
demonstrated. Q.E.D. is used in mathematical proofs to show that what
was ...ancienthistory.about.com/od/qterms/g/quoderatdemonst.htm
8. Dictionary of Difficult Words - quod erat demonstrandumquod erat
demonstrandum. 'which was to be demonstrated' (abbr. Q.E.D.). quod
erat faciendum, 'which was to be done' (abbr. Q.E.F.). © RM
2009. ...www.tiscali.co.uk/reference/dictionaries/.../d0010815.html
9. quod erat demonstrandum - Wiktionary 3 Aug 2009 ...
quod (nominative neuter singular of quī) + erat (third-person singular
imperfect active indicative of sum) + dēmonstrandum
(nominative ...en.wiktionary.org/wiki/quod_erat_demonstrandum
10. Quoderat demonstrandum? The mystery of experimental validation
of ...Quod erat demonstrandum? The mystery of experimental validation
of apparently erroneous computational analyses of protein sequences.
Iyer LM, Aravind L, ...www.ncbi.nlm.nih.gov/pubmed/11790254 -
Similarby LM Iyer - 2001 - Cited by 42 - Related articles - All 20
versions
Searches related to: quod erat demonstrandum quod erat demonstrandum
meaning qed quod erat demonstrandum quod erat demonstrandum
pronunciation quod erat demonstratum quad erat
demonstratum  1  2  3  4  5  6  7  8  9  10  Next  Personalized based
on your web history. Cop::y(C)MeAmI wrote: 2009. M. Michael
Musatov. http://www.meami.org space:
'Search for the for: People!'. Lease
I would like to thank Alexandra Scott and all of my other angels who
have helped me. God blesses me with you and it is time we bless the
children dying of cancer and award the Millennium Prize for The P
VERSUS NP problem for donation based on this work (in part of a much
larger body of work) with all one-million of each dollars will be
donated to http://alexslemonade.org
> > [code elided]
> > > You don't say why you are posting this.
> >
For this: P y\NP> > I'm puzzled by that one also that either.
>D
Silencer
M
R
you psycho loser.
Nabble - gcc - fortran - [PATCH, gfortran testsuite]: Do not ...
20 posts - 7 authors - Last post: Aug 20
current gcc/gfortran with iso_c_interop for the i386/x64 platforms.
I ....
if IEEE_SET_UNDERFLOW is written according to iso_c_interop ? ...
www.nabble.com/-PATCH,-gfortran-testsuite-:-Do-not-generate-denormals-in-gfortran.dg-boz_9.f90-td25082961.html
Tim Prince - Re: [PATCH, gfortran testsuite]: Do not load ...
Aug 24, 2009 ...
straightforward in current gcc/gfortran with iso_c_interop for the
i386/x64 platforms.
I don't have any other targets available to test ...
gcc.gnu.org/ml/fortran/2009-08/msg00356.html
J3/04-230 Date: 02 February 2004 To: J3 From: Aleksandar Donev ......
which one should be able to write an equivalent module: MODULE GL USE
ISO_C_INTEROP TYPEDEF :: GLenum=>INTEGER(C_INT), & GLboolean=>CHARACTER
(C_CHAR), . ...
www.j3-fortran.org/doc/year/04/04-230.txt
J3/04-230
Date: 02 February 2004
To: J3
From: Aleksandar Donev
Subject: TYPEDEF facility
Title: TYPEDEF facility
Submitted by: J3
Status: For Consideration
References: J3/03-256
Basic Functionality:
Bring back the previously removed TYPEALIAS facility as a TYPEDEF
facility to be used in C Interop. This facility would allow one to
give
alias names to other types. This facility may be restricted to
Index: gfortran.dg/boz_9.f90
===================================================================
--- gfortran.dg/boz_9.f90 (revision 150988)
+++ gfortran.dg/boz_9.f90 (working copy)
@@ -1,6 +1,5 @@
! { dg-do run }
! { dg-options "-fno-range-check" }
-! { dg-options "-fno-range-check -mieee" { target alpha*-*-* } }
!
! PR fortran/34342
!
@@ -10,40 +9,40 @@
implicit none
real,parameter :: r2c = real(int(z'3333'))
-real,parameter :: rc = real(z'3333')
+real,parameter :: rc = real(z'50CB9F09')
double precision,parameter :: dc = dble(Z'3FD34413509F79FF')
-complex,parameter :: z1c = cmplx(b'10101',-4.0)
-complex,parameter :: z2c = cmplx(5.0, o'01245')
+complex,parameter :: z1c = cmplx
(b'11000001010001101101110110000011', 3.049426e-10)
+complex,parameter :: z2c = cmplx(4.160326e16,
o'6503667306')
real :: r2 = real(int(z'3333'))
-real :: r = real(z'3333')
+real :: r = real(z'50CB9F09')
double precision :: d = dble(Z'3FD34413509F79FF')
-complex :: z1 = cmplx(b'10101',-4.0)
-complex :: z2 = cmplx(5.0, o'01245')
+complex :: z1 = cmplx(b'11000001010001101101110110000011',
3.049426e-10)
+complex :: z2 = cmplx(4.160326e16, o'6503667306')
if (r2c /= 13107.0) call abort()
-if (rc /= 1.83668190E-41) call abort()
+if (rc /= 2.732958e10) call abort()
if (dc /= 0.30102999566398120d0) call abort()
-if (real(z1c) /= 2.94272678E-44 .or. aimag(z1c) /= -4.0) call abort
()
-if (real(z2c) /= 5.0 .or. aimag(z2c) /= 9.48679060E-43) call abort()
+if (real(z1c) /= -1.242908e1 .or. aimag(z1c) /= 3.049426e-10) call
abort()
+if (real(z2c) /= 4.160326e16 .or. aimag(z2c) /= 5.343285e-7) call
abort()
if (r2 /= 13107.0) call abort()
-if (r /= 1.83668190E-41) call abort()
+if (r /= 2.732958e10) call abort()
if (d /= 0.30102999566398120d0) call abort()
-if (real(z1) /= 2.94272678E-44 .or. aimag(z1) /= -4.0) call abort()
-if (real(z2) /= 5.0 .or. aimag(z2) /= 9.48679060E-43) call abort()
+if (real(z1) /= -1.242908e1 .or. aimag(z1) /= 3.049426e-10) call abort
()
+if (real(z2) /= 4.160326e16 .or. aimag(z2) /= 5.343285e-7) call abort
()
r2 = dble(int(z'3333'))
-r = real(z'3333')
+r = real(z'50CB9F09')
d = dble(Z'3FD34413509F79FF')
-z1 = cmplx(b'10101',-4.0)
-z2 = cmplx(5.0, o'01245')
+z1 = cmplx(b'11000001010001101101110110000011', 3.049426e-10)
+z2 = cmplx(4.160326e16, o'6503667306')
-if (r2 /= 13107.0) call abort()
-if (r /= 1.83668190E-41) call abort()
+if (r2 /= 13107d0) call abort()
+if (r /= 2.732958e10) call abort()
if (d /= 0.30102999566398120d0) call abort()
-if (real(z1) /= 2.94272678E-44 .or. aimag(z1) /= -4.0) call abort()
-if (real(z2) /= 5.0 .or. aimag(z2) /= 9.48679060E-43) call abort()
+if (real(z1) /= -1.242908e1 .or. aimag(z1) /= 3.049426e-10) call abort
()
+if (real(z2) /= 4.160326e16 .or. aimag(z2) /= 5.343285e-7) call abort
()
call test4()
call test8()
@@ -52,34 +51,34 @@
subroutine test4
real,parameter :: r2c = real(int(z'3333', kind=4),
kind=4)
-real,parameter :: rc = real(z'3333', kind=4)
-complex,parameter :: z1c = cmplx(b'10101',-4.0, kind=4)
-complex,parameter :: z2c = cmplx(5.0, o'01245', kind=4)
+real,parameter :: rc = real(z'50CB9F09', kind=4)
+complex,parameter :: z1c = cmplx
(b'11000001010001101101110110000011', 3.049426e-10, kind=4)
+complex,parameter :: z2c = cmplx(4.160326e16, o'6503667306',
kind=4)
real :: r2 = real(int(z'3333', kind=4), kind=4)
-real :: r = real(z'3333', kind=4)
-complex :: z1 = cmplx(b'10101',-4.0, kind=4)
-complex :: z2 = cmplx(5.0, o'01245', kind=4)
+real :: r = real(z'50CB9F09', kind=4)
+complex :: z1 = cmplx(b'11000001010001101101110110000011',
3.049426e-10, kind=4)
+complex :: z2 = cmplx(4.160326e16, o'6503667306', kind=4)
if (r2c /= 13107.0) call abort()
-if (rc /= 1.83668190E-41) call abort()
-if (real(z1c) /= 2.94272678E-44 .or. aimag(z1c) /= -4.0) call abort
()
-if (real(z2c) /= 5.0 .or. aimag(z2c) /= 9.48679060E-43) call abort()
+if (rc /= 2.732958e10) call abort()
+if (real(z1) /= -1.242908e1 .or. aimag(z1) /= 3.049426e-10) call abort
()
+if (real(z2) /= 4.160326e16 .or. aimag(z2) /= 5.343285e-7) call abort
()
if (r2 /= 13107.0) call abort()
-if (r /= 1.83668190E-41) call abort()
-if (real(z1) /= 2.94272678E-44 .or. aimag(z1) /= -4.0) call abort()
-if (real(z2) /= 5.0 .or. aimag(z2) /= 9.48679060E-43) call abort()
+if (r /= 2.732958e10) call abort()
+if (real(z1) /= -1.242908e1 .or. aimag(z1) /= 3.049426e-10) call abort
()
+if (real(z2) /= 4.160326e16 .or. aimag(z2) /= 5.343285e-7) call abort
()
r2 = real(int(z'3333'), kind=4)
-r = real(z'3333', kind=4)
-z1 = cmplx(b'10101',-4.0, kind=4)
-z2 = cmplx(5.0, o'01245', kind=4)
+r = real(z'50CB9F09', kind=4)
+z1 = cmplx(b'11000001010001101101110110000011', 3.049426e-10,
kind=4)
+z2 = cmplx(4.160326e16, o'6503667306', kind=4)
if (r2 /= 13107.0) call abort()
-if (r /= 1.83668190E-41) call abort()
-if (real(z1) /= 2.94272678E-44 .or. aimag(z1) /= -4.0) call abort()
-if (real(z2) /= 5.0 .or. aimag(z2) /= 9.48679060E-43) call abort()
+if (r /= 2.732958e10) call abort()
+if (real(z1) /= -1.242908e1 .or. aimag(z1) /= 3.049426e-10) call abort
()
+if (real(z2) /= 4.160326e16 .or. aimag(z2) /= 5.343285e-7) call abort
()
end subroutine test4
interoperable types only if so desired. A proper and better designed
Fortran TYPEALIAS facility may then be developed separately, not
compromising the very urgent C Interop void created by the removal of
TYPEALIAS.
Rationale:
C has a typedef facility which allows one to give aliases to existing
C
types. These are used very widely, and therefore we need to provide a
way to easily interface with libraries that use them. Header files are
central to "portability" in C, and they usually mainly consist of
preprocessing directives like #define's, function prototypes and
typedef's. The first seem outside the scope of Interop and can be
emulated with PARAMETERs, the second are
handled already, and now we also need to be able to handle typedef's.
The goal is to be able to write a Fortran module that emulates a
header
file provided by the vendor of a library, and that someday it may
even
be possible to do an automatic .h->.f90 conversion (at least an
initial
stage thereof). Here is a short piece of the include file GL/gl.h
which
is the central include file for the OpenGL library:
...
typedef unsigned int GLenum;
typedef unsigned char GLboolean;
typedef unsigned int GLbitfield;
typedef void GLvoid;
...
#define GL_BYTE 0x1400
#define GL_UNSIGNED_BYTE 0x1401
...
GLAPI void GLAPIENTRY glTexCoord1dv( const GLdouble *v );
GLAPI void GLAPIENTRY glTexCoord1fv( const GLfloat *v );
for which one should be able to write an equivalent module:
MODULE GL
USE ISO_C_INTEROP
TYPEDEF :: GLenum=>INTEGER(C_INT), &
GLboolean=>CHARACTER(C_CHAR), ...
TYPE(GLenum), PARAMETER :: GL_BYTE=...
INTERFACE
SUBROUTINE glTexCoord1dv(v)
TYPE(GLdouble), DIMENSION(*), INTENT(IN) :: v
END SUBROUTINE
...
END INTERFACE
END MODULE GL
As another example, consider writing a Fortran interface to the MPICH
implementation of MPI. Typically, this will be a module that contains
various constants, type parameters, and interfaces. MPI uses many
typealiases, which are needed when writing interfaces. For example,
MPI_Datatype is typically an alias for int. But one cannot assume
this,
nor that it is indeed an integer. One cannot get away with our untyped
C_PTR, since arguments of type MPI_Datatype are passed by value, not
by
reference. It is necessary for any kind of portability that one be
able
to write:
TYPEDEF :: MPI_Datatype=>INTEGER(KIND=C_INT)
in the module for the interface to MPICH.
Estimated Impact:
The complication is mostly syntactic. The same issues we had with
TYPEALIAS remain. If we restrict this facility only to interoperable
types some of the problems may go away (for example, no more
parameterized types).
Detailed Specification:
Same as previous TYPEALIAS, but called TYPEDEF, possibly restricted
only to interoperable types.
History:
Support a cure for childhood cancer: Alex's Lemonade
©2009 MeAmI.org "Search for the People!"
V = 180i + 450j/*
Covert Tunnelling in ICMP 0x00 ECHO REPLY messages
Many thanks to FuSyS and Richard Stevens ^_^
Dark Schneider X1999
*/
#include <winsock2.h>
#include <ws2tcpip.h>
#include <stdio.h>
#define ICMP_ECHOREPLY 0
#define ICMP_ECHO 8
// definizione di alcune costanti
#define IP_HDR 20
#define ICMP_HDR 8
#define ICMP_MINLEN 8
#define MAXMESG 4096
#define MAXPACKET 5004
#define LAST 1
#define REPLY 1
#define ECHO_TAG 0xF001
#define ECHO_LAST 0xF002
// Strutture e Variabili
// Lancio un doveroso Porko D*io liberatorio... dopo ore ho trovato
come fare
// a togliermi dalle palle la fottuta icmp.dll (winsock maledette)
// IP Header
struct ip
{
unsigned char Hlen:4;
unsigned char Version:4;
unsigned char Tos;
unsigned short LungTot;
unsigned short Id;
unsigned short Flags;
unsigned char Ttl;
unsigned char Proto;
unsigned short Checksum;
unsigned int SourceIP;
unsigned int DestIP;
};
// ICMP Header
struct icmp {
BYTE Type;
BYTE Code;
USHORT CheckSum;
USHORT Id;
USHORT Seq;
ULONG Dati;
};
SOCKET sockfd;
u_int icmp_init =1;
struct sockaddr_in clisrc;
// Funzione di checksum
USHORT checksum(USHORT *buffer, int size)
{
unsigned long cksum=0;
while(size >1)
{
cksum+=*buffer++;
size -=sizeof(USHORT);
}
if(size )
{
cksum += *(UCHAR*)buffer;
}
cksum = (cksum >> 16) + (cksum & 0xffff);
cksum += (cksum >>16);
return (USHORT)(~cksum);
}
// Reimplemento bcopy e bzero... Ma perche' cavolo windows non le
// rende disponibili?
void bzero(char *pnt, int dim_pnt )
{
memset((char *)&pnt, 0, dim_pnt);
};
void bcopy(char *src, char *dest, int dim_src)
{
memmove((char *)&dest, (char *)&src, dim_src);
};
// Micro$oft Sucks
// Funzioni di gestione dei pacchetti ICMP
// Fankulo a quegli stronzi maledetti che si sono inventati la
icmp.dll
// Brutti bastardi pezzi di merda, la compatibilita' ve la siete
ficcata su
// per il culo?
// Micro$oft Sucks
void ICMP_init(void)
{
if(icmp_init)
{
if((sockfd = socket(AF_INET, SOCK_RAW, IPPROTO_ICMP)) ==
INVALID_SOCKET)
{
fprintf(stderr, "impossibile creare raw ICMP socket");
exit(0);
}
}
icmp_init = 0;
};
void ICMP_reset(void)
{
closesocket(sockfd);
icmp_init = 1;
};
int ICMP_send(char *send_mesg, size_t mesglen, ULONG dest_ip, int
echo, int last)
{
int sparato;
struct tunnel
{
struct icmp icmp;
UCHAR data[MAXMESG];
} icmp_pk;
int icmplen = sizeof(struct icmp);
int pack_dim;
struct sockaddr_in dest;
int destlen;
if(mesglen > MAXMESG) return(-1);
if(icmp_init) ICMP_init();
destlen = sizeof(dest);
bzero((char *)&dest, destlen);
dest.sin_family = AF_INET;
dest.sin_addr.s_addr = dest_ip;
pack_dim = mesglen + sizeof(struct icmp);
memset(&icmp_pk, 0, pack_dim);
icmp_pk.icmp.Type = ICMP_ECHOREPLY;
bcopy(send_mesg, (char *)&icmp_pk.icmp.Dati, mesglen);
icmp_pk.icmp.CheckSum = checksum((USHORT *) &icmp_pk.icmp, (sizeof
(struct icmp) + mesglen));
if(echo) icmp_pk.icmp.Seq = ECHO_TAG;
if(last) icmp_pk.icmp.Seq = ECHO_LAST;
if(sparato = sendto(sockfd, (char *)&icmp_pk, pack_dim, 0, (struct
sockaddr *)&dest, destlen) < 0)
{
perror("RAW ICMP SendTo: ");
return(-1);
}
else if(sparato != pack_dim)
{
perror("dimensioni pacchetto IP errate: ");
return(-1);
}
return(sparato);
};
int ICMP_recv(char *recv_mesg, size_t mesglen, int echo)
{
struct recv
{
struct ip ip;
struct icmp icmp;
char data[MAXMESG];
} rcv_pk;
int pack_dim;
int accolto;
int iphdrlen;
int clilen = sizeof(clisrc);
if(icmp_init) ICMP_init();
while(1)
{
pack_dim = mesglen + sizeof(struct ip) + sizeof(struct icmp);
memset(&rcv_pk, 0, pack_dim);
if((accolto = recvfrom(sockfd, (char *)&rcv_pk, pack_dim, 0, (struct
sockaddr *) &clisrc, &clilen)) < 0) continue;
iphdrlen = rcv_pk.ip.Hlen << 2;
if(accolto < (iphdrlen + ICMP_MINLEN)) continue;
accolto -= iphdrlen;
if(!echo)
{
if(!rcv_pk.icmp.Id && !rcv_pk.icmp.Code && rcv_pk.icmp.Type ==
ICMP_ECHOREPLY && rcv_pk.icmp.Seq != ECHO_TAG && rcv_pk.icmp.Seq !=
ECHO_LAST) break;
}
if(echo)
{
if(!rcv_pk.icmp.Id && !rcv_pk.icmp.Code && rcv_pk.icmp.Type ==
ICMP_ECHOREPLY && (rcv_pk.icmp.Seq == ECHO_TAG || rcv_pk.icmp.Seq ==
ECHO_LAST)) break;
}
}
if(!echo)
{
accolto -= ICMP_HDR;
bcopy((char *)&rcv_pk.icmp.Dati, recv_mesg, accolto);
n = 200819. 2. 782
therefore increment t
s = a * sp + b * sq mod N. where. sp=mdp mod (p) and sq=mdq mod (q)
a/b = (rp. m+1. ± sp. m )/(rq. m+1. ± sq. m )
[a, b, s:1] presents a plot of the mod-4 prime race over the interval
[a, b]
Let it ride: http://meami.org (C) All Rights Reserved.
+
http://buildasearch.com/meami
c=Copywrite(C) 2009. [at] (c) Copyright: http://MEAMI.ORG
! Programma FINALE
! versione risc1 IBM
! versione FINALE nel caso assiale
! versione FINALE nel caso rombico (non funziona in assenza di ZFS)
C PARAMAGNETIC ENHANCEMENT IN NMRD PROFILE
C
C THIS PROGRAM REQUIRES THE INPUT FILE PAR.DAT
C 0 BEFORE A PARAMETER MEANS IT HAS TO BE ASSUMED AS CONSTANT
C 1 BEFORE A PARAMATER MEANS IT HAS TO BE CHANGED IN THE FITTING
C
C
C INTERNAL FITTING IS PERFORMED IN THE PARAMETERS:
C d , Ddiff , RK , A/h , MOLAR FRACTION , TAUM
C
C SUBROUTINES:
C FUNCZFS: SET PARAMETERS IN FITTING PROCEDUREC
C FUNCINT: SET PARAMETERS IN INTERNAL FITTING PROCEDURE
C DIAG: WRITE ENERGY MATRIX AND CALCULATE EXPECTATION VALUES
C POWELL: FITTING PROCEDURE
C CONNECTED SUBROUTINES:
C LINMIN
C F1DIM
C MNBRAK
C MRENT
C XISTEP
C POWELLINT: INTERNAL FITTING PROCEDURE
C F01BCF....X04BAF: CALCULATE EIGENVALUES AND EIGENFUNCTIONS
C GAUINT: PERFORME INTEGRATION ON ANGLES
C TUNO: PERFORME CALCULATION OF T1M
C TDUE: PERFORME CALCULATION OF T2M
C TUNOISO: PERFORME CALCULATION OF T1M IN AXIAL CASE
C
C THETA AND PHI ARE THE POLAR ANGLES DEFINING THE WATER PROTON
DIRECTION
C IN THE MOLECULAR FRAME
C
C BETA AND GAMMA ARE THE EULER ANGLES DEFINING THE MOLECULAR FRAME
WITH
C RESPECT TO THE HTTP://MEAMI.ORG/ FRAME
IMPLICIT REAL*8(A-H,O-Z)
PARAMETER(PI2 = 6.2831853, VL = 2.9979D+10)
COMMON /SET/SET
COMMON /RK10/ SPIN, SI
COMMON /GAMMAH/ GAMMAI
COMMON /B31/ TPUNO(500)
COMMON /B32/ PP(500),Z(500)
COMMON /TAUDELTA/ TAUDELTA
COMMON /B4/ NVMEM,NPT(10),NPTOT
COMMON /WATER/ ACQ
COMMON /T1T2/ IREL
COMMON /INDEX/ INDEX
COMMON /ALFASTEP/ ALFASTEP
COMMON /STAMPA/INDEXSTAMPA
COMMON /TOT/ DPARATOT,EPARATOT,APARTOT,APERTOT,APERTOT2,ACONIND
COMMON/INDA/ INDDPARA(10),INDAPAR(10),INDAPER(10),INDEPARA(10)
& ,INDED(10),INDAMOLFRA(10),INDS4(10),INDAPER2(10)
COMMON /C1M/DM(10),DDM(10),CONCM(10)
COMMON /IPERFM/AZM(10),AYM(10),AXM(10),THETAM(10),RKM(10),
& TAUCM(10),DPARAM(10),EPARAM(10),PHIM(10),S4M(10),
& GXM(10),GYM(10),GZM(10)
COMMON /TAU1M/ TAUS0M(10,10),TAURM(10,10),TAUVM(10,10),
& TAUMM(10,10)
COMMON /MOLFRAZM/ AMOLFRAM(10)
COMMON /CONTATM/ ACONTM(10)
COMMON /CICLE/ NVEST
COMMON /BPARA/B1(10),B2(10),B3(10),B4(10),B5(10),B6(10),B7(10),
& B8(10),B9(10),B10(10),B11(10),B12(10),B13(10),B14(10),B15(10),
& B16(10),B17(10),B18(10),B19(10),B20(10),B21(10)
CHARACTER*20 FILENAME
COMMON/TEMPERATURE/ TEMP(10)
COMMON /TMSTART/ TM11(10),TM21(10)
C DIMENSION=MAX NUMBER OF PARAMETERS (21)
DIMENSION P(21)
DIMENSION P1(21)
COMMON /PPAR/ P2(21)
DIMENSION XI(21,21)
INDEX=1
INDEXSTAMPA=0
C CONSTANTS READ IN FILE PAR.DAT
OPEN (1, STATUS = 'OLD', FILE = 'PARC.DAT')
OPEN (4, FILE="PAR.OUT")
C OUTPUT FILE
READ(1,'(A)')FILENAME
C NUCLEAR MOLECULAR SPIN
READ(1,*)SI
C GAMMA OF THE INVESTIGATED PARTICLE
READ(1,*)GAMMAI
C ELECTRON SPIN
READ(1,*)SPIN
C T1 OR T2 CALCULATION
READ(1,*)IREL
C LIMITS OF THE FIELD
READ(1,*)X1,X2,X3
IF(X3.EQ.1)THEN
XMIN=X1
XMAX=X2
ELSE
XMIN=LOG10(GAMMAI*X1/6.283)
XMAX=LOG10(GAMMAI*X2/6.283)
ENDIF
C NUMBER OF POINTS TO BE CALCULATED
READ(1,*)NUMPUN
IF(XMIN.EQ.XMAX)NUMPUN=1
C NUMBER OF SETS OF DATA FOR FITTING
READ(1,*)SET
IF(SET.EQ.0) SET=1
C TEMPERATURE
READ(1,*)(TEMP(K),K=1,SET)
J=1
IND=1
IND1=1
IND2=1
NV=0
NVEST=0
NPLUS=0
NPLUS2=0
C CORRELATION TIMES
READ(1,*)B1(J), (TAUS0M(J,K),K=1,2),TAUDELTA
IF(B1(J).GE.2)THEN
TS1=TAUS0M(J,1)
TS2=TAUS0M(J,2)
DO K=1,SET
TAUS0M(J,K)=TS1*EXP(TS2/TEMP(K))
END DO
IF(B1(J).EQ.3)NPLUS=NPLUS+1
ENDIF
IF(B1(J).EQ.1)THEN
P(IND)=TAUS0M(J,1)
P1(IND1)=TAUS0M(J,1)
WRITE(4,'(2X,4(E10.4,2X))') P(IND)
IND=IND+1
IND1=IND1+1
ENDIF
IF(B1(J).EQ.2)THEN
P(IND)=TS1
P1(IND1)=TS1
P(IND+1)=TS2
P1(IND1+1)=TS2
WRITE(4,'(2X,4(E10.4,2X))') TS1,TS2
IND=IND+2
IND1=IND1+2
ENDIF
READ(1,*)B2(J), (TAURM(J,K),K=1,2)
IF(B2(J).GE.2)THEN
TR1=TAURM(J,1)
TR2=TAURM(J,2)
DO K=1,SET
TAURM(J,K)=TR1*EXP(TR2/TEMP(K)) !Stokes
END DO
IF(B2(J).EQ.3) NPLUS=NPLUS+1
ENDIF
IF(B2(J).EQ.1)THEN
P(IND)=TAURM(J,1)
P1(IND1)=TAURM(J,1)
WRITE(4,'(2X,4(E10.4,2X))') P(IND)
IND=IND+1
IND1=IND1+1
ENDIF
IF(B2(J).EQ.2)THEN
P(IND)=TR1
P1(IND1)=TR1
P(IND+1)=TR2
P1(IND1+1)=TR2
WRITE(4,'(2X,4(E10.4,2X))') TR1,TR2
IND=IND+2
IND1=IND1+2
ENDIF
READ(1,*)B3(J), (TAUVM(J,K),K=1,2)
IF(B3(J).GE.2)THEN
TV1=TAUVM(J,1)
TV2=TAUVM(J,2)
DO K=1,SET
TAUVM(J,K)=TV1*EXP(TV2/TEMP(K))
END DO
IF(B3(J).EQ.3) NPLUS=NPLUS+1
ENDIF
IF(B3(J).EQ.1)THEN
P(IND)=TAUVM(J,1)
P1(IND1)=TAUVM(J,1)
WRITE(4,'(2X,4(E10.4,2X))') P(IND)
IND=IND+1
IND1=IND1+1
ENDIF
IF(B3(J).EQ.2)THEN
P(IND)=TV1
P1(IND1)=TV1
P(IND+1)=TV2
P1(IND1+1)=TV2
WRITE(4,'(2X,4(E10.4,2X))') TV1,TV2
IND=IND+2
IND1=IND1+2
ENDIF
IF (TAURM(J,1).EQ.0.AND.TAUVM(J,1).EQ.0)THEN
TAUCM(J)=TAUS0M(J,1)
IF(B1(J).EQ.1) P(IND-1)=TAUS0M(J,K)
IF(B1(J).EQ.1) P1(IND1-1)=TAUS0M(J,K)
ENDIF
C PARAMETERS OF ZERO FIELD SPLITTING
READ(1,*)B11(J), DPARAM(J)
DPARAM(J) = PI2*VL*DPARAM(J)
IF(B11(J).EQ.1)THEN
P(IND)=DPARAM(J)
P1(IND1)=DPARAM(J)
WRITE(4,'(2X,4(E10.4,2X))') P(IND)/VL/PI2
INDDPARA(J)=IND
IND=IND+1
IND1=IND1+1
ENDIF
READ(1,*)B12(J), EPARAM(J)
EPARAM(J) = PI2*VL*EPARAM(J)
IF(B12(J).EQ.1)THEN
P(IND)=EPARAM(J)
P1(IND1)=EPARAM(J)
WRITE(4,'(2X,4(E10.4,2X))') P(IND)/VL/PI2
INDEPARA(J)=IND
IND=IND+1
IND1=IND1+1
ENDIF
READ(1,*)B19(J), S4M(J)
S4M(J) = PI2*VL*S4M(J)
IF(B19(J).EQ.1)THEN
P(IND)=S4M(J)
P1(IND1)=S4M(J)
WRITE(4,'(2X,4(E10.4,2X))') P(IND)/VL/PI2
INDS4(J)=IND
IND=IND+1
IND1=IND1+1
ENDIF
C PARAMETER OF G-TENSOR
READ(1,*)B17(J), GSER
GXM(J)=GSER/2.003
IF(B17(J).EQ.1)THEN
P(IND)=GXM(J)
P1(IND1)=GXM(J)
WRITE(4,'(2X,4(E10.4,2X))') P(IND)
IND=IND+1
IND1=IND1+1
ENDIF
READ(1,*)B18(J), GSER
GYM(J)=GSER/2.003
IF(B18(J).EQ.1)THEN
P(IND)=GYM(J)
P1(IND1)=GYM(J)
WRITE(4,'(2X,4(E10.4,2X))') P(IND)
IND=IND+1
IND1=IND1+1
ENDIF
READ(1,*)B20(J), GSER
GZM(J)=GSER/2.003
IF(B20(J).EQ.1)THEN
P(IND)=GZM(J)
P1(IND1)=GZM(J)
WRITE(4,'(2X,4(E10.4,2X))') P(IND)
IND=IND+1
IND1=IND1+1
ENDIF
C PARAMETERS OF HYPERFINE COUPLING
READ(1,*)B9(J), AXM(J)
C CONVERTION FROM CM-1 TO S-1.RAD
AXM(J)=PI2*VL*AXM(J)
IF(B9(J).EQ.1)THEN
P(IND)=AXM(J)
P1(IND1)=AXM(J)
WRITE(4,'(2X,4(E10.4,2X))') P(IND)/VL/PI2
INDAPAR(J)=IND
IND=IND+1
IND1=IND1+1
ENDIF
READ(1,*)B10(J), AYM(J)
AYM(J)=PI2*VL*AYM(J)
IF(B10(J).EQ.1)THEN
P(IND)=AYM(J)
P1(IND1)=AYM(J)
WRITE(4,'(2X,4(E10.4,2X))') P(IND)/VL/PI2
INDAPER(J)=IND
IND=IND+1
IND1=IND1+1
ENDIF
READ(1,*)B21(J), AZM(J)
AZM(J)=PI2*VL*AZM(J)
IF(B21(J).EQ.1)THEN
P(IND)=AZM(J)
P1(IND1)=AZM(J)
WRITE(4,'(2X,4(E10.4,2X))') P(IND)/VL/PI2
INDAPER2(J)=IND
IND=IND+1
IND1=IND1+1
ENDIF
C PARAMETERS OF OUTER-SPHERE
READ(1,*)B13(J), DM(J)
DM(J)=DM(J)*1.E-8
IF(B13(J).EQ.1)THEN
P(IND)=DM(J)
P2(IND2)=DM(J)
WRITE(4,'(2X,4(E10.4,2X))') P(IND)/1.E-8
INDED(J)=IND
IND=IND+1
IND2=IND2+1
ENDIF
READ(1,*)B14(J), DDM(J)
IF(B14(J).EQ.1)THEN
P(IND)=DDM(J)
P2(IND2)=DDM(J)
WRITE(4,'(2X,4(E10.4,2X))') P(IND)
IND=IND+1
IND2=IND2+1
ENDIF
READ(1,*)B15(J), CONCM(J)
IF(B15(J).EQ.1)THEN
P(IND)=CONCM(J)
P2(IND2)=CONCM(J)
WRITE(4,'(2X,4(E10.4,2X))') P(IND)
IND=IND+1
IND2=IND2+1
ENDIF
C NUMBER OF DIFFERENT SITES
READ(1,*) ACQ
C PARAMETERS FOR DIFFERENT SITES
DO J=1,ACQ
READ(1,*)B4(J), (TAUMM(J,K),K=1,2)
IF(B4(J).GE.2)THEN
TM1=TAUMM(J,1)
TM2=TAUMM(J,2)
TM11(J)=TAUMM(J,1)
TM21(J)=TAUMM(J,2)
DO K=1,SET
TAUMM(J,K)=TM1*EXP(TM2/TEMP(K))
END DO
IF(B4(J).EQ.3) NPLUS2=1
ENDIF
IF(B4(J).EQ.1)THEN
P(IND)=TAUMM(J,1)
P2(IND2)=TAUMM(J,1)
WRITE(4,'(2X,4(E10.4,2X))') P(IND)
IND=IND+1
IND2=IND2+1
ENDIF
IF(B4(J).EQ.2)THEN
P(IND)=TM1
P2(IND2)=TM1
P(IND+1)=TM2
P2(IND2+1)=TM2
WRITE(4,'(2X,4(E10.4,2X))') TM1,TM2
IND=IND+2
IND2=IND2+2
ENDIF
READ(1,*)B5(J), AMOLFRAM(J)
AMOLFRAM(J)=AMOLFRAM(J)*1.E-3/111.
IF(B5(J).EQ.1)THEN
P(IND)=AMOLFRAM(J)
P2(IND2)=AMOLFRAM(J)
WRITE(4,'(2X,4(E10.4,2X))') P(IND)*111./1.E-3
INDAMOLFRA(J)=IND
IND=IND+1
IND2=IND2+1
ENDIF
READ(1,*)B6(J), RKM(J)
IF(B6(J).EQ.1)THEN
P(IND)=RKM(J)
P2(IND2)=RKM(J)
WRITE(4,'(2X,4(E10.4,2X))') P(IND)
IND=IND+1
IND2=IND2+1
ENDIF
READ(1,*)B16(J), ACONTM(J)
IF(B16(J).EQ.1)THEN
P(IND)=ACONTM(J)
P2(IND2)=ACONTM(J)
WRITE(4,'(2X,4(E10.4,2X))') P(IND)
IND=IND+1
IND2=IND2+1
ENDIF
READ(1,*)B7(J), THETAM(J)
IF(B7(J).EQ.1)THEN
P(IND)=THETAM(J)
P1(IND1)=THETAM(J)
WRITE(4,'(2X,4(E10.4,2X))') P(IND)
IND=IND+1
IND1=IND1+1
ENDIF
READ(1,*)B8(J), PHIM(J)
IF(B8(J).EQ.1)THEN
P(IND)=PHIM(J)
P1(IND1)=PHIM(J)
WRITE(4,'(2X,4(E10.4,2X))') P(IND)
IND=IND+1
IND1=IND1+1
ENDIF
C NUMBER OF FITTING PARAMETERS
NV=NV+B1(J)+B2(J)+B3(J)+B4(J)+B5(J)+B6(J)+B7(J)+B8(J)+B9(J)+
& B10(J)+B11(J)+B12(J)+B13(J)+B14(J)+B15(J)+B16(J)+B17(J)
& +B18(J)+B19(J)+B20(J)+B21(J)-NPLUS*3-NPLUS2*3
C NUMBER OF FITTING PARAMETERS IN EXTERNAL CICLE
NVEST=NVEST+B1(J)+B2(J)+B3(J)+B7(J)+B8(J)+B9(J)+
& B10(J)+B11(J)+B12(J)+B17(J)
& +B18(J)+B19(J)+B20(J)+B21(J)-NPLUS*3
NPLUS=0
NPLUS2=0
END DO
NVMEM=NV
DPARATOT=0.
EPARATOT=0.
APERTOT=0.
APERTOT2=0.
APARTOT=0.
ACONTOT=0.
ACONIND=0.
DO J=1,ACQ
DPARATOT=DPARAM(J)+DPARATOT
EPARATOT=EPARAM(J)+EPARATOT
APERTOT=AZM(J)+APERTOT
APERTOT2=AXM(J)+APERTOT2
APARTOT=AYM(J)+APARTOT
ACONTOT=ACONTM(J)+ACONTOT
END DO
C DEFINITION OF NMX: DIMENSION OF ENERGY MATRIX
IF(DPARATOT.EQ.0..AND.EPARATOT.EQ.0.
& AND.GX.EQ.GZ.AND.GX.EQ.GY.AND.SPIN.EQ.0.5)THEN
NMX=2.*(2*SI+1.)
ELSE
IF(APERTOT.EQ.0.AND.APARTOT.EQ.0.AND.APERTOT2.EQ.0.AND.
& GX.EQ.GZ.AND.GX.EQ.GY.AND.EPARATOT.EQ.0)THEN
SI=0.5
NMX=2*SPIN+1.
ELSE
NMX = (2*SI + 1)*(2*SPIN+1)
ENDIF
ENDIF
IF(IREL.NE.1) NMX = (2*SI + 1)*(2*SPIN+1)
IF(ACONTOT.NE.0) THEN
IF(APERTOT.NE.0.OR.APARTOT.NE.0.OR.APERTOT2.NE.0.OR.
& GX.NE.GZ.OR.GX.NE.GY.OR.EPARATOT.NE.0.OR.DPARATOT.NE.0)THEN
NMX = (2*SI + 1)*(2*SPIN+1)
ACONIND=1.
ENDIF
ENDIF
READ(1,*) (NPT(K),K=1, SET)
READ(1,*) FTOL
READ(1,*) ALFASTEP
C Z: FREQUENCIES, PP: RATE
NPTOT=0
DO K=1,SET
NPTOT=NPTOT+NPT(K)
END DO
DO 11
I=1,NPTOT
READ(1,*) Z(I),PP(I)
11 CONTINUE
CLOSE(1)
OO=10
IF (NPT(1).EQ.0)GOTO 250
C STARTING FITTING PROCEDURE
IF(NVEST.EQ.0)THEN
CALL FUNCZFS(P2,FUNC,NMX,NV)
NP=NV
N=NV
ITER=1000
CALL POWELLINT(P2,XI,N,NP,FTOL,ITER,FRET,NMX)
DO J=1,NP
WRITE(6,'(2X,4(E10.4,2X))') P2(J),2
END DO
ELSE
NP=NVEST
N=NVEST
ITER=1000
CALL POWELL(P1,XI,N,NP,FTOL,ITER,FRET,NMX)
ENDIF
WRITE(4,*) 'ERROR=', FRET/(NPTOT-NV)
DO J=1,ACQ
IF(B5(J).EQ.1)P2(INDAMOLFRA(J))=P2(INDAMOLFRA(J))*111/1.E-3
IF(B9(J).EQ.1)P1(INDAPAR(J))=P1(INDAPAR(J))/PI2/VL
IF(B10(J).EQ.1)P1(INDAPER(J))=P1(INDAPER(J))/PI2/VL
IF(B21(J).EQ.1)P1(INDAPER2(J))=P1(INDAPER2(J))/PI2/VL
IF(B11(J).EQ.1)P1(INDDPARA(J))=P1(INDDPARA(J))/PI2/VL
IF(B12(J).EQ.1)P1(INDEPARA(J))=P1(INDEPARA(J))/PI2/VL
IF(B13(J).EQ.1)P2(INDED(J))=P2(INDED(J))/1.E-8
IF(B19(J).EQ.1)P1(INDS4(J))=P1(INDS4(J))/PI2/VL
END DO
C WRITE RESULTS OF FITTING PROCEDURE
IF(NVEST.NE.0)THEN
JI=1
IF(B1(1).EQ.2)THEN
DO K=1,SET
WRITE(4,'(2X,4(E10.4,2X))')P1(JI)*EXP(P1(JI+1)/TEMP(K))
END DO
JI=JI+2
ENDIF
IF(B2(1).EQ.2)THEN
DO K=1,SET
WRITE(4,'(2X,4(E10.4,2X))')P1(JI)*EXP(P1(JI+1)/TEMP(K))
END DO
JI=JI+2
ENDIF
IF(B3(1).EQ.2)THEN
DO K=1,SET
WRITE(4,'(2X,4(E10.4,2X))')P1(JI)*EXP(P1(JI+1)/TEMP(K))
END DO
JI=JI+2
ENDIF
DO J=JI,NVEST
WRITE(4,'(2X,4(E10.4,2X))') P1(J)
END DO
ENDIF
DO JJ=1,ACQ
IF(NV-NVEST.NE.0)THEN
JI=1
IF(B4(JJ).EQ.2)THEN
DO K=1,SET
WRITE(4,'(2X,4(E10.4,2X))')P2(JI)*EXP(P2(JI+1)/TEMP(K))
END DO
JI=JI+2
ENDIF
DO J=JI,NV-NVEST
WRITE(4,'(2X,4(E10.4,2X))') P2(J)
END DO
ENDIF
END DO
WRITE(4,*) 'MAGN. FIELD, OBSED., CAL. '
DO 1 I=1,NPTOT
WRITE(4,'(2X,3(F8.3,2X))') Z(I),PP(I)/CONCM(1)*0.001,
& TPUNO(I)/CONCM(1)*0.001
1 CONTINUE
CLOSE(4)
DO J=1,ACQ
IF(B5(J).EQ.1)P2(INDAMOLFRA(J))=P2(INDAMOLFRA(J))/111*1.E-3
IF(B9(J).EQ.1)P1(INDAPAR(J))=P1(INDAPAR(J))*PI2*VL
IF(B10(J).EQ.1)P1(INDAPER(J))=P1(INDAPER(J))*PI2*VL
IF(B21(J).EQ.1)P1(INDAPER2(J))=P1(INDAPER2(J))*PI2*VL
IF(B11(J).EQ.1)P1(INDDPARA(J))=P1(INDDPARA(J))*PI2*VL
IF(B12(J).EQ.1)P1(INDEPARA(J))=P1(INDEPARA(J))*PI2*VL
IF(B13(J).EQ.1)P2(INDED(J))=P2(INDED(J))*1.E-8
IF(B19(J).EQ.1)P1(INDS4(J))=P1(INDS4(J))*PI2*VL
END DO
OO=0
250 CONTINUE
C CALCULATION OF THE CURVE
OPEN (44, FILE=FILENAME)
DO K=1,SET
NPT(K)=NUMPUN
END DO
INDEXSTAMPA=1
DO K=1,SET
DO I=1,NPT(K)
ZE = XMIN + (XMAX - XMIN)*FLOAT(I)/FLOAT(NPT(K))
ADD=0
DO IJK=1,K-1
ADD=ADD+NPT(IJK)
END DO
PP(I+1+ADD)=PP(I+ADD)
Z(I+ADD) = 10.**ZE/1000000.
END DO
END DO
NVEST=30+OO
CALL FUNCZFS(P1,FUNC,NMX,NV)
CLOSE(44)
STOP
END
SUBROUTINE FUNCZFS(P,FUNC,NMX,NV)
IMPLICIT REAL*8(A-H,O-Z)
DIMENSION P(NV)
DIMENSION XI2(21,21)
PARAMETER(PI2 = 6.2831853, VL = 2.9979D+10)
COMMON /SET/SET
COMMON /PPAR/ P2(21)
COMMON /RK10/ SPIN, SI
COMMON /WATER/ ACQ
COMMON /STAMPA/INDEXSTAMPA
COMMON /STEPGAMMA/ STEPGAMMA
COMMON /B4/ NVMEM,NPT(10),NPTOT
COMMON /B31/ TPUNO(500)
COMMON /B32/ PP(500),Z(500)
COMMON /C1M/DM(10),DDM(10),CONCM(10)
COMMON /IPERFM/AZM(10),AYM(10),AXM(10),THETAM(10),RKM(10),
& TAUCM(10),DPARAM(10),EPARAM(10),PHIM(10),S4M(10),GXM(10),
& GYM(10),GZM(10)
COMMON /TAU1M/ TAUS0M(10,10),TAURM(10,10),TAUVM(10,10),
& TAUMM(10,10)
COMMON /MOLFRAZM/ AMOLFRAM(10)
COMMON /CONTATM/ ACONTM(10)
COMMON /INDA/ INDDPARA(10),INDAPAR(10),INDAPER(10),INDEPARA(10)
& ,INDED(10),INDAMOLFRA(10),INDS4(10),INDAPER2(10)
COMMON /C1/D,DD,CONC
COMMON /IPERF/AZ,AY,AX,THETA,RK,TAUC,DPARA,EPARA,PHI,S4
COMMON /GTENSOR/ GX,GY,GZ
COMMON /TAU1/ TAUS0
COMMON /TAU/ TAUR,TAUV,TAUM
COMMON /MOLFRAZ/ AMOLFRA
COMMON /CONTAT/ ACONT
COMMON /GAMMAH/ GAMMAI
COMMON /TM/ TMUNO,TMUNOCONT,TMUNOCROSS
COMMON /CICLE/ NVEST
COMMON /TMAT/ TMAT(500,10),TMATCONT(500,10),TMATCROSS(500,10)
COMMON /BPARA/B1(10),B2(10),B3(10),B4(10),B5(10),B6(10),B7(10),
& B8(10),B9(10),B10(10),B11(10),B12(10),B13(10),B14(10),B15(10),
& B16(10),B17(10),B18(10),B19(10),B20(10),B21(10)
COMMON/TEMPERATURE/ TEMP(10)
DIMENSION TS1(10),TS2(10),TR1(10),TR2(10),TV1(10),TV2(10)
COMMON /TMSTART/ TM11(10),TM21(10)
C SET PARAMETERS
FB=0.
FBW=0.
IF(NVEST.LT.30)THEN
WRITE(6,*)'PARAMETERS: '
DO J=1,ACQ
IF(B5(J).EQ.1)P(INDAMOLFRA(J))=P(INDAMOLFRA(J))*111/1.E-3
IF(B9(J).EQ.1)P(INDAPAR(J))=P(INDAPAR(J))/PI2/VL
IF(B10(J).EQ.1)P(INDAPER(J))=P(INDAPER(J))/PI2/VL
IF(B21(J).EQ.1)P(INDAPER2(J))=P(INDAPER2(J))/PI2/VL
IF(B11(J).EQ.1)P(INDDPARA(J))=P(INDDPARA(J))/PI2/VL
IF(B12(J).EQ.1)P(INDEPARA(J))=P(INDEPARA(J))/PI2/VL
IF(B13(J).EQ.1)P(INDED(J))=P(INDED(J))/1.E-8
IF(B19(J).EQ.1)P(INDS4(J))=P(INDS4(J))/PI2/VL
END DO
DO I=1,NV
WRITE(6,'(2X,E10.4)')P(I)
END DO
DO J=1,ACQ
IF(B5(J).EQ.1)P(INDAMOLFRA(J))=P(INDAMOLFRA(J))/111*1.E-3
IF(B9(J).EQ.1)P(INDAPAR(J))=P(INDAPAR(J))*PI2*VL
IF(B10(J).EQ.1)P(INDAPER(J))=P(INDAPER(J))*PI2*VL
IF(B21(J).EQ.1)P(INDAPER2(J))=P(INDAPER2(J))*PI2*VL
IF(B11(J).EQ.1)P(INDDPARA(J))=P(INDDPARA(J))*PI2*VL
IF(B12(J).EQ.1)P(INDEPARA(J))=P(INDEPARA(J))*PI2*VL
IF(B13(J).EQ.1)P(INDED(J))=P(INDED(J))*1.E-8
IF(B19(J).EQ.1)P(INDS4(J))=P(INDS4(J))*PI2*VL
END DO
ENDIF
DO 223 K=1,SET
DO I=1, NPT(K)
IND=1
TPUNOTOT=0
J=1
TAUS0=TAUS0M(J,1)
TAUR=TAURM(J,1)
TAUV=TAUVM(J,1)
IF(TAUR.EQ.0.AND.TAUV.EQ.0.)TAUC=TAUCM(J)
AX=AXM(J)
AY=AYM(J)
AZ=AZM(J)
DPARA=DPARAM(J)
EPARA=EPARAM(J)
GX=GXM(J)
GY=GYM(J)
GZ=GZM(J)
S4=S4M(J)
IF(INDEXSTAMPA.EQ.0.OR.NVEST.EQ.40)THEN
D=DM(J)
DD=DDM(J)
CONC=CONCM(J)
ENDIF
C
PARAMETERS**************************************************************
IF(B1(J).EQ.1)THEN
TAUS0=P(IND)
IND=IND+1
ENDIF
IF(B1(J).EQ.2)THEN
TS1(K)=P(IND)
TS2(K)=P(IND+1)
IND=IND+2
ENDIF
IF(B2(J).EQ.1)THEN
TAUR=P(IND)
IND=IND+1
ENDIF
IF(B2(J).EQ.2)THEN
TR1(K)=P(IND)
TR2(K)=P(IND+1)
IND=IND+2
ENDIF
IF(B3(J).EQ.1)THEN
TAUV=P(IND)
IND=IND+1
ENDIF
IF(B3(J).EQ.2)THEN
TV1(K)=P(IND)
TV2(K)=P(IND+1)
IND=IND+2
ENDIF
IF(B1(J).EQ.1.AND.TAUR.EQ.0.AND.TAUV.EQ.0)THEN
TAUC=P(IND-1)
ENDIF
IF(B11(J).EQ.1)THEN
DPARA=P(IND)
IND=IND+1
ENDIF
IF(B12(J).EQ.1)THEN
EPARA=P(IND)
IND=IND+1
ENDIF
IF(B19(J).EQ.1)THEN
S4=P(IND)
IND=IND+1
ENDIF
IF(B17(J).EQ.1)THEN
GX=P(IND)
IND=IND+1
ENDIF
IF(B18(J).EQ.1)THEN
GY=P(IND)
IND=IND+1
ENDIF
IF(B20(J).EQ.1)THEN
GZ=P(IND)
IND=IND+1
ENDIF
IF(B9(J).EQ.1)THEN
AX=P(IND)
IND=IND+1
ENDIF
IF(B10(J).EQ.1)THEN
AY=P(IND)
IND=IND+1
ENDIF
IF(B21(J).EQ.1)THEN
AZ=P(IND)
IND=IND+1
ENDIF
IND2=1 !!!!!!!!!!!!!!!!!!!!!!!!!
IF(B13(J).EQ.1) IND2=IND2+1
IF(B14(J).EQ.1) IND2=IND2+1
IF(B15(J).EQ.1) IND2=IND2+1 !!!!!!!!!!!!!!!!!
DO J=1,ACQ
IF(INDEXSTAMPA.EQ.0.OR.NVEST.EQ.40)THEN
TAUM=TAUMM(J,1)
AMOLFRA=AMOLFRAM(J)
RK=RKM(J)
ACONT=ACONTM(J)
ENDIF
IF(INDEXSTAMPA.EQ.1)THEN !!!!!!!!!!!!!!!!!!!!!!!!!!!!!
IF(B4(J).EQ.0)TAUM=TAUMM(J,1)
IF(B4(J).EQ.1)THEN
TAUM=P2(IND2)
IND2=IND2+1
ENDIF
IF(B4(J).EQ.2)THEN
TM11(J)=P2(IND2)
TM21(J)=P2(IND2+1)
IND2=IND2+2
ENDIF
IF(B5(J).EQ.0)AMOLFRA=AMOLFRAM(J)
IF(B5(J).EQ.1)THEN
AMOLFRA=P2(IND2)
IND2=IND2+1
ENDIF
IF(B6(J).EQ.0)RK=RKM(J)
IF(B6(J).EQ.1)THEN
RK=P2(IND2)
IND2=IND2+1
ENDIF
IF(B16(J).EQ.0)ACONT=ACONTM(J)
IF(B16(J).EQ.1)THEN
ACONT=P2(IND2)
IND2=IND2+1
ENDIF
ENDIF !!!!!!!!!!!!!!!!!!!!!!!!!!!!
THETA=THETAM(J)
PHI=PHIM(J)
IF(B7(J).EQ.1)THEN
THETA=P(IND)
IND=IND+1
ENDIF
IF(B8(J).EQ.1)THEN
PHI=P(IND)
IND=IND+1
ENDIF
IF(B1(J).EQ.2)TAUS0=TS1(1)*EXP(TS2(1)/TEMP(K))
IF(B2(J).EQ.2)TAUR=TR1(1)*EXP(TR2(1)/TEMP(K)) !stokes
IF(B3(J).EQ.2)TAUV=TV1(1)*EXP(TV2(1)/TEMP(K))
IF(B4(J).EQ.2)TAUM=TM11(J)*EXP(TM21(J)/TEMP(K))
IF(B1(J).EQ.3)TAUS0=TAUS0M(J,K)
IF(B2(J).EQ.3)TAUR=TAURM(J,K)
IF(B3(J).EQ.3)TAUV=TAUVM(J,K)
IF(B4(J).EQ.3)TAUM=TAUMM(J,K)
C*******************************************************************************
ADD=0
DO IJK=1,K-1
ADD=ADD+NPT(IJK)
END DO
IPLUS=I+ADD
C PROTON LARMOR FREQUENCY
BZ=Z(IPLUS)*1000000
C CONSTANTS IN DIPOLAR RELAXATION
RK1=10.*2.46502D-52/(2.6752D8)**2*GAMMAI**2*(RK*1.E-10)**(-6)
IF (AMOLFRA.EQ.0.)RK1=10./(SPIN*(SPIN+1.)*2./15.)/TAUS0*1.E-9
IF(RK.EQ.0) RK1=0.
C CONSTANT OF CONTACT RELAXATION
CONA=(ACONT*6.28*1.E6*1.0546E-34)
IF(TAUC.EQ.0.AND.TAUS0.EQ.0) GOTO 56
CALL GAUINT (BZ,TAUM,NMX)
C STORE CONTRIBUTIONS FOR TUNO
TMAT(IPLUS,J)=TMUNO
TMATCONT(IPLUS,J)=TMUNOCONT
TMATCROSS(IPLUS,J)=TMUNOCROSS
C CALCULATION OF TPUNO
TMUNO=TMUNO*RK1+CONA*CONA*TMUNOCONT+SQRT(RK1)*CONA*TMUNOCROSS
TPUNO1=1./(1./TMUNO+TAUM)*AMOLFRA
IF(AMOLFRA.EQ.0.)TPUNO1=1./(1./TMUNO+TAUM)
TPUNO(IPLUS)=TPUNO1/CONC*0.001
TPUNOTOT=TPUNOTOT+TPUNO(IPLUS)
END DO
56 CONTINUE
C CALCULATION OF OUTER-SPHERE CONTRIBUTION
TERM=0
IF(DD.NE.0)THEN
V=BZ*2.*3.14
TAUD=D**2/DD
CZ=SQRT(2*V*TAUD)
ZZ=SQRT(2*V*TAUD*658.)
GEI=(1.+5.*CZ/8.+CZ**2/8.)/(1.+CZ+CZ**2/2.+CZ**3/6.+4.*
& CZ**4/81.+CZ**5/81.+CZ**6/648.)
GES=(1.+5.*ZZ/8.+ZZ**2/8.)/(1.+ZZ+ZZ**2/2.+ZZ**3/6.+4.*
& ZZ**4/81.+ZZ**5/81.+ZZ**6/648.)
PRIMO=32.*3.14/405.*(2.6752E4*1.7608E7*1.0546E-27)**2*6.022E20
SECONDO=SPIN*(SPIN+1)*CONC/(D*DD)
TERZO=(3.*GEI+7.*GES)
TERM=(PRIMO*SECONDO*TERZO)
ENDIF
TPUNOTOT=TPUNOTOT+TERM
TPUNO(IPLUS)=TPUNOTOT
C DIFFERENCE BETWEEN EXPERIMENTAL AND FITTING VALUES
FB=((PP(IPLUS)-TPUNO(IPLUS))**2)/PP(IPLUS)+FB
FBW=SQRT((PP(IPLUS)-TPUNO(IPLUS))**2)/PP(IPLUS)+FBW
IF (STEPGAMMA.NE.1) WRITE(6,'(2X,2(E10.4))')Z(IPLUS),TPUNO
(IPLUS)
IF(INDEXSTAMPA.EQ.1)WRITE(44,'(2X,2(E10.4))')6.283d+6*Z(IPLUS)/
& (GAMMAI) ,TPUNO(IPLUS)
END DO
223 CONTINUE
FUNC=FBW
IF(INDEXSTAMPA.EQ.0)WRITE(6,*)' ERROR',
& '=',FBW/(NPTOT-NVMEM),'**********'
IF(NV.NE.NVMEM)THEN
C INTERNAL FITTING PROCEDURE
NP2=NVMEM-NV
N2=NVMEM-NV
ITER2=1000
FRET2=0.
CALL POWELLINT(P2,XI2,N2,NP2,FTOL,ITER2,FRET2,NMX)
WRITE(6,*)' ERROR',
& '=', FRET2/(NPTOT-NVMEM)
DO J=1,NP2
WRITE(6,'(2X,E10.4)')P2(J)
END DO
FUNC=FRET2
ENDIF
RETURN
END
SUBROUTINE FUNCINT(P,FUNC,NMX,NV)
IMPLICIT REAL*8(A-H,O-Z)
DIMENSION P(NV)
PARAMETER(PI2 = 6.2831853, VL = 2.9979D+10)
COMMON /SET/SET
COMMON /RK10/ SPIN, SI
COMMON /WATER/ ACQ
COMMON /STAMPA/INDEXSTAMPA
COMMON /STEPGAMMA/ STEPGAMMA
COMMON /B4/ NVMEM,NPT(10),NPTOT
COMMON /B31/ TPUNO(500)
COMMON /B32/ PP(500),Z(500)
COMMON /C1M/DM(10),DDM(10),CONCM(10)
COMMON /IPERFM/AZM(10),AYM(10),AXM(10),THETAM(10),RKM(10),
& TAUCM(10),DPARAM(10),EPARAM(10),PHIM(10),S4M(10),GXM(10),
& GYM(10),GZM(10)
COMMON /TAU1M/ TAUS0M(10,10),TAURM(10,10),TAUVM(10,10),
& TAUMM(10,10)
COMMON /MOLFRAZM/ AMOLFRAM(10)
COMMON /CONTATM/ ACONTM(10)
COMMON/INDA/ INDDPARA(10),INDAPAR(10),INDAPER(10),INDEPARA(10)
& ,INDED(10),INDAMOLFRA(10),INDS4(10),INDAPER2(10)
COMMON /C1/D,DD,CONC
COMMON /IPERF/AZ,AY,AX,THETA,RK,TAUC,DPARA,EPARA,PHI,S4
COMMON /GTENSOR/ GX,GY,GZ
COMMON /TAU1/ TAUS0
COMMON /TAU/ TAUR,TAUV,TAUM
COMMON /MOLFRAZ/ AMOLFRA
COMMON /CONTAT/ ACONT
COMMON /GAMMAH/ GAMMAI
COMMON /TM/ TMUNO,TMUNOCONT,TMUNOCROSS
COMMON /CICLE/ NVEST
COMMON /TMAT/ TMAT(500,10),TMATCONT(500,10),TMATCROSS(500,10)
COMMON /BPARA/B1(10),B2(10),B3(10),B4(10),B5(10),B6(10),B7(10),
& B8(10),B9(10),B10(10),B11(10),B12(10),B13(10),B14(10),B15(10),
& B16(10),B17(10),B18(10),B19(10),B20(10),B21(10)
COMMON/TEMPERATURE/ TEMP(10)
DIMENSION TM1(10),TM2(10)
C SET PARAMETERS
FB=0.
FBW=0.
DO 223 K=1,SET
DO I=1, NPT(K)
IND=1
TPUNOTOT=0
J=1
C PARAMETERS OF INTERNAL
FITTING*********************************************
IF(B13(J).EQ.1)THEN
D=P(IND)
IND=IND+1
ELSE
D=DM(J)
ENDIF
IF(B14(J).EQ.1)THEN
DD=P(IND)
IND=IND+1
ELSE
DD=DDM(J)
ENDIF
IF(B15(J).EQ.1)THEN
CONC=P(IND)
IND=IND+1
ELSE
CONC=CONCM(J)
ENDIF
DO J=1,ACQ
IF(B4(J).EQ.1)THEN
TAUM=P(IND)
IND=IND+1
ENDIF
IF(B4(J).EQ.2)THEN
TM1(1)=P(IND)
TM2(1)=P(IND+1)
IND=IND+2
ENDIF
IF(B4(J).EQ.0) TAUM=TAUMM(J,1)
IF(B5(J).EQ.1)THEN
AMOLFRA=P(IND)
IND=IND+1
ELSE
AMOLFRA=AMOLFRAM(J)
ENDIF
IF(B6(J).EQ.1)THEN
RK=P(IND)
IND=IND+1
ELSE
RK=RKM(J)
ENDIF
IF(B16(J).EQ.1)THEN
ACONT=P(IND)
IND=IND+1
ELSE
ACONT=ACONTM(J)
ENDIF
IF(B4(J).EQ.2)TAUM=TM1(1)*EXP(TM2(1)/TEMP(K))
IF(B4(J).EQ.3)TAUM=TAUMM(J,K)
C*******************************************************************************
ADD=0
DO IJK=1,K-1
ADD=ADD+NPT(IJK)
END DO
IPLUS=I+ADD
BZ=Z(IPLUS)*1000000.
RK1=10.*2.46502D-52/(2.6752D8)**2*GAMMAI**2*(RK*1.E-10)**(-6)
IF (AMOLFRA.EQ.0.)RK1=10./(SPIN*(SPIN+1.)*2./15.)/TAUS0*1.E-9
IF(RK.EQ.0) RK1=0.
CONA=(ACONT*6.28*1.E6*1.0546E-34)
C READ CALCULATED CONTRIBUTIONS TO TMUNO
TMUNO=TMAT(IPLUS,J)
TMUNOCONT=TMATCONT(IPLUS,J)
TMUNOCROSS=TMATCROSS(IPLUS,J)
C CALCULATION OF TMUNO
TMUNO=TMUNO*RK1+CONA*CONA*TMUNOCONT+SQRT(RK1)*CONA*TMUNOCROSS
TPUNO1=1./(1./TMUNO+TAUM)*AMOLFRA
IF(AMOLFRA.EQ.0.)TPUNO1=1./(1./TMUNO+TAUM)
TPUNO(IPLUS)=TPUNO1/CONC*0.001
TPUNOTOT=TPUNOTOT+TPUNO(IPLUS)
END DO
56 CONTINUE
C CALCULATION OF OUTER-SPHERE CONTRIBUTION
TERM=0
IF(DD.NE.0)THEN
V=BZ*2.*3.14
TAUD=D**2/DD
CZ=SQRT(2*V*TAUD)
ZZ=SQRT(2*V*TAUD*658.)
GEI=(1.+5.*CZ/8.+CZ**2/8.)/(1.+CZ+CZ**2/2.+CZ**3/6.+4.*
& CZ**4/81.+CZ**5/81.+CZ**6/648.)
GES=(1.+5.*ZZ/8.+ZZ**2/8.)/(1.+ZZ+ZZ**2/2.+ZZ**3/6.+4.*
& ZZ**4/81.+ZZ**5/81.+ZZ**6/648.)
PRIMO=32.*3.14/405.*(2.6752E4*1.7608E7*1.0546E-27)**2*6.022E20
SECONDO=SPIN*(SPIN+1)*CONC/(D*DD)
TERZO=(3.*GEI+7.*GES)
TERM=(PRIMO*SECONDO*TERZO)
ENDIF
TPUNOTOT=TPUNOTOT+TERM
C DIFFERENCES BETWEEN EXPERIMENTAL AND FITTED VALUES
TPUNO(IPLUS)=TPUNOTOT
FB=((PP(IPLUS)-TPUNO(IPLUS))**2)/PP(IPLUS)+FB
FBW=SQRT((PP(IPLUS)-TPUNO(IPLUS))**2)/PP(IPLUS)+FBW
!IF (STEPGAMMA.NE.1) WRITE(6,'(2X,2(E10.4))')Z(IPLUS),TPUNO
(IPLUS)
IF(INDEXSTAMPA.EQ.1) WRITE(44,'(2X,2(E10.4))') 6.283d+6*Z
(IPLUS)/
& (GAMMAI), TPUNO(IPLUS)
END DO
223 CONTINUE
FUNC=FBW
RETURN
END
FUNCTION E(BZ,BETA,THETA,TAUC,NMX,PHI,GAMMA)
IMPLICIT REAL*8(A-H,O-Z)
COMMON /A3/ T11,T12,T13
COMMON /T1T2/ IREL
COMMON /ECOM/ ECONT,ECROSS
COMMON /STEPGAMMA/ STEPGAMMA
COMMON /TOT/ DPARATOT,EPARATOT
OMI=BZ*6.2831853
CALL DIAG(BETA,GAMMA,BZ,NMX)
CCCCCCCCCCCCCCCCCCCCCCCC modificare CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
IF(IREL.EQ.1.AND.STEPGAMMA.GT.1)CALL TUNO(BETA,OMI,THETA,
& TAUC,NMX,PHI,GAMMA)
IF(IREL.EQ.1.AND.STEPGAMMA.EQ.1)CALL TUNOISO(BETA,OMI,THETA,
& TAUC,NMX)
! IF(IREL.EQ.1)CALL TUNO(BETA,OMI,THETA,TAUC,NMX,PHI,GAMMA)
IF(IREL.EQ.2)CALL TDUE(BETA,OMI,THETA,TAUC,NMX,PHI,GAMMA)
E=T11*SIN(BETA)
ECONT=T12*SIN(BETA)
ECROSS=T13*SIN(BETA)
RETURN
END
SUBROUTINE DIAG(BETA,GAMMA,BZ,NMX)
IMPLICIT REAL*8(A-H,O-Z)
PARAMETER(PI2 = 6.2831853, VL = 2.9979D+10)
common /butto/ taurb,tausb
COMMON /RK10/ SPIN, SI
COMMON /T1T2/ IREL
COMMON /INDEX/ INDEX
COMMON /STEPGAMMA/ STEPGAMMA
COMMON /TAUDELTA/ TAUDELTA
COMMON /TAUE/ TAUE
COMMON /CONTAT/ ACONT
COMMON /TOT/ DPARATOT,EPARATOT,APARTOT,APERTOT,APERTOT2,ACONIND
COMMON /IPERF/AZ,AY,AX,THETA,RK,TAUC,DPARA,EPARA,PHI,S4
COMMON /GTENSOR/ GX,GY,GZ
COMMON /TAU/ TAUR,TAUV,TAUM
COMMON /TAU1/ TAUS0
COMMON /AOLD/ OMOLD(10000),COLD(10000,4)
dimension cold1(10000,4),cold2(10000,4)
COMPLEX*16 C(100,100,19)
COMPLEX*16 Cn1(100,100,19)
COMPLEX*16 Cn2(100,100,19)
COMMON /A1/ OM(1000,1000),C
PARAMETER(MBRANC=90)
COMPLEX*16 SZ(MBRANC,MBRANC),SP(MBRANC,MBRANC),SM(MBRANC,MBRANC)
DIMENSION CO1(MBRANC),CO2(MBRANC), CO3(MBRANC)
COMPLEX*16 SZT(MBRANC,MBRANC),SPT(MBRANC,MBRANC)
COMPLEX*16 SMT(MBRANC,MBRANC)
COMPLEX*16 TZ,TP,TM
COMPLEX*16 cpp,cpm,cpz,cmp,cmm,cmz,czz,czp,czm
COMPLEX*16 ammcp1,ammcp2,ammcp3,ammcp4,ammcp5,ammcp6,ammcp7,
& ammcp8,ammcp9
COMPLEX*16 ammcm1,ammcm2,ammcm3,ammcm4,ammcm5,ammcm6,ammcm7,
& ammcm8,ammcm9
COMPLEX*16 ammcz1,ammcz2,ammcz3,ammcz4,ammcz5,ammcz6,ammcz7,
& ammcz8,ammcz9
COMPLEX*16 primo,secondo,terzo,quarto,quinto,sesto,settimo,
& ottavo,onono
COMPLEX*16 a1,b1,c1,d1,e1,f1,g1,h1,ai1,al1,trhoa,trhob,trhoc
COMPLEX*16 a,b,cp,a11,a12,a13,a21,a22,a23,a31,a32,a33
COMPLEX*16 a44,a55,a66,a77,a88,a99,a45,a54,a78,a87
COMPLEX*16 t10,t21,t3,t4,t5,t6,t78,t8,t11,t12
COMPLEX*16 rpzpz,rzmzm,r1223,rpmpm
complex*16 ak11,ak12,ak13,ak21,ak22,ak23,ak31,ak32,ak33
complex*16 cz1,cz2,cz3,cz4,cz5,cz6,cz7,cz8,cz9
complex*16 cp1,cp2,cp3,cp4,cp5,cp6,cp7,cp8,cp9
complex*16 cm1,cm2,cm3,cm4,cm5,cm6,cm7,cm8,cm9
COMPLEX*16 rhoa,rhob,rhoc,rt1
common /stoccolma/ disp
DIMENSION WR( MBRANC )
COMPLEX*16 AR(MBRANC,MBRANC),ZR(MBRANC,MBRANC )
COMPLEX*16 aaa(3,3)
DIMENSION ARR(MBRANC,MBRANC),ARI(MBRANC,MBRANC)
DIMENSION ZRR(MBRANC,MBRANC),ZRI(MBRANC,MBRANC)
DIMENSION WK1(MBRANC),WK2(MBRANC),WK3(MBRANC)
integer lda,num,ipvt(3),info,job
complex*16 det(2),work(3),zrn(3,3)
COMMON /GAMMAH/GAMMAI
complex a44zfs,a55zfs,a44zee,a55zee
CT=COS(BETA)
ST=SIN(BETA)
C CALCULATION OF CORRELATION TIME
WI=2*3.1416*BZ
WS=658.2*WI
IF(TAUV.EQ.0.AND.TAUR.EQ.0)THEN
TAUC=TAUS0
IF(TAUM.EQ.0)THEN
TAUE=TAUS0
ELSE
TAUE=1./(1./TAUS0+1./TAUM)
ENDIF
ELSE
STI=WS**2*TAUV**2
IF (TAUDELTA.EQ.2)THEN
delta=taus0*VL*PI2
RTAUS=2.*(TAUS0*VL*PI2)**2*(4.*SPIN*(SPIN+1)-3)/50.*
& (TAUV/(1+STI)+4*TAUV/(1+4*STI))
! write(6,*)rtaus
! stop !cance
ELSE
delta=0
RTAUS=(0.2/TAUS0)*(1./(1.+STI)+4./(1.+4.*STI))
ENDIF
IF(TAUR.EQ.0)THEN
RTAUC=RTAUS
ELSE
RTAUC=RTAUS+1./TAUR
ENDIF
IF(TAUM.NE.0)THEN
RTAUC=RTAUC+1./TAUM
ENDIF
TAUC=1./RTAUC
IF (ACONT.NE.0)THEN
IF (TAUM.EQ.0)THEN
RTAUE=RTAUS
ELSE
RTAUE=RTAUS+1./TAUM
ENDIF
TAUE=1./RTAUE
ENDIF
ENDIF
IF (STEPGAMMA.GT.1)THEN
COEFFH=-1.
ELSE
COEFFH=1.
ENDIF
IF(ACONIND.EQ.1.)GO TO 456
IF (DPARATOT.EQ.0..AND.EPARATOT.EQ.0..AND.SPIN.EQ.0.5.AND.
& GX.EQ.GY.AND.GX.EQ.GZ.AND.IREL.EQ.1)THEN
C MATRIX OF ENERGY FOR HYPERFINE COUPLING
X=BZ*3.1415927*658.2
ZC=X*CT
ZS=X*ST
DO 200 I=1,(2.*SI+1.)*2.
DO 200 J=1,(2.*SI+1.)
*2.
200 AR(I,J)
=0.
SSI = SI*(SI + 1.)
DO I = 1, (2.*SI+1.)*2., 2
COEF = SI - (I - 1)/2
AR(I,I) = ZC*GZ + (SI-I/2)*AZ/2.
AR(I+1,I+1) = -(ZC*GZ + (SI-I/2)*AZ/2.)
AR(I,I+1) = COEFFH*ZS*GY
AR(I+1,I) = COEFFH*ZS*GY
AR(I+1,I+2) = 0.5*(AX+AY)/2.*SQRT(SSI-(COEF-1.)*COEF)
AR(I+2,I+1) = 0.5*(AX+AY)/2.*SQRT(SSI-(COEF-1.)*COEF)
END DO
IF (INDEX.EQ.1)THEN
WRITE(6,*)'DIM. MATRIX', NMX
OPEN(UNIT=17,FILE='MAT')
DO I=1,(2.*SI+1)*(2.*SPIN+1)
DO J=1,(2.*SI+1)*(2.*SPIN+1)
WRITE(17,*)AR(I,J)
END DO
WRITE(17,*)' '
END DO
CLOSE(17)
ENDIF
INDEX=INDEX+1
C DIAGONALISATION OF THE MATRIX OF ENERGY
DO 45 I=1,NMX
DO 45 J=1,NMX
ARR(I,J)=REAL(AR(I,J))
ARI(I,J)=IMAG(AR(I,J))
45 CONTINUE
CALL F02AXF(ARR,MBRANC,ARI,MBRANC,NMX,WR,ZRR,MBRANC,ZRI,MBRANC
$ ,WK1,WK2,WK3,0)
DO 46 I=1,NMX
DO 46 J=1,NMX
ZR(I,J)=CMPLX(ZRR(I,J),ZRI(I,J))
46 CONTINUE
I=1
OM(1,1)=0.
OMOLD(1)=0.
DO 700 K=1,NMX
DO 700 L=1,NMX
IF (K.EQ.L)GO TO 700
I=I+1
OM(K,L)=WR(K)-WR(L)
C DIFFERENCES IN ENERGY LEVELS
OMOLD(I)=WR(K)-WR(L)
700 CONTINUE
C CALCULATION OF CORRELATION FUNCTIONS
DO 400
K=1,NMX
DO 400
L=1,NMX
TZ=0
DO 1500
J=1,NMX
TZ=-((-1.)**J)*ZR(J,K)*CONJG(ZR(J,L))
+TZ
1500
CONTINUE
SZ(K,L)=TZ/
2.
SZT(K,L)=SZ
(K,L)
TP=0
DO J=1,NMX,2
TP=ZR(J,K)*CONJG(ZR(J+1,L))+TP
END DO
SP(K,L)
=TP
SPT(K,L)=SP
(K,L)
TM=0
DO J=2,NMX,2
TM=ZR(J,K)*CONJG(ZR(J-1,L))
+TM
END DO
SM(K,L)
=TM
SMT(K,L)=SM
(K,L)
400
CONTINUE
GO TO 567
ENDIF
IF (APARTOT.EQ.0.AND.APERTOT.EQ.0.AND.APERTOT2.EQ.0.
& AND.EPARATOT.EQ.0.AND.GX.EQ.GY.AND.GX.EQ.GZ.AND.IREL.EQ.1)THEN
IF (INDEX.EQ.1)THEN
WRITE(6,*)'DIM. MATRIX', NMX
ENDIF
INDEX=INDEX+1
C MATRIX OF ENERGY IN ZERO FIELD SPLITTING
X=BZ*2*3.1415927*658.2
ZC=X*CT
ZS=X*ST
DO 5200 I=1,NMX
DO 5200 J=1,NMX
5200 AR(I,J)=0.
S = FLOAT(NMX - 1)/2.
SS = S*(S + 1.)
DO I = 1, NMX
COEF = S - DFLOAT(I - 1)
AR(I,I) = COEF*ZC*GZ + DPARA*(COEF**2 - SS/3.)
IF(I.LT.NMX) THEN
AR(I,I+1) = COEFFH*0.5*ZS*GY*SQRT(SS-(COEF-1.)*COEF)
AR(I+1,I) = AR(I,I+1)
END IF
END DO
IF (INDEX.EQ.2)THEN
OPEN(UNIT=17,FILE='MAT')
DO I=1,(2.*SI+1)*(2.*SPIN+1)
DO J=1,(2.*SI+1)*(2.*SPIN+1)
WRITE(17,*)AR(I,J)
END DO
WRITE(17,*)' '
END DO
CLOSE(17)
ENDIF
INDEX=INDEX+1
DO 145 I=1,NMX
DO 145 J=1,NMX
ARR(I,J)=REAL(AR(I,J))
ARI(I,J)=IMAG(AR(I,J))
145 CONTINUE
CALL F02AXF(ARR,MBRANC,ARI,MBRANC,NMX,WR,ZRR,MBRANC,ZRI,MBRANC
$ ,WK1,WK2,WK3,0)
DO 146 I=1,NMX
DO 146 J=1,NMX
ZR(I,J)=CMPLX(ZRR(I,J),ZRI(I,J))
146 CONTINUE
I=1
OM(1,1)=0.
OMOLD(1)=0.
DO 570 K=1,NMX
DO 570 L=1,NMX
IF (K.EQ.L)GO TO 570
I=I+1
OM(K,L)=WR(K)-WR(L)
OMOLD(I)=WR(K)-WR(L)
570 CONTINUE
C PER SPIN (SZ) DIVERSI DA 1/2
DO J=1,NMX
CO1(J)=(NMX-(2*J-1))/2.
END DO
DO J=1,NMX-1
CO2(J)=SQRT(SS-CO1(J+1)*(CO1(J+1)+1.))
END DO
DO J=2,NMX
CO3(J)=SQRT(SS-CO1(J-1)*(CO1(J-1)-1.))
END DO
omega1=abs(wr(2)-wr(3)) !omold(5)) !ZFS
omega2=abs(wr(1)-wr(2)) !omold(3)) !ZFS
omega3=abs(wr(3)-wr(1)) !omold(2)) !ZFS
cmp=zrr(3,3)
czp=zrr(2,3)
cpp=zrr(1,3)
cmm=zrr(3,1) !solo Zeeman
czm=zrr(2,1)
cpm=zrr(1,1)
cmz=zrr(3,2)
czz=zrr(2,2)
cpz=zrr(1,2)
bfield=bz*6.283/GAMMAI
C*****************************************************************************
C*********** CACLULATION OF ELECTRONIC R2
************************************
C*********** IN AXIAL ROUTINE
************************************
C*****************************************************************************
!nuovo calcolo R2
!++00
a1=cpp*conjg(cpp)-2*czp*conjg(czp)+cmp*conjg(cmp)
b1=cpp*conjg(czp)-czp*conjg(cmp)
c1=conjg(czp)*cmp-czp*conjg(cpp)
d1=cpp*conjg(cmp)
e1=cmp*conjg(cpp)
f1=cpz*conjg(cpz)-2*czz*conjg(czz)+cmz*conjg(cmz)
g1=cpz*conjg(czz)-czz*conjg(cmz)
h1=cmz*conjg(czz)-czz*conjg(cpz)
ai1=cmz*conjg(cpz)
al1=cpz*conjg(cmz)
t10=a1*f1/6.-(c1*g1+b1*h1)/2.+e1*ai1+d1*al1
!0++0
a1=cpz*conjg(cpp)-2*czz*conjg(czp)+cmz*conjg(cmp)
b1=cpz*conjg(czp)-czz*conjg(cmp)
c1=cmz*conjg(czp)-czz*conjg(cpp)
d1=cpz*conjg(cmp)
e1=cmz*conjg(cpp)
f1=cpp*conjg(cpz)-2*czp*conjg(czz)+cmp*conjg(cmz)
g1=cpp*conjg(czz)-czp*conjg(cmz)
h1=cmp*conjg(czz)-czp*conjg(cpz)
ai1=cpp*conjg(cmz)
al1=cmp*conjg(cpz)
t3=a1*f1/6.-(c1*g1+b1*h1)/2.+e1*ai1+d1*al1
!0000
a1=cpz*conjg(cpz)-2*czz*conjg(czz)+cmz*conjg(cmz)
b1=cpz*conjg(czz)-czz*conjg(cmz)
c1=cmz*conjg(czz)-czz*conjg(cpz)
d1=cmz*conjg(cpz)
e1=cpz*conjg(cmz)
f1=cpz*conjg(cpz)-2*czz*conjg(czz)+cmz*conjg(cmz)
g1=cpz*conjg(czz)-czz*conjg(cmz)
h1=cmz*conjg(czz)-czz*conjg(cpz)
ai1=cmz*conjg(cpz)
al1=cpz*conjg(cmz)
t4=a1*f1/6.-(c1*g1+b1*h1)/2.+e1*ai1+d1*al1
!0--0
a1=cpz*conjg(cpm)-2*czz*conjg(czm)+cmz*conjg(cmm)
b1=cpz*conjg(czm)-czz*conjg(cmm)
c1=cmz*conjg(czm)-czz*conjg(cpm)
d1=cpz*conjg(cmm)
e1=cmz*conjg(cpm)
f1=cpm*conjg(cpz)-2*czm*conjg(czz)+cmm*conjg(cmz)
g1=cpm*conjg(czz)-czm*conjg(cmz)
h1=cmm*conjg(czz)-czm*conjg(cpz)
ai1=cpm*conjg(cmz)
al1=cmm*conjg(cpz)
t5=a1*f1/6.-(c1*g1+b1*h1)/2.+e1*ai1+d1*al1
!++++
a1=cpp*conjg(cpp)-2*czp*conjg(czp)+cmp*conjg(cmp)
b1=cpp*conjg(czp)-czp*conjg(cmp)
c1=conjg(czp)*cmp-czp*conjg(cpp)
d1=cpp*conjg(cmp)
e1=cmp*conjg(cpp)
f1=cpp*conjg(cpp)-2*czp*conjg(czp)+cmp*conjg(cmp)
g1=cpp*conjg(czp)-czp*conjg(cmp)
h1=conjg(czp)*cmp-czp*conjg(cpp)
ai1=cpp*conjg(cmp)
al1=cmp*conjg(cpp)
t6=a1*f1/6.-(c1*g1+b1*h1)/2.+e1*ai1+d1*al1
!-++-
a1=cpp*conjg(cpm)-2*czp*conjg(czm)+cmp*conjg(cmm)
b1=cpp*conjg(czm)-czp*conjg(cmm)
c1=cmp*conjg(czm)-czp*conjg(cpm)
d1=cpp*conjg(cmm)
e1=cmp*conjg(cpm)
f1=cpm*conjg(cpp)-2*czm*conjg(czp)+cmm*conjg(cmp)
g1=cpm*conjg(czp)-czm*conjg(cmp)
h1=cmm*conjg(czp)-czm*conjg(cpp)
ai1=cpm*conjg(cmp)
al1=cmm*conjg(cpp)
t78=a1*f1/6.-(c1*g1+b1*h1)/2.+e1*ai1+d1*al1
rpzpz=-(2*t10*tauv
& -t6*tauv
& -2*t3*
& tauv/(1+omega1**2*tauv**2)
& -t78*
& tauv/(1+omega3**2*tauv**2)
& -t4*tauv
& -t5*
& tauv/(1+omega2**2*tauv**2))
!----
a1=cpm*conjg(cpm)-2*czm*conjg(czm)+cmm*conjg(cmm)
b1=cpm*conjg(czm)-czm*conjg(cmm)
c1=cmm*conjg(czm)-czm*conjg(cpm)
d1=cpm*conjg(cmm)
e1=cmm*conjg(cpm)
f1=cpm*conjg(cpm)-2*czm*conjg(czm)+cmm*conjg(cmm)
g1=cpm*conjg(czm)-czm*conjg(cmm)
h1=cmm*conjg(czm)-czm*conjg(cpm)
ai1=cpm*conjg(cmm)
al1=cmm*conjg(cpm)
t12=a1*f1/6.-(c1*g1+b1*h1)/2.+e1*ai1+d1*al1
!00--
a1=cpz*conjg(cpz)-2*czz*conjg(czz)+cmz*conjg(cmz)
b1=cpz*conjg(czz)-czz*conjg(cmz)
c1=cmz*conjg(czz)-czz*conjg(cpz)
d1=cmz*conjg(cpz)
e1=cpz*conjg(cmz)
f1=cpm*conjg(cpm)-2*czm*conjg(czm)+cmm*conjg(cmm)
g1=cpm*conjg(czm)-czm*conjg(cmm)
h1=cmm*conjg(czm)-czm*conjg(cpm)
ai1=cpm*conjg(cmm)
al1=cmm*conjg(cpm)
t11=a1*f1/6.-(c1*g1+b1*h1)/2.+e1*ai1+d1*al1
rzmzm=-(2*t11*tauv
& -t3*
& tauv/(1+omega1**2*tauv**2)
& -t4*tauv
& -2*t5*
& tauv/(1+omega2**2*tauv**2)
& -t78*
& tauv/(1+omega3**2*tauv**2)
& -t12*tauv)
a1=cpp*conjg(cpz)-2*czp*conjg(czz)+cmp*conjg(cmz)
b1=cpp*conjg(czz)-czp*conjg(cmz)
c1=cmp*conjg(czz)-czp*conjg(cpz)
d1=cpp*conjg(cmz)
e1=cmp*conjg(cpz)
f1=cpm*conjg(cpz)-2*czm*conjg(czz)+cmm*conjg(cmz)
g1=cpm*conjg(czz)-czm*conjg(cmz)
h1=cmm*conjg(czz)-czm*conjg(cpz)
ai1=cpm*conjg(cmz)
al1=cmm*conjg(cpz)
t8=a1*f1/6.-(c1*g1+b1*h1)/2.+e1*ai1+d1*al1
r1223zee=-t8*
& (tauv/(1+omega1**2*tauv**2)+
& tauv/(1+omega2**2*tauv**2))
c******************************************************************************
c******************** R1223zfs
************************************************
if ((dpara.gt.ws).or.(beta.gt.0.07)) then
cmp=zrr(3,3) !ZFS
czp=zrr(2,3)
cpp=zrr(1,3)
cmm=zrr(3,2)
czm=zrr(2,2)
cpm=zrr(1,2)
cmz=zrr(3,1)
czz=zrr(2,1)
cpz=zrr(1,1)
omega1=abs(wr(1)-wr(3)) !omold(5)) !ZFS
omega2=abs(wr(1)-wr(2)) !omold(3)) !ZFS
omega3=abs(wr(3)-wr(2)) !omold(2)) !ZFS
a1=cpp*conjg(cpz)-2*czp*conjg(czz)+cmp*conjg(cmz)
b1=cpp*conjg(czz)-czp*conjg(cmz)
c1=cmp*conjg(czz)-czp*conjg(cpz)
d1=cpp*conjg(cmz)
e1=cmp*conjg(cpz)
f1=cpm*conjg(cpz)-2*czm*conjg(czz)+cmm*conjg(cmz)
g1=cpm*conjg(czz)-czm*conjg(cmz)
h1=cmm*conjg(czz)-czm*conjg(cpz)
ai1=cpm*conjg(cmz)
al1=cmm*conjg(cpz)
t8=a1*f1/6.-(c1*g1+b1*h1)/2.+e1*ai1+d1*al1
r1223zfs=-t8*
& (tauv/(1+omega1**2*tauv**2)+
& tauv/(1+omega2**2*tauv**2))
endif
omega1=abs(wr(2)-wr(3)) !omold(5)) !ZFS
omega2=abs(wr(1)-wr(2)) !omold(3)) !ZFS
omega3=abs(wr(3)-wr(1)) !omold(2)) !ZFS
cmp=zrr(3,3)
czp=zrr(2,3)
cpp=zrr(1,3)
cmm=zrr(3,1) !solo Zeeman
czm=zrr(2,1)
cpm=zrr(1,1)
cmz=zrr(3,2)
czz=zrr(2,2)
cpz=zrr(1,2)
c******************************************************************************
! ++-0
a1=cpp*conjg(cpp)-2*czp*conjg(czp)+cmp*conjg(cmp)
b1=cpp*conjg(czp)-czp*conjg(cmp)
c1=conjg(czp)*cmp-czp*conjg(cpp)
d1=cpp*conjg(cmp)
e1=cmp*conjg(cpp)
f1=cpm*conjg(cpz)-2*czm*conjg(czz)+cmm*conjg(cmz)
g1=cpm*conjg(czz)-czm*conjg(cmz)
h1=cmm*conjg(czz)-czm*conjg(cpz)
ai1=cpm*conjg(cmz)
al1=cmm*conjg(cpz)
t51=a1*f1/6.-(c1*g1+b1*h1)/2.+e1*ai1+d1*al1
! +0-+
a1=cpp*conjg(cpz)-2*czp*conjg(czz)+cmp*conjg(cmz)
b1=cpp*conjg(czz)-czp*conjg(cmz)
c1=cmp*conjg(czz)-czp*conjg(cpz)
d1=cpp*conjg(cmz)
e1=cmp*conjg(cpz)
f1=cpp*conjg(cpm)-2*czp*conjg(czm)+cmp*conjg(cmm)
g1=cpp*conjg(czm)-czp*conjg(cmm)
h1=cmp*conjg(czm)-czp*conjg(cpm)
ai1=cpp*conjg(cmm)
al1=cmp*conjg(cpm)
t52=a1*f1/6.-(c1*g1+b1*h1)/2.+e1*ai1+d1*al1
! 00-0
a1=cpz*conjg(cpz)-2*czz*conjg(czz)+cmz*conjg(cmz)
b1=cpz*conjg(czz)-czz*conjg(cmz)
c1=cmz*conjg(czz)-czz*conjg(cpz)
d1=cmz*conjg(cpz)
e1=cpz*conjg(cmz)
f1=cpm*conjg(cpz)-2*czm*conjg(czz)+cmm*conjg(cmz)
g1=cpm*conjg(czz)-czm*conjg(cmz)
h1=cmm*conjg(czz)-czm*conjg(cpz)
ai1=cpm*conjg(cmz)
al1=cmm*conjg(cpz)
t53=a1*f1/6.-(c1*g1+b1*h1)/2.+e1*ai1+d1*al1
! -0--
a1=cpm*conjg(cpz)-2*czm*conjg(czz)+cmm*conjg(cmz)
b1=cpm*conjg(czz)-czm*conjg(cmz)
c1=cmm*conjg(czz)-czm*conjg(cpz)
d1=cpm*conjg(cmz)
e1=cmm*conjg(cpz)
f1=cpm*conjg(cpm)-2*czm*conjg(czm)+cmm*conjg(cmm)
g1=cpm*conjg(czm)-czm*conjg(cmm)
h1=cmm*conjg(czm)-czm*conjg(cpm)
ai1=cpm*conjg(cmm)
al1=cmm*conjg(cpm)
t54=a1*f1/6.-(c1*g1+b1*h1)/2.+e1*ai1+d1*al1
r1213=-(t51*tauv+t51*
& tauv/(1+omega2**2*tauv**2)
& -t52*
& tauv/(1+omega3**2*tauv**2)
& -t53*
& tauv/(1+omega2**2*tauv**2)
& -t54*tauv)
! 0-+-
a1=cpz*conjg(cpm)-2*czz*conjg(czm)+cmz*conjg(cmm)
b1=cpz*conjg(czm)-czz*conjg(cmm)
c1=cmz*conjg(czm)-czz*conjg(cpm)
d1=cpz*conjg(cmm)
e1=cmz*conjg(cpm)
f1=cpm*conjg(cpp)-2*czm*conjg(czp)+cmm*conjg(cmp)
g1=cpm*conjg(czp)-czm*conjg(cmp)
h1=cmm*conjg(czp)-czm*conjg(cpp)
ai1=cpm*conjg(cmp)
al1=cmm*conjg(cpp)
t41=a1*f1/6.-(c1*g1+b1*h1)/2.+e1*ai1+d1*al1
r2331=-t41*
& (tauv/(1+omega3**2*tauv**2)+
& tauv/(1+omega2**2*tauv**2))
rt2a=(rpzpz+rzmzm+sqrt((rpzpz-rzmzm)**2+4*r1223**2))/2
rt2b=(rpzpz+rzmzm-sqrt((rpzpz-rzmzm)**2+4*r1223**2))/2
!++--
a1=cpp*conjg(cpp)-2*czp*conjg(czp)+cmp*conjg(cmp)
b1=cpp*conjg(czp)-czp*conjg(cmp)
c1=conjg(czp)*cmp-czp*conjg(cpp)
d1=cpp*conjg(cmp)
e1=cmp*conjg(cpp)
f1=cpm*conjg(cpm)-2*czm*conjg(czm)+cmm*conjg(cmm)
g1=cpm*conjg(czm)-czm*conjg(cmm)
h1=cmm*conjg(czm)-czm*conjg(cpm)
ai1=cpm*conjg(cmm)
al1=cmm*conjg(cpm)
t21=a1*f1/6.-(c1*g1+b1*h1)/2.+e1*ai1+d1*al1
rpmpm=-(2*t21*tauv
& -t6*tauv
& -t3*
& tauv/(1+omega1**2*tauv**2)
& -2*t78*
& tauv/(1+omega3**2*tauv**2)
& -t5*
& tauv/(1+omega2**2*tauv**2)
& -t12*tauv)
C*****************************************************************************
C*********** CALCULATION OF ELECTRONIC R1
************************************
C*********** IN AXIAL ROUTINE
************************************
C*****************************************************************************
! originali
a1=cpz*cpp-2*czz*czp+cmz*cmp
b1=cpz*czp-czz*cmp
c1=cmz*czp-czz*cpp
d1=cpz*cmp
e1=cmz*cpp
f1=cpp*cpz-2*czp*czz+cmp*cmz
g1=cpp*czz-czp*cmz
h1=cmp*czz-czp*cpz
ai1=cpp*cmz
al1=cmp*cpz
trhoa=a1*f1/6.-(c1*g1+b1*h1)/2.+e1*ai1+d1*al1
rhoa=2*trhoa*
& tauv/(1+omega1**2*tauv**2)
a1=cpp*cpm-2*czp*czm+cmp*cmm
b1=cpp*czm-czp*cmm
c1=cmp*czm-czp*cpm
d1=cpp*cmm
e1=cmp*cpm
f1=cpm*cpp-2*czm*czp+cmm*cmp
g1=cpm*czp-czm*cmp
h1=cmm*czp-czm*cpp
ai1=cpm*cmp
al1=cmm*cpp
trhob=a1*f1/6.-(c1*g1+b1*h1)/2.+e1*ai1+d1*al1
rhob=2*trhob*
& tauv/(1+omega3**2*tauv**2)
a1=cpz*cpm-2*czz*czm+cmz*cmm
b1=cpz*czm-czz*cmm
c1=cmz*czm-czz*cpm
d1=cpz*cmm
e1=cmz*cpm
f1=cpm*cpz-2*czm*czz+cmm*cmz
g1=cpm*czz-czm*cmz
h1=cmm*czz-czm*cpz
ai1=cpm*cmz
al1=cmm*cpz
trhoc=a1*f1/6.-(c1*g1+b1*h1)/2.+e1*ai1+d1*al1
rhoc=2*trhoc*
& tauv/(1+omega2**2*tauv**2)
C******************************************************************************
C********** WRITING RELAXATION RATES IN AXIAL TOUTINE
*************************
C******************************************************************************
part=(rhoa+rhob+rhoc)
if(abs(rhoa/2.+rhoc/2.).gt.abs(rhob))then
rt1=(part-sqrt(part**2-3.*(rhoa*rhob+rhoa*rhoc+rhob*rhoc)))
else
rt1=(part+sqrt(part**2-3.*(rhoa*rhob+rhoa*rhoc+rhob*rhoc)))
endif
rt11=(part-sqrt(part**2-3.*(rhoa*rhob+rhoa*rhoc+rhob*rhoc)))
rt12=(part+sqrt(part**2-3.*(rhoa*rhob+rhoa*rhoc+rhob*rhoc)))
rt1=rt1*sqrt(abs(dpara-ws)/(dpara+ws))+(rt11+rt12)/2*
^ (1-sqrt(abs(dpara-ws)/(dpara+ws)))
bfield=bz*6.283/GAMMAI
if (beta.gt.1.55) then
endif
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
C*****************************************************************************
C*************** DEFINITION OF SUPERMATRIX
***********************************
C*************** IN AXIAL ROUTINE
***********************************
C*****************************************************************************
a=delta**2/5.*rhoa
b=delta**2/5.*rhoc
cp=delta**2/5.*rhob
d=1./taur
e=wi
! write(6,*)a,b,cp,d,e
a11=(d*(3*b*cp+3*a*(b+cp)+2*(a+b+cp)*d+d*d)*(d*(cp+d)+a*
& (b+cp+d)+b*(cp+2*d))+(4*b*b*(cp+d)+2*d*(cp+d)**2+
& a*a*(b+cp+2*d)+b*(cp*cp+6*cp*d+4*d*d)+a*(4*b*b+6*b*
& (cp+d)+(cp+2*d)**2))*e*e+(a+cp+d)*e**4)/
& ((d**2+e**2)*((3*b*cp+3*a*(b+cp)+2*(a+b+cp)*d+d*d)**2+
& 2*(2*a*a+a*b+2*b*b+a*cp+b*cp+2*cp*cp+2*(a+b+cp)*d+d*d)*
& e*e+e**4))
a12=(d*(3*b*cp+3*a*(b+cp)+2*(a+b+cp)*d+d*d)*(b*cp+a*
& (b+cp+d))+(a*a*(b+cp+2*d)-b*cp*(2*(b+cp)+3*d)-
& a*(2*b*b+cp*(2*cp+d)+b*(3*cp+d)))*e*e-a*e**4)/
& ((d**2+e**2)*((3*b*cp+3*a*(b+cp)+2*(a+b+cp)*d+d*d)**2+
& 2*(2*a*a+a*b+2*b*b+a*cp+b*cp+2*cp*cp+2*(a+b+cp)*d+d*d)*
& e*e+e**4))
a13=(d*(3*b*cp+3*a*(b+cp)+2*(a+b+cp)*d+d*d)*(a*(b+cp)+cp*(b+d))-
& (2*a*a*(b+cp)+cp*(2*b*b-2*cp*d+b*(-cp+d))+a*(2*b*b+cp*
& (-cp+d)+3*b*(cp+d)))*e*e-cp*e**4)/
& ((d**2+e**2)*((3*b*cp+3*a*(b+cp)+2*(a+b+cp)*d+d*d)**2+
& 2*(2*a*a+a*b+2*b*b+a*cp+b*cp+2*cp*cp+2*(a+b+cp)*d+d*d)*
& e*e+e**4))
a21=(d*(3*b*cp+3*a*(b+cp)+2*(a+b+cp)*d+d*d)*(b*cp+a*
& (b+cp+d))+(a*a*(b+cp+2*d)-b*cp*(2*(b+cp)+3*d)-
& a*(2*b*b+cp*(2*cp+d)+b*(3*cp+d)))*e*e-a*e**4)/
& ((d**2+e**2)*((3*b*cp+3*a*(b+cp)+2*(a+b+cp)*d+d*d)**2+
& 2*(2*a*a+a*b+2*b*b+a*cp+b*cp+2*cp*cp+2*(a+b+cp)*d+d*d)*
& e*e+e**4))
a22=(d*(3*b*cp+3*a*(b+cp)+2*(a+b+cp)*d+d*d)*(b*(cp+d)+a*
& (b+cp+d)+d*(2*cp+d))+(a*a*b+a*b*b+a*a*cp+6*a*b*cp+b*b*cp+
& 4*a*cp*cp+4*b*cp*cp+2*(a+b+cp)*(a+b+2*cp)*d+4*(a+b+cp)*
& d*d+2*d**3)*e*e+(a+b+d)*e**4)/
& ((d**2+e**2)*((3*b*cp+3*a*(b+cp)+2*(a+b+cp)*d+d*d)**2+
& 2*(2*a*a+a*b+2*b*b+a*cp+b*cp+2*cp*cp+2*(a+b+cp)*d+d*d)*
& e*e+e**4))
a23=(d*(3*b*cp+3*a*(b+cp)+2*(a+b+cp)*d+d*d)*(a*(b+cp)+b*(cp+d))-
& (2*a*a*(b+cp)+b*(cp*(2*cp+d)-b*(cp+2*d))+a*(-b*b+b*(3*cp
+d)
+
& cp*(2*cp+3*d)))*e*e-b*e**4)/
& ((d**2+e**2)*((3*b*cp+3*a*(b+cp)+2*(a+b+cp)*d+d*d)**2+
& 2*(2*a*a+a*b+2*b*b+a*cp+b*cp+2*cp*cp+2*(a+b+cp)*d+d*d)*
& e*e+e**4))
a31=(d*(3*b*cp+3*a*(b+cp)+2*(a+b+cp)*d+d*d)*(a*(b+cp)+cp*(b+d))-
& (2*a*a*(b+cp)+cp*(2*b*b-2*cp*d+b*(-cp+d))+a*(2*b*b+cp*
& (-cp+d)+3*b*(cp+d)))*e*e-cp*e**4)/
& ((d**2+e**2)*((3*b*cp+3*a*(b+cp)+2*(a+b+cp)*d+d*d)**2+
& 2*(2*a*a+a*b+2*b*b+a*cp+b*cp+2*cp*cp+2*(a+b+cp)*d+d*d)*
& e*e+e**4))
a32=(d*(3*b*cp+3*a*(b+cp)+2*(a+b+cp)*d+d*d)*(a*(b+cp)+b*(cp+d))-
& (2*a*a*(b+cp)+b*(cp*(2*cp+d)-b*(cp+2*d))+a*(-b*b+b*(3*cp
+d)
+
& cp*(2*cp+3*d)))*e*e-b*e**4)/
& ((d**2+e**2)*((3*b*cp+3*a*(b+cp)+2*(a+b+cp)*d+d*d)**2+
& 2*(2*a*a+a*b+2*b*b+a*cp+b*cp+2*cp*cp+2*(a+b+cp)*d+d*d)*
& e*e+e**4))
a33=(d*(3*b*cp+3*a*(b+cp)+2*(a+b+cp)*d+d*d)*((b+d)*(cp+d)+a*
& (b+cp+2*d))+(2*d*(cp+d)**2+4*a*a*(b+cp+d)+b*b*(cp+2*d)+b
& *(cp+2*d)**2+a*(b**2+cp**2+6*cp*d+4*d**2+6*b*(cp+d)))*
& e*e+(b+cp+d)*e**4)/
& ((d**2+e**2)*((3*b*cp+3*a*(b+cp)+2*(a+b+cp)*d+d*d)**2+
& 2*(2*a*a+a*b+2*b*b+a*cp+b*cp+2*cp*cp+2*(a+b+cp)*d+d*d)*
& e*e+e**4))
c*************** defined in Zeeman limit
***********************************
a=(delta**2/5.)*rpzpz
b=(delta**2/5.)*rzmzm
f=(delta**2/5.)*rpmpm
cp=(delta**2/5.)*r1223
d=1.0/taur
e=wi
ea=omega1
eb=omega2
cpzfs=(delta**2/5.)*r1223zfs
cpzee=(delta**2/5.)*r1223zee
c******************************************************************************
c**************** Check secular approximation for R1223
***********************
c******************************************************************************
if (dpara.gt.ws) then
w1zfs=omega3+e
w2zfs=omega2+e
w3zfs=omega1+e
a44r1=real(f)+d
a55r1=real(b)+d
a44zfs=dcmplx(a44r1,w1zfs)
a55zfs=dcmplx(a55r1,w2zfs)
cpzfs1=cdabs(cpzfs/(a44zfs-a55zfs))
else
w1zee=omega1
w2zee=omega2
w3zee=omega3
a44r2=real(a)+d
a55r2=real(b)+d
a44zee=dcmplx(a44r2,w1zee)
a55zee=dcmplx(a55r2,w2zee)
cpzee1=cdabs(cpzee/(a44zee-a55zee))
endif
c*********** COMMENT: I FOUND THE DC APPROXIMATION UNNECESSARY,
MEAMI.*********
c******************************************************************************
c******************************************************************************
C************ defined in zfs or Zeeman
****************************************
C******************************************************************************
if ((dpara.gt.ws).or.(beta.gt.0.07)) then
a=(delta**2/5.)*rpmpm
b=(delta**2/5.)*rzmzm
f=(delta**2/5.)*rpzpz
cp=cpzfs
d=1.0/taur
e=wi
ea=omega3
eb=omega2
ec=omega1
else
a=(delta**2/5.)*rpzpz
b=(delta**2/5.)*rzmzm
f=(delta**2/5.)*rpmpm
cp=cpzee
d=1.0/taur
e=wi
ea=omega1
eb=omega2
ec=omega3
endif
c*************************** writing out eigenvalues for R2e
*****************
rt2a=(a+b+sqrt((a-b)**2+4*cp**2))/2
rt2b=(a+b-sqrt((a-b)**2+4*cp**2))/2
rt2c=f
if (beta.lt.0.06) then
c WRITE(6,*)bfield,rt2a
c WRITE(6,*)bfield,rt2b
c WRITE(6,*)bfield,rt2c
endif
c******************************************************************************
a44=(b**2*d-b*(cp**2-2*d**2)+d*(-(cp**2)+d**2+(e+eb)**2)+
& a*((b+d)**2+(e+eb)**2))/(cp**4+2*cp**2*(-d*(b+d)+
& (e+ea)*(e+eb))+a**2*((b+d)**2+(e+eb)**2)+(d**2+(e+ea)**2)*
& ((b+d)**2+(e+eb)**2)+2*a*(b**2*d-b*(cp**2-2*d**2)+
& d*(-(cp**2)+d**2+(e+eb)**2)))
a55=(-(cp**2)*d+a**2*(b+d)+a*(-(cp**2)+2*d*(b+d))+(b+d)*
& (d**2+(e+ea)**2))/(cp**4+2*cp**2*(-d*(b+d)+
& (e+ea)*(e+eb))+a**2*((b+d)**2+(e+eb)**2)+(d**2+(e+ea)**2)*
& ((b+d)**2+(e+eb)**2)+2*a*(b**2*d-b*(cp**2-2*d**2)+
& d*(-(cp**2)+d**2+(e+eb)**2)))
a45=-(cp*(-(cp**2)+(a+d)*(b+d)-(e+ea)*(e+eb)))/
& ((cp**2-(a+d)*(b+d)+(e+ea)*(e+eb))**2+
& (b*(e+ea)+a*(e+eb)+d*(2*e+ea+eb))**2)
a54=a45
taus3=1./(d+f)
a66=taus3/(1+(wi+ec)**2*taus3**2)
a77=a44
a88=a55
a78=a45
a87=a78
a99=a66
a65=0.0d0
a56=0.0d0
a89=0.0d0
a98=0.0d0
if (dpara.gt.ws) then
astore1=a44
astore2=a66
a44=a99
a66=a77
a99=astore1
a77=astore2
a65=a45
a56=a65
a89=a65
a98=a65
a45=0.0d0
a54=0.0d0
a78=0.0d0
a87=0.0d0
endif
C******************************************************************************
C******** DEFINITION OF COEFF. MATRIX FOR DIAGONALIZED H0
*********************
C******** IN AXIAL ROUTINE
*********************
C******************************************************************************
zrn(1,1)=zr(1,3) !Zeeman
zrn(1,2)=zr(2,3)
zrn(1,3)=zr(3,3)
zrn(2,1)=zr(1,2)
zrn(2,2)=zr(2,2)
zrn(2,3)=zr(3,2)
zrn(3,1)=zr(1,1)
zrn(3,2)=zr(2,1)
zrn(3,3)=zr(3,1)
C******************************************************************************
C************** DIAGONALIZATION OF COEFF. MATRIX
******************************
C************** IN AXIAL ROUTINE
******************************
C******************************************************************************
lda=3
num=3
call zgefa(zrn,lda,num,ipvt,info)
job=01
call zgedi(zrn,lda,num,ipvt,det,work,job)
ak11=zrn(1,1)
ak12=zrn(1,2)
ak13=zrn(1,3)
ak21=zrn(2,1)
ak22=zrn(2,2)
ak23=zrn(2,3)
ak31=zrn(3,1)
ak32=zrn(3,2)
ak33=zrn(3,3)
C******************************************************************************
C************* LIOUVILLE SPACE REPRESENTATION OF S-OPERATORS
******************
C************* PROJECTION VECTORS IN AXIAL ROUTINE
******************
C******************************************************************************
cz1=ak11**2-ak31**2
cz2=ak12**2-ak32**2
cz3=ak13**2-ak33**2
cz4=ak11*ak12-ak31*ak32
cz5=ak12*ak13-ak32*ak33
cz6=ak11*ak13-ak31*ak33
cz7=ak12*ak11-ak32*ak31
cz8=ak13*ak12-ak33*ak32
cz9=ak13*ak11-ak33*ak31
cp1=-sqrt(2.)*(ak21*ak31+ak11*ak21)
cp2=-sqrt(2.)*(ak22*ak32+ak12*ak22)
cp3=-sqrt(2.)*(ak23*ak33+ak13*ak23)
cp4=-sqrt(2.)*(ak21*ak32+ak11*ak22)
cp5=-sqrt(2.)*(ak22*ak33+ak12*ak23)
cp6=-sqrt(2.)*(ak21*ak33+ak11*ak23)
cp7=-sqrt(2.)*(ak22*ak31+ak12*ak21)
cp8=-sqrt(2.)*(ak23*ak32+ak13*ak22)
cp9=-sqrt(2.)*(ak23*ak31+ak13*ak21)
cm1=sqrt(2.)*(ak31*ak21+ak21*ak11)
cm2=sqrt(2.)*(ak32*ak22+ak22*ak12)
cm3=sqrt(2.)*(ak33*ak23+ak23*ak13)
cm4=sqrt(2.)*(ak31*ak22+ak21*ak12)
cm5=sqrt(2.)*(ak32*ak23+ak22*ak13)
cm6=sqrt(2.)*(ak31*ak23+ak21*ak13)
cm7=sqrt(2.)*(ak32*ak21+ak22*ak11)
cm8=sqrt(2.)*(ak33*ak22+ak23*ak12)
cm9=sqrt(2.)*(ak33*ak21+ak23*ak11)
C******************************************************************************
C**************** PROJECTION OF SUPERMATRIX
***********************************
C**************** IN AXIAL ROUTINE
***********************************
C******************************************************************************
ammcz1=a11*1e9*cz1+a12*1e9*cz2+a13*1e9*cz3
ammcz2=a21*1e9*cz1+a22*1e9*cz2+a23*1e9*cz3
ammcz3=a31*1e9*cz1+a32*1e9*cz2+a33*1e9*cz3
ammcz4=a44*1e9*cz4+a45*1e9*cz5
ammcz5=a55*1e9*cz5+a54*1e9*cz4+a56*1e9*cz6
ammcz6=a66*1e9*cz6+a65*1e9*cz5
ammcz7=a77*1e9*cz7+a78*1e9*cz8
ammcz8=a88*1e9*cz8+a87*1e9*cz7+a89*1e9*cz9
ammcz9=a99*1e9*cz9+a98*1e9*cz8
ammcp1=a11*1e9*cp1+a12*1e9*cp2+a13*1e9*cp3
ammcp2=a21*1e9*cp1+a22*1e9*cp2+a23*1e9*cp3
ammcp3=a31*1e9*cp1+a32*1e9*cp2+a33*1e9*cp3
ammcp4=a44*1e9*cp4+a45*1e9*cp5
ammcp5=a55*1e9*cp5+a54*1e9*cp4+a56*1e9*cp6
ammcp6=a66*1e9*cp6+a65*1e9*cp5
ammcp7=a77*1e9*cp7+a78*1e9*cp8
ammcp8=a88*1e9*cp8+a87*1e9*cp7+a89*1e9*cp9
ammcp9=a99*1e9*cp9+a98*1e9*cp8
ammcm1=a11*1e9*cm1+a12*1e9*cm2+a13*1e9*cm3
ammcm2=a21*1e9*cm1+a22*1e9*cm2+a23*1e9*cm3
ammcm3=a31*1e9*cm1+a32*1e9*cm2+a33*1e9*cm3
ammcm4=a44*1e9*cm4+a45*1e9*cm5
ammcm5=a55*1e9*cm5+a54*1e9*cm4+a56*1e9*cm6
ammcm6=a66*1e9*cm6+a65*1e9*cm5
ammcm7=a77*1e9*cm7+a78*1e9*cm8
ammcm8=a88*1e9*cm8+a87*1e9*cm7+a89*1e9*cm9
ammcm9=a99*1e9*cm9+a98*1e9*cm8
C******************************************************************************
C**************** SPECTRAL DENSITIES OF S
*************************************
C**************** IN AXIAL ROUTINE
*************************************
C******************************************************************************
!cz M cz
primo=abs(cz1*ammcz1+cz2*ammcz2+cz3*ammcz3)+abs(cz4*ammcz4+
& cz5*ammcz5+cz6*ammcz6+cz7*ammcz7+cz8*ammcz8+cz9*ammcz9)
!cp M cp
secondo=abs(cp1*ammcp1+cp2*ammcp2+cp3*ammcp3)+abs(cp4*ammcp4+
& cp5*ammcp5+cp6*ammcp6+cp7*ammcp7+cp8*ammcp8+cp9*ammcp9)
!cp M cm
terzo=abs(cp1*ammcm1+cp2*ammcm2+cp3*ammcm3)+abs(cp4*ammcm4+
& cp5*ammcm5+cp6*ammcm6+cp7*ammcm7+cp8*ammcm8+cp9*ammcm9)
!cp M cz
quarto=abs(cp1*ammcz1+cp2*ammcz2+cp3*ammcz3)+abs(cp4*ammcz4+
& cp5*ammcz5+cp6*ammcz6+cp7*ammcz7+cp8*ammcz8+cp9*ammcz9)
COLD(1,1)=gz**2*abs(secondo)
COLD(1,2)=gz**2*abs(quarto)
COLD(1,3)=gz**2*abs(terzo)
COLD(1,4)=gz**2*abs(primo)
!cm M cp
quinto=cm1*ammcp1+cm2*ammcp2+cm3*ammcp3+cm4*ammcp4+cm5*ammcp5+
& cm6*ammcp6+cm7*ammcp7+cm8*ammcp8+cm9*ammcp9
!cm M cm
sesto=cm1*ammcm1+cm2*ammcm2+cm3*ammcm3+cm4*ammcm4+cm5*ammcm5+
& cm6*ammcm6+cm7*ammcm7+cm8*ammcm8+cm9*ammcm9
!cm M cz
settimo=cm1*ammcz1+cm2*ammcz2+cm3*ammcz3+cm4*ammcz4+cm5*ammcz5+
& cm6*ammcz6+cm7*ammcz7+cm8*ammcz8+cm9*ammcz9
!cz M cp
ottavo=cz1*ammcp1+cz2*ammcp2+cz3*ammcp3+cz4*ammcp4+cz5*ammcp5+
& cz6*ammcp6+cz7*ammcp7+cz8*ammcp8+cz9*ammcp9
!cz M cm
onono=cz1*ammcm1+cz2*ammcm2+cz3*ammcm3+cz4*ammcm4+cz5*ammcm5+
& cz6*ammcm6+cz7*ammcm7+cz8*ammcm8+cz9*ammcm9
C******************************************************************************
C**************** SPECTRAL DENSITIES OF G
*************************************
C**************** IN AXIAL ROUTINE
*************************************
C******************************************************************************
C AUTOCORRELATION FUNCTIONS
C(1,1,11)=secondo
C(1,1,12)=sesto
C(1,1,13)=quinto
C(1,1,14)= terzo
C(1,1,15)=quarto
C(1,1,16)= ottavo
C(1,1,17)=settimo
C(1,1,18)=onono
C(1,1,19)=primo
540 CONTINUE
GO TO 567
ENDIF
456 CONTINUE
C GENERAL MATRIX OF ENERGY
X=BZ*3.1415927*658.2
ZC=X*CT
ZS=X*ST
DO I=1,(2.*SI+1)*(2.*SPIN+1)*2
DO J=1,(2.*SI+1)*(2.*SPIN+1)*2
AR(I,J)=0.
END DO
END DO
SSI = SI*(SI + 1.)
SS=SPIN*(SPIN+1.)
ISMS=2.*SPIN+1
K=0
DO I=1,(2.*SI+1)*(2.*SPIN+1),2.*SPIN+1
DO J=0,2.*SPIN
COEF = SI - K
COEF2 = SPIN-J
AR(I+J,I+J)=2.*(SPIN-J)*ZC*GZ + (SPIN-J)*AZ*(SI-K)+DPARA*
$ (COEF2**2-SS/3.)+S4/120.*(35.*COEF2**4-30.*SS*COEF2**2+
$ 25.*COEF2**2-6.*SS+3.*SS**2)
IF (J+1.LT.2.*SPIN+1) THEN
AR(I+J,I+J+1)=CMPLX(COEFFH*ZS*GX*COS(GAMMA)*SQRT(SS-
$ (COEF2-1.)*COEF2),-ZS*GY*SIN(GAMMA)*SQRT(SS-(COEF2-1.)*COEF2))
AR(I+J+1,I+J)=CMPLX(COEFFH*ZS*GX*COS(GAMMA)*SQRT(SS-(COEF2-1.)*
$ COEF2),ZS*GY*SIN(GAMMA)*SQRT(SS-(COEF2-1.)*COEF2))
IF (J+2.LT.2.*SPIN+1) THEN
AR(I+J,I+J+2)=EPARA*
& SQRT(SS-COEF2*(COEF2-1.))*
& SQRT(SS-(COEF2-1.)*(COEF2-2.))/2.
AR(I+J+2,I+J)=AR(I+J,I+J+2)
ENDIF
IF (J+4.LT.2.*SPIN+1) THEN
AR(I+J,I+J+4)=S4/48.*
& SQRT(SS-COEF2*(COEF2-1.))*
& SQRT(SS-(COEF2-1.)*(COEF2-2.))*
& SQRT(SS-(COEF2-2.)*(COEF2-3.))*
& SQRT(SS-(COEF2-3.)*(COEF2-4.))
AR(I+J+4,I+J)=AR(I+J,I+J+4)
ENDIF
IF (I+(2.*SPIN+1)+J.LE.(2.*SPIN+1)*(2.*SI+1))THEN
AR(I+(ISMS)+J,(I+1+J))=.5*(AX+AY)/2.*SQRT(SSI-(COEF-1.)*COEF)*
$ SQRT(SS-(COEF2-1.)*COEF2)
AR((I+1+J),I+(ISMS)+J)= AR(I+(ISMS)+J,(I+1+J))
ENDIF
IF (I+(2.*SPIN+1)+J+1.LE.(2.*SPIN+1)*(2.*SI+1))THEN
AR(I+(ISMS)+J+1,(I+J))=.5*(AX-AY)/2.*SQRT(SSI-(COEF-1.)*COEF)*
$ SQRT(SS-(COEF2-1.)*COEF2)
AR((I+J),I+(ISMS)+J+1)= AR(I+(ISMS)+J+1,(I+J))
ENDIF
ENDIF
END DO
K=K+1
END DO
IF (INDEX.EQ.1)THEN
WRITE(6,*)'DIM. MATRIX', NMX
OPEN(UNIT=17,FILE='MAT')
DO I=1,(2.*SI+1)*(2.*SPIN+1)
DO J=1,(2.*SI+1)*(2.*SPIN+1)
WRITE(17,*)AR(I,J)
END DO
WRITE(17,*)' '
END DO
CLOSE(17)
ENDIF
INDEX=INDEX+1
C WR EIGENVALUES ZR EIGENVECTORS
DO 245 I=1,NMX
DO 245 J=1,NMX
ARR(I,J)=REAL(AR(I,J))
ARI(I,J)=IMAG(AR(I,J))
245 CONTINUE
CALL F02AXF(ARR,MBRANC,ARI,MBRANC,NMX,WR,ZRR,MBRANC,ZRI,MBRANC
$ ,WK1,WK2,WK3,0)
DO 246 I=1,NMX
DO 246 J=1,NMX
ZR(I,J)=CMPLX(ZRR(I,J),ZRI(I,J))
! write(70,*)i,j,zr(i,j),wr(j) !scrivi
246 CONTINUE
I=1
OM(1,1)=0.
OMOLD(1)=0.
DO 70 K=1,NMX
DO 70 L=1,NMX
IF (K.EQ.L)GO TO 70
I=I+1
OM(K,L)=WR(K)-WR(L)
OMOLD(I)=WR(K)-WR(L)
70 CONTINUE
J=0
DO I=1,(2.*SI+1)*(2.*SPIN+1)
J=J+1
IF (J.GT.(2*SPIN+1)) J=1
CO1(I)=(2*SPIN+1-(2*J-1))/2.
C CO1(I)=-((-1.)**J)
END DO
DO I=1,(2.*SI+1)*(2.*SPIN+1)-1
CO2(I)=SQRT(SPIN*(SPIN+1)-CO1(I+1)*(CO1(I+1)+1.))
END DO
DO I=2,(2.*SI+1)*(2.*SPIN+1)
CO3(I)=SQRT(SPIN*(SPIN+1)-CO1(I-1)*(CO1(I-1)-1.))
END DO
IF (INDEX.EQ.2)THEN
OPEN(UNIT=18,FILE='COE')
DO I=1,NMX
WRITE(18,*)CO1(I),CO2(I), CO3(I)
END DO
CLOSE(18)
ENDIF
C******************************************************************************
C******************************************************************************
C*************** RHOMBIC ROUTINE
*********************************************
C******************************************************************************
C******************************************************************************
!
C*****************************************************************************
!C************* DEFINITION OF COEFF. MATRIX
***********************************
!C************* IN RHOMBIC ROUTINE
***********************************
!
C*****************************************************************************
omega1=abs(wr(2)-wr(3)) !omold(5)) !ZFS
omega2=abs(wr(1)-wr(2)) !omold(3)) !ZFS
omega3=abs(wr(3)-wr(1)) !omold(2)) !ZFS
cmp=zr(3,3)
czp=zr(2,3)
cpp=zr(1,3)
cmm=zr(3,1) !solo Zeeman
czm=zr(2,1)
cpm=zr(1,1)
cmz=zr(3,2)
czz=zr(2,2)
cpz=zr(1,2)
C*****************************************************************************
C*********** CACLULATION OF ELECTRONIC R1
************************************
C*********** IN RHOMBIC ROUTINE
************************************
C*****************************************************************************
!calcolo R1
a1=cpz*conjg(cpp)-2*czz*conjg(czp)+cmz*conjg(cmp)
b1=cpz*conjg(czp)-czz*conjg(cmp)
c1=cmz*conjg(czp)-czz*conjg(cpp)
d1=cpz*conjg(cmp)
e1=cmz*conjg(cpp)
f1=cpp*conjg(cpz)-2*czp*conjg(czz)+cmp*conjg(cmz)
g1=cpp*conjg(czz)-czp*conjg(cmz)
h1=cmp*conjg(czz)-czp*conjg(cpz)
ai1=cpp*conjg(cmz)
al1=cmp*conjg(cpz)
trhoa=a1*f1/6.-(c1*g1+b1*h1)/2.+e1*ai1+d1*al1
rhoa=2*trhoa*
& tauv/(1+omega1**2*tauv**2)
a1=cpp*conjg(cpm)-2*czp*conjg(czm)+cmp*conjg(cmm)
b1=cpp*conjg(czm)-czp*conjg(cmm)
c1=cmp*conjg(czm)-czp*conjg(cpm)
d1=cpp*conjg(cmm)
e1=cmp*conjg(cpm)
f1=cpm*conjg(cpp)-2*czm*conjg(czp)+cmm*conjg(cmp)
g1=cpm*conjg(czp)-czm*conjg(cmp)
h1=cmm*conjg(czp)-czm*conjg(cpp)
ai1=cpm*conjg(cmp)
al1=cmm*conjg(cpp)
trhob=a1*f1/6.-(c1*g1+b1*h1)/2.+e1*ai1+d1*al1
rhob=2*trhob*
& tauv/(1+omega3**2*tauv**2)
a1=cpz*conjg(cpm)-2*czz*conjg(czm)+cmz*conjg(cmm)
b1=cpz*conjg(czm)-czz*conjg(cmm)
c1=cmz*conjg(czm)-czz*conjg(cpm)
d1=cpz*conjg(cmm)
e1=cmz*conjg(cpm)
f1=cpm*conjg(cpz)-2*czm*conjg(czz)+cmm*conjg(cmz)
g1=cpm*conjg(czz)-czm*conjg(cmz)
h1=cmm*conjg(czz)-czm*conjg(cpz)
ai1=cpm*conjg(cmz)
al1=cmm*conjg(cpz)
trhoc=a1*f1/6.-(c1*g1+b1*h1)/2.+e1*ai1+d1*al1
rhoc=2*trhoc*
& tauv/(1+omega2**2*tauv**2)
part=(rhoa+rhob+rhoc)
if(abs(rhoa/2.+rhoc/2.).gt.abs(rhob))then
rt1=(part-sqrt(part**2-3.*(rhoa*rhob+rhoa*rhoc+rhob*rhoc)))
else
rt1=(part+sqrt(part**2-3.*(rhoa*rhob+rhoa*rhoc+rhob*rhoc)))
endif
rt11=(part-sqrt(part**2-3.*(rhoa*rhob+rhoa*rhoc+rhob*rhoc)))
rt12=(part+sqrt(part**2-3.*(rhoa*rhob+rhoa*rhoc+rhob*rhoc)))
rt1=rt1*sqrt(abs(dpara-ws)/(dpara+ws))+(rt11+rt12)/2*
^ (1-sqrt(abs(dpara-ws)/(dpara
+ws))) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
c write(6,*)'rhoa=',rhoa !cancellare
c write(6,*)'rhob=',rhob !cancellare
c write(6,*)'rhoc=',rhoc !cancellare
!
C******************************************************************************
C**************** CALCULATION OF ELECTRONIC R2
********************************
C**************** IN RHOMBIC ROUTINE
********************************
C******************************************************************************
!nuovo calcolo R2
!++00
a1=cpp*conjg(cpp)-2*czp*conjg(czp)+cmp*conjg(cmp)
b1=cpp*conjg(czp)-czp*conjg(cmp)
c1=conjg(czp)*cmp-czp*conjg(cpp)
d1=cpp*conjg(cmp)
e1=cmp*conjg(cpp)
f1=cpz*conjg(cpz)-2*czz*conjg(czz)+cmz*conjg(cmz)
g1=cpz*conjg(czz)-czz*conjg(cmz)
h1=cmz*conjg(czz)-czz*conjg(cpz)
ai1=cmz*conjg(cpz)
al1=cpz*conjg(cmz)
t10=a1*f1/6.-(c1*g1+b1*h1)/2.+e1*ai1+d1*al1
!0++0
a1=cpz*conjg(cpp)-2*czz*conjg(czp)+cmz*conjg(cmp)
b1=cpz*conjg(czp)-czz*conjg(cmp)
c1=cmz*conjg(czp)-czz*conjg(cpp)
d1=cpz*conjg(cmp)
e1=cmz*conjg(cpp)
f1=cpp*conjg(cpz)-2*czp*conjg(czz)+cmp*conjg(cmz)
g1=cpp*conjg(czz)-czp*conjg(cmz)
h1=cmp*conjg(czz)-czp*conjg(cpz)
ai1=cpp*conjg(cmz)
al1=cmp*conjg(cpz)
t3=a1*f1/6.-(c1*g1+b1*h1)/2.+e1*ai1+d1*al1
!0000
a1=cpz*conjg(cpz)-2*czz*conjg(czz)+cmz*conjg(cmz)
b1=cpz*conjg(czz)-czz*conjg(cmz)
c1=cmz*conjg(czz)-czz*conjg(cpz)
d1=cmz*conjg(cpz)
e1=cpz*conjg(cmz)
f1=cpz*conjg(cpz)-2*czz*conjg(czz)+cmz*conjg(cmz)
g1=cpz*conjg(czz)-czz*conjg(cmz)
h1=cmz*conjg(czz)-czz*conjg(cpz)
ai1=cmz*conjg(cpz)
al1=cpz*conjg(cmz)
t4=a1*f1/6.-(c1*g1+b1*h1)/2.+e1*ai1+d1*al1
!0--0
a1=cpz*conjg(cpm)-2*czz*conjg(czm)+cmz*conjg(cmm)
b1=cpz*conjg(czm)-czz*conjg(cmm)
c1=cmz*conjg(czm)-czz*conjg(cpm)
d1=cpz*conjg(cmm)
e1=cmz*conjg(cpm)
f1=cpm*conjg(cpz)-2*czm*conjg(czz)+cmm*conjg(cmz)
g1=cpm*conjg(czz)-czm*conjg(cmz)
h1=cmm*conjg(czz)-czm*conjg(cpz)
ai1=cpm*conjg(cmz)
al1=cmm*conjg(cpz)
t5=a1*f1/6.-(c1*g1+b1*h1)/2.+e1*ai1+d1*al1
!++++
a1=cpp*conjg(cpp)-2*czp*conjg(czp)+cmp*conjg(cmp)
b1=cpp*conjg(czp)-czp*conjg(cmp)
c1=conjg(czp)*cmp-czp*conjg(cpp)
d1=cpp*conjg(cmp)
e1=cmp*conjg(cpp)
f1=cpp*conjg(cpp)-2*czp*conjg(czp)+cmp*conjg(cmp)
g1=cpp*conjg(czp)-czp*conjg(cmp)
h1=conjg(czp)*cmp-czp*conjg(cpp)
ai1=cpp*conjg(cmp)
al1=cmp*conjg(cpp)
t6=a1*f1/6.-(c1*g1+b1*h1)/2.+e1*ai1+d1*al1
!-++-
a1=cpp*conjg(cpm)-2*czp*conjg(czm)+cmp*conjg(cmm)
b1=cpp*conjg(czm)-czp*conjg(cmm)
c1=cmp*conjg(czm)-czp*conjg(cpm)
d1=cpp*conjg(cmm)
e1=cmp*conjg(cpm)
f1=cpm*conjg(cpp)-2*czm*conjg(czp)+cmm*conjg(cmp)
g1=cpm*conjg(czp)-czm*conjg(cmp)
h1=cmm*conjg(czp)-czm*conjg(cpp)
ai1=cpm*conjg(cmp)
al1=cmm*conjg(cpp)
t78=a1*f1/6.-(c1*g1+b1*h1)/2.+e1*ai1+d1*al1
rpzpz=-(2*t10*tauv
& -t6*tauv
& -2*t3*
& tauv/(1+omega1**2*tauv**2)
& -t78*
& tauv/(1+omega3**2*tauv**2)
& -t4*tauv
& -t5*
& tauv/(1+omega2**2*tauv**2))
!----
a1=cpm*conjg(cpm)-2*czm*conjg(czm)+cmm*conjg(cmm)
b1=cpm*conjg(czm)-czm*conjg(cmm)
c1=cmm*conjg(czm)-czm*conjg(cpm)
d1=cpm*conjg(cmm)
e1=cmm*conjg(cpm)
f1=cpm*conjg(cpm)-2*czm*conjg(czm)+cmm*conjg(cmm)
g1=cpm*conjg(czm)-czm*conjg(cmm)
h1=cmm*conjg(czm)-czm*conjg(cpm)
ai1=cpm*conjg(cmm)
al1=cmm*conjg(cpm)
t12=a1*f1/6.-(c1*g1+b1*h1)/2.+e1*ai1+d1*al1
!00--
a1=cpz*conjg(cpz)-2*czz*conjg(czz)+cmz*conjg(cmz)
b1=cpz*conjg(czz)-czz*conjg(cmz)
c1=cmz*conjg(czz)-czz*conjg(cpz)
d1=cmz*conjg(cpz)
e1=cpz*conjg(cmz)
f1=cpm*conjg(cpm)-2*czm*conjg(czm)+cmm*conjg(cmm)
g1=cpm*conjg(czm)-czm*conjg(cmm)
h1=cmm*conjg(czm)-czm*conjg(cpm)
ai1=cpm*conjg(cmm)
al1=cmm*conjg(cpm)
t11=a1*f1/6.-(c1*g1+b1*h1)/2.+e1*ai1+d1*al1
rzmzm=-(2*t11*tauv
& -t3*
& tauv/(1+omega1**2*tauv**2)
& -t4*tauv
& -2*t5*
& tauv/(1+omega2**2*tauv**2)
& -t78*
& tauv/(1+omega3**2*tauv**2)
& -t12*tauv)
a1=cpp*conjg(cpz)-2*czp*conjg(czz)+cmp*conjg(cmz)
b1=cpp*conjg(czz)-czp*conjg(cmz)
c1=cmp*conjg(czz)-czp*conjg(cpz)
d1=cpp*conjg(cmz)
e1=cmp*conjg(cpz)
f1=cpm*conjg(cpz)-2*czm*conjg(czz)+cmm*conjg(cmz)
g1=cpm*conjg(czz)-czm*conjg(cmz)
h1=cmm*conjg(czz)-czm*conjg(cpz)
ai1=cpm*conjg(cmz)
al1=cmm*conjg(cpz)
t8=a1*f1/6.-(c1*g1+b1*h1)/2.+e1*ai1+d1*al1
c r1223=-t8*
c & (tauv/(1+omega1**2*tauv**2)+
c & tauv/(1+omega2**2*tauv**2))
r1223zee=-t8*
& (tauv/(1+omega1**2*tauv**2)+
& tauv/(1+omega2**2*tauv**2))
c******************************************************************************
c******************** R1223zfs
************************************************
if ((dpara.gt.ws))then
cmp=zr(3,3) !ZFS
czp=zr(2,3)
cpp=zr(1,3)
cmm=zr(3,2)
czm=zr(2,2)
cpm=zr(1,2)
cmz=zr(3,1)
czz=zr(2,1)
cpz=zr(1,1)
omega1=abs(wr(1)-wr(3)) !omold(5)) !ZFS
omega2=abs(wr(1)-wr(2)) !omold(3)) !ZFS
omega3=abs(wr(3)-wr(2)) !omold(2)) !ZFS
a1=cpp*conjg(cpz)-2*czp*conjg(czz)+cmp*conjg(cmz)
b1=cpp*conjg(czz)-czp*conjg(cmz)
c1=cmp*conjg(czz)-czp*conjg(cpz)
d1=cpp*conjg(cmz)
e1=cmp*conjg(cpz)
f1=cpm*conjg(cpz)-2*czm*conjg(czz)+cmm*conjg(cmz)
g1=cpm*conjg(czz)-czm*conjg(cmz)
h1=cmm*conjg(czz)-czm*conjg(cpz)
ai1=cpm*conjg(cmz)
al1=cmm*conjg(cpz)
t8=a1*f1/6.-(c1*g1+b1*h1)/2.+e1*ai1+d1*al1
r1223zfs=-t8*
& (tauv/(1+omega1**2*tauv**2)+
& tauv/(1+omega2**2*tauv**2))
endif
omega1=abs(wr(2)-wr(3)) !omold(5)) !ZFS
omega2=abs(wr(1)-wr(2)) !omold(3)) !ZFS
omega3=abs(wr(3)-wr(1)) !omold(2)) !ZFS
cmp=zr(3,3)
czp=zr(2,3)
cpp=zr(1,3)
cmm=zr(3,1) !solo Zeeman
czm=zr(2,1)
cpm=zr(1,1)
cmz=zr(3,2)
czz=zr(2,2)
cpz=zr(1,2)
rt2a=(rpzpz+rzmzm+sqrt((rpzpz-rzmzm)**2+4*r1223**2))/2
rt2b=(rpzpz+rzmzm-sqrt((rpzpz-rzmzm)**2+4*r1223**2))/2
!++--
a1=cpp*conjg(cpp)-2*czp*conjg(czp)+cmp*conjg(cmp)
b1=cpp*conjg(czp)-czp*conjg(cmp)
c1=conjg(czp)*cmp-czp*conjg(cpp)
d1=cpp*conjg(cmp)
e1=cmp*conjg(cpp)
f1=cpm*conjg(cpm)-2*czm*conjg(czm)+cmm*conjg(cmm)
g1=cpm*conjg(czm)-czm*conjg(cmm)
h1=cmm*conjg(czm)-czm*conjg(cpm)
ai1=cpm*conjg(cmm)
al1=cmm*conjg(cpm)
t21=a1*f1/6.-(c1*g1+b1*h1)/2.+e1*ai1+d1*al1
rpmpm=-(2*t21*tauv
& -t6*tauv
& -t3*
& tauv/(1+omega1**2*tauv**2)
& -2*t78*
& tauv/(1+omega3**2*tauv**2)
& -t5*
& tauv/(1+omega2**2*tauv**2)
& -t12*tauv)
rhombf=2*t78
C****************************************************************************
C***************** WRITING RELAXATION RATES IN RHOMBIC ROUTINE
**************
C****************************************************************************
c WRITE(6,*)'RHOMBIC ROUTINE'
c WRITE(6,*)'rpzpz=',rpzpz
c WRITE(6,*)'rzmzm=',rzmzm
c WRITE(6,*)'rpmpm=',rpmpm
c WRITE(6,*)'rhoa=',rhoa
c WRITE(6,*)'rhob=',rhob
c WRITE(6,*)'rhoc=',rhoc
C STOP
C*****************************************************************************
C*************** DEFINITION OF SUPERMATRIX
***********************************
C*************** IN RHOMBIC ROUTINE
***********************************
C*****************************************************************************
!!!! DA DEFINIRE CON COMPLESSI CONIUGATI E DA CONSIDERARE MATRICE
DI DIMENSIONE DOPPIA !!!!
a=delta**2/5*rhoa
b=delta**2/5*rhoc
cp=delta**2/5*rhob
d=1./taur
e=wi
a11=(d*(3*b*cp+3*a*(b+cp)+2*(a+b+cp)*d+d*d)*(d*(cp+d)+a*
& (b+cp+d)+b*(cp+2*d))+(4*b*b*(cp+d)+2*d*(cp+d)**2+
& a*a*(b+cp+2*d)+b*(cp*cp+6*cp*d+4*d*d)+a*(4*b*b+6*b*
& (cp+d)+(cp+2*d)**2))*e*e+(a+cp+d)*e**4)/
& ((d**2+e**2)*((3*b*cp+3*a*(b+cp)+2*(a+b+cp)*d+d*d)**2+
& 2*(2*a*a+a*b+2*b*b+a*cp+b*cp+2*cp*cp+2*(a+b+cp)*d+d*d)*
& e*e+e**4))
a12=(d*(3*b*cp+3*a*(b+cp)+2*(a+b+cp)*d+d*d)*(b*cp+a*
& (b+cp+d))+(a*a*(b+cp+2*d)-b*cp*(2*(b+cp)+3*d)-
& a*(2*b*b+cp*(2*cp+d)+b*(3*cp+d)))*e*e-a*e**4)/
& ((d**2+e**2)*((3*b*cp+3*a*(b+cp)+2*(a+b+cp)*d+d*d)**2+
& 2*(2*a*a+a*b+2*b*b+a*cp+b*cp+2*cp*cp+2*(a+b+cp)*d+d*d)*
& e*e+e**4))
a13=(d*(3*b*cp+3*a*(b+cp)+2*(a+b+cp)*d+d*d)*(a*(b+cp)+cp*(b+d))-
& (2*a*a*(b+cp)+cp*(2*b*b-2*cp*d+b*(-cp+d))+a*(2*b*b+cp*
& (-cp+d)+3*b*(cp+d)))*e*e-cp*e**4)/
& ((d**2+e**2)*((3*b*cp+3*a*(b+cp)+2*(a+b+cp)*d+d*d)**2+
& 2*(2*a*a+a*b+2*b*b+a*cp+b*cp+2*cp*cp+2*(a+b+cp)*d+d*d)*
& e*e+e**4))
a21=(d*(3*b*cp+3*a*(b+cp)+2*(a+b+cp)*d+d*d)*(b*cp+a*
& (b+cp+d))+(a*a*(b+cp+2*d)-b*cp*(2*(b+cp)+3*d)-
& a*(2*b*b+cp*(2*cp+d)+b*(3*cp+d)))*e*e-a*e**4)/
& ((d**2+e**2)*((3*b*cp+3*a*(b+cp)+2*(a+b+cp)*d+d*d)**2+
& 2*(2*a*a+a*b+2*b*b+a*cp+b*cp+2*cp*cp+2*(a+b+cp)*d+d*d)*
& e*e+e**4))
a22=(d*(3*b*cp+3*a*(b+cp)+2*(a+b+cp)*d+d*d)*(b*(cp+d)+a*
& (b+cp+d)+d*(2*cp+d))+(a*a*b+a*b*b+a*a*cp+6*a*b*cp+b*b*cp+
& 4*a*cp*cp+4*b*cp*cp+2*(a+b+cp)*(a+b+2*cp)*d+4*(a+b+cp)*
& d*d+2*d**3)*e*e+(a+b+d)*e**4)/
& ((d**2+e**2)*((3*b*cp+3*a*(b+cp)+2*(a+b+cp)*d+d*d)**2+
& 2*(2*a*a+a*b+2*b*b+a*cp+b*cp+2*cp*cp+2*(a+b+cp)*d+d*d)*
& e*e+e**4))
a23=(d*(3*b*cp+3*a*(b+cp)+2*(a+b+cp)*d+d*d)*(a*(b+cp)+b*(cp+d))-
& (2*a*a*(b+cp)+b*(cp*(2*cp+d)-b*(cp+2*d))+a*(-b*b+b*(3*cp
+d)
+
& cp*(2*cp+3*d)))*e*e-b*e**4)/
& ((d**2+e**2)*((3*b*cp+3*a*(b+cp)+2*(a+b+cp)*d+d*d)**2+
& 2*(2*a*a+a*b+2*b*b+a*cp+b*cp+2*cp*cp+2*(a+b+cp)*d+d*d)*
& e*e+e**4))
a31=(d*(3*b*cp+3*a*(b+cp)+2*(a+b+cp)*d+d*d)*(a*(b+cp)+cp*(b+d))-
& (2*a*a*(b+cp)+cp*(2*b*b-2*cp*d+b*(-cp+d))+a*(2*b*b+cp*
& (-cp+d)+3*b*(cp+d)))*e*e-cp*e**4)/
& ((d**2+e**2)*((3*b*cp+3*a*(b+cp)+2*(a+b+cp)*d+d*d)**2+
& 2*(2*a*a+a*b+2*b*b+a*cp+b*cp+2*cp*cp+2*(a+b+cp)*d+d*d)*
& e*e+e**4))
a32=(d*(3*b*cp+3*a*(b+cp)+2*(a+b+cp)*d+d*d)*(a*(b+cp)+b*(cp+d))-
& (2*a*a*(b+cp)+b*(cp*(2*cp+d)-b*(cp+2*d))+a*(-b*b+b*(3*cp
+d)
+
& cp*(2*cp+3*d)))*e*e-b*e**4)/
& ((d**2+e**2)*((3*b*cp+3*a*(b+cp)+2*(a+b+cp)*d+d*d)**2+
& 2*(2*a*a+a*b+2*b*b+a*cp+b*cp+2*cp*cp+2*(a+b+cp)*d+d*d)*
& e*e+e**4))
a33=(d*(3*b*cp+3*a*(b+cp)+2*(a+b+cp)*d+d*d)*((b+d)*(cp+d)+a*
& (b+cp+2*d))+(2*d*(cp+d)**2+4*a*a*(b+cp+d)+b*b*(cp+2*d)+b
& *(cp+2*d)**2+a*(b**2+cp**2+6*cp*d+4*d**2+6*b*(cp+d)))*
& e*e+(b+cp+d)*e**4)/
& ((d**2+e**2)*((3*b*cp+3*a*(b+cp)+2*(a+b+cp)*d+d*d)**2+
& 2*(2*a*a+a*b+2*b*b+a*cp+b*cp+2*cp*cp+2*(a+b+cp)*d+d*d)*
& e*e+e**4))
cpzfs=(delta**2/5.)*r1223zfs
cpzee=(delta**2/5.)*r1223zee
if ((dpara.gt.ws).or.(beta.gt.0.07)) then
a=(delta**2/5.)*rpmpm
b=(delta**2/5.)*rzmzm
f=(delta**2/5.)*rpzpz
cp=cpzfs
d=1.0/taur
e=wi
ea=omega3
eb=omega2
ec=omega1
else
a=(delta**2/5.)*rpzpz
b=(delta**2/5.)*rzmzm
f=(delta**2/5.)*rpmpm
cp=cpzee
d=1.0/taur
e=wi
ea=omega1
eb=omega2
ec=omega3
endif
c*************************** writing out eigenvalues for R2e
*****************
c rt2a=(a+b+sqrt((a-b)**2+4*cp**2))/2
c rt2b=(a+b-sqrt((a-b)**2+4*cp**2))/2
c rt2c=f
c******************************************************************************
a44=(b**2*d-b*(cp**2-2*d**2)+d*(-(cp**2)+d**2+(e+eb)**2)+
& a*((b+d)**2+(e+eb)**2))/(cp**4+2*cp**2*(-d*(b+d)+
& (e+ea)*(e+eb))+a**2*((b+d)**2+(e+eb)**2)+(d**2+(e+ea)**2)*
& ((b+d)**2+(e+eb)**2)+2*a*(b**2*d-b*(cp**2-2*d**2)+
& d*(-(cp**2)+d**2+(e+eb)**2)))
a55=(-(cp**2)*d+a**2*(b+d)+a*(-(cp**2)+2*d*(b+d))+(b+d)*
& (d**2+(e+ea)**2))/(cp**4+2*cp**2*(-d*(b+d)+
& (e+ea)*(e+eb))+a**2*((b+d)**2+(e+eb)**2)+(d**2+(e+ea)**2)*
& ((b+d)**2+(e+eb)**2)+2*a*(b**2*d-b*(cp**2-2*d**2)+
& d*(-(cp**2)+d**2+(e+eb)**2)))
a45=-(cp*(-(cp**2)+(a+d)*(b+d)-(e+ea)*(e+eb)))/
& ((cp**2-(a+d)*(b+d)+(e+ea)*(e+eb))**2+
& (b*(e+ea)+a*(e+eb)+d*(2*e+ea+eb))**2)
a54=a45
taus3=1./(d+f)
a66=taus3/(1+(wi+ec)**2*taus3**2)
a77=a44
a88=a55
a78=a45
a87=a78
a99=a66
a65=0.0d0
a56=0.0d0
a89=0.0d0
a98=0.0d0
if (dpara.gt.ws) then
astore1=a44
astore2=a66
a44=a99
a66=a77
a99=astore1
a77=astore2
a65=a45
a56=a45
a89=a45
a98=a65
a45=0.0d0
a54=0.0d0
a78=0.0d0
a87=0.0d0
endif
c
C******************************************************************************
C******** DEFINITION OF COEFF. MATRIX FOR DIAGONALIZED H0
*********************
C******** IN RHOMBIC ROUTINE
*********************
C******************************************************************************
zrn(1,1)=zr(1,3) !Zeeman
zrn(1,2)=zr(2,3)
zrn(1,3)=zr(3,3)
zrn(2,1)=zr(1,2)
zrn(2,2)=zr(2,2)
zrn(2,3)=zr(3,2)
zrn(3,1)=zr(1,1)
zrn(3,2)=zr(2,1)
zrn(3,3)=zr(3,1)
C******************************************************************************
C************** DIAGONALIZATION OF COEFF. MATRIX
******************************
C************** IN RHOMBIC ROUTINE
******************************
C******************************************************************************
lda=3
num=3
call zgefa(zrn,lda,num,ipvt,info)
job=01
call zgedi(zrn,lda,num,ipvt,det,work,job)
ak11=zrn(1,1)
ak12=zrn(1,2)
ak13=zrn(1,3)
ak21=zrn(2,1)
ak22=zrn(2,2)
ak23=zrn(2,3)
ak31=zrn(3,1)
ak32=zrn(3,2)
ak33=zrn(3,3)
C******************************************************************************
C************* LIOUVILLE SPACE REPRESENTATION OF S-OPERATORS
******************
C************* PROJECTION VECTORS IN RHOMBIC ROUTINE
******************
C******************************************************************************
cz1=ak11*conjg(ak11)-ak31*conjg(ak31)
cz2=ak12*conjg(ak12)-ak32*conjg(ak32)
cz3=ak13*conjg(ak13)-ak33*conjg(ak33)
cz4=ak11*conjg(ak12)-ak31*conjg(ak32)
cz5=ak12*conjg(ak13)-ak32*conjg(ak33)
cz6=ak11*conjg(ak13)-ak31*conjg(ak33)
cz7=ak12*conjg(ak11)-ak32*conjg(ak31)
cz8=ak13*conjg(ak12)-ak33*conjg(ak32)
cz9=ak13*conjg(ak11)-ak33*conjg(ak31)
cp1=-sqrt(2.)*(ak21*conjg(ak31)+ak11*conjg(ak21))
cp2=-sqrt(2.)*(ak22*conjg(ak32)+ak12*conjg(ak22))
cp3=-sqrt(2.)*(ak23*conjg(ak33)+ak13*conjg(ak23))
cp4=-sqrt(2.)*(ak21*conjg(ak32)+ak11*conjg(ak22))
cp5=-sqrt(2.)*(ak22*conjg(ak33)+ak12*conjg(ak23))
cp6=-sqrt(2.)*(ak21*conjg(ak33)+ak11*conjg(ak23))
cp7=-sqrt(2.)*(ak22*conjg(ak31)+ak12*conjg(ak21))
cp8=-sqrt(2.)*(ak23*conjg(ak32)+ak13*conjg(ak22))
cp9=-sqrt(2.)*(ak23*conjg(ak31)+ak13*conjg(ak21))
cm1=sqrt(2.)*(ak31*conjg(ak21)+ak21*conjg(ak11))
cm2=sqrt(2.)*(ak32*conjg(ak22)+ak22*conjg(ak12))
cm3=sqrt(2.)*(ak33*conjg(ak23)+ak23*conjg(ak13))
cm4=sqrt(2.)*(ak31*conjg(ak22)+ak21*conjg(ak12))
cm5=sqrt(2.)*(ak32*conjg(ak23)+ak22*conjg(ak13))
cm6=sqrt(2.)*(ak31*conjg(ak23)+ak21*conjg(ak13))
cm7=sqrt(2.)*(ak32*conjg(ak21)+ak22*conjg(ak11))
cm8=sqrt(2.)*(ak33*conjg(ak22)+ak23*conjg(ak12))
cm9=sqrt(2.)*(ak33*conjg(ak21)+ak23*conjg(ak11))
C******************************************************************************
C**************** PROJECTION OF SUPERMATRIX
***********************************
C**************** IN RHOMBIC ROUTINE
***********************************
C******************************************************************************
ammcz1=a11*1e9*cz1+a12*1e9*cz2+a13*1e9*cz3
ammcz2=a21*1e9*cz1+a22*1e9*cz2+a23*1e9*cz3
ammcz3=a31*1e9*cz1+a32*1e9*cz2+a33*1e9*cz3
ammcz4=a44*1e9*cz4+a45*1e9*cz5
ammcz5=a55*1e9*cz5+a54*1e9*cz4+a56*1e9*cz6
ammcz6=a66*1e9*cz6+a65*1e9*cz5
ammcz7=a77*1e9*cz7+a78*1e9*cz8
ammcz8=a88*1e9*cz8+a87*1e9*cz7+a89*1e9*cz9
ammcz9=a99*1e9*cz9+a98*1e9*cz8
ammcp1=a11*1e9*cp1+a12*1e9*cp2+a13*1e9*cp3
ammcp2=a21*1e9*cp1+a22*1e9*cp2+a23*1e9*cp3
ammcp3=a31*1e9*cp1+a32*1e9*cp2+a33*1e9*cp3
ammcp4=a44*1e9*cp4+a45*1e9*cp5
ammcp5=a55*1e9*cp5+a54*1e9*cp4+a56*1e9*cp6
ammcp6=a66*1e9*cp6+a65*1e9*cp5
ammcp7=a77*1e9*cp7+a78*1e9*cp8
ammcp8=a88*1e9*cp8+a87*1e9*cp7+a89*1e9*cp9
ammcp9=a99*1e9*cp9+a98*1e9*cp8
ammcm1=a11*1e9*cm1+a12*1e9*cm2+a13*1e9*cm3
ammcm2=a21*1e9*cm1+a22*1e9*cm2+a23*1e9*cm3
ammcm3=a31*1e9*cm1+a32*1e9*cm2+a33*1e9*cm3
ammcm4=a44*1e9*cm4+a45*1e9*cm5
ammcm5=a55*1e9*cm5+a54*1e9*cm4+a56*1e9*cm6
ammcm6=a66*1e9*cm6+a65*1e9*cm5
ammcm7=a77*1e9*cm7+a78*1e9*cm8
ammcm8=a88*1e9*cm8+a87*1e9*cm7+a89*1e9*cm9
ammcm9=a99*1e9*cm9+a98*1e9*cm8
C******************************************************************************
C*************** SPECTRAL DENSITIES OF S
**************************************
C*************** IN RHOMBIC ROUTINE
**************************************
C******************************************************************************
!cz M cz
primo=(conjg(cz1)*ammcz1+conjg(cz2)*ammcz2+conjg(cz3)*ammcz3)+
& (conjg(cz4)*ammcz4+conjg(cz5)*ammcz5+conjg(cz6)*ammcz6+
& conjg(cz7)*ammcz7+conjg(cz8)*ammcz8+conjg(cz9)*ammcz9)
!cp M cp
secondo=(conjg(cp1)*ammcp1+conjg(cp2)*ammcp2+conjg(cp3)
& *ammcp3)+
& (conjg(cp4)*ammcp4+conjg(cp5)*ammcp5+conjg(cp6)*ammcp6+
& conjg(cp7)*ammcp7+conjg(cp8)*ammcp8+conjg(cp9)*ammcp9)
!cp M cm
terzo=(conjg(cp1)*ammcm1+conjg(cp2)*ammcm2+conjg(cp3)*ammcm3)+
& (conjg(cp4)*ammcm4+conjg(cp5)*ammcm5+conjg(cp6)*ammcm6+
& conjg(cp7)*ammcm7+conjg(cp8)*ammcm8+conjg(cp9)*ammcm9)
!cp M cz
quarto=(conjg(cp1)*ammcz1+conjg(cp2)*ammcz2+conjg(cp3)*
& ammcz3)+
& (conjg(cp4)*ammcz4+conjg(cp5)*ammcz5+conjg(cp6)*ammcz6+
& conjg(cp7)*ammcz7+conjg(cp8)*ammcz8+conjg(cp9)*ammcz9)
!cm M cp
quinto=(conjg(cm1)*ammcp1+conjg(cm2)*ammcp2+conjg(cm3)*
& ammcp3)+
& (conjg(cm4)*ammcp4+conjg(cm5)*ammcp5+conjg(cm6)*ammcp6+
& conjg(cm7)*ammcp7+conjg(cm8)*ammcp8+conjg(cm9)*ammcp9)
!cm M cm
sesto=(conjg(cm1)*ammcm1+conjg(cm2)*ammcm2+conjg(cm3)*ammcm3)+
& (conjg(cm4)*ammcm4+conjg(cm5)*ammcm5+conjg(cm6)*ammcm6+
& conjg(cm7)*ammcm7+conjg(cm8)*ammcm8+conjg(cm9)*ammcm9)
!cm M cz
settimo=(conjg(cm1)*ammcz1+conjg(cm2)*ammcz2+conjg(cm3)*
& ammcz3)+
& (conjg(cm4)*ammcz4+conjg(cm5)*ammcz5+conjg(cm6)*ammcz6+
& conjg(cm7)*ammcz7+conjg(cm8)*ammcz8+conjg(cm9)*ammcz9)
!cz M cp
ottavo=(conjg(cz1)*ammcp1+conjg(cz2)*ammcp2+conjg(cz3)*
& ammcp3)+
& (conjg(cz4)*ammcp4+conjg(cz5)*ammcp5+conjg(cz6)*ammcp6+
& conjg(cz7)*ammcp7+conjg(cz8)*ammcp8+conjg(cz9)*ammcp9)
!cz M cm
onono=(conjg(cz1)*ammcm1+conjg(cz2)*ammcm2+conjg(cz3)
& *ammcm3)+
& (conjg(cz4)*ammcm4+conjg(cz5)*ammcm5+conjg(cz6)*ammcm6+
& conjg(cz7)*ammcm7+conjg(cz8)*ammcm8+conjg(cz9)*ammcm9)
C******************************************************************************
C**************** SPECTRAL DENSITIES OF G
*************************************
C**************** IN RHOMBIC ROUTINE
*************************************
C******************************************************************************
C AUTOCORRELATION FUNCTIONS
C(1,1,11)=secondo
C(1,1,12)=sesto
C(1,1,13)=quinto
C(1,1,14)= terzo
C(1,1,15)=quarto
!
C(1,1,16)= ottavo
C(1,1,17)=settimo
C(1,1,18)=onono
C(1,1,19)=primo
567 CONTINUE
disp=1
RETURN
END
SUBROUTINE POWELL(P,XI,N,NP,FTOL,ITER,FRET,NMX)
C PERFORMES THE FITTING PROCEDURE
IMPLICIT REAL*8(A-H,O-Z)
PARAMETER (NMAX=20,ITMAX=2000)
DIMENSION P(NP),XI(NP,NP),PT(NMAX),PTT(NMAX),XIT(NMAX)
COMMON /NMXCOM/ NMX1
COMMON /NVNP/ NV
NMX1=NMX
NV=NP
CALL XISTEP(NV,XI,P)
CALL FUNCZFS(P,FUNC,NMX,NP)
FRET=FUNC
DO 11 J=1,N
PT(J)=P(J)
11 CONTINUE
ITER=0
1 ITER=ITER+1
FP=FRET
IBIG=0
DEL=0.
DO 13 I=1,N
DO 12 J=1,N
XIT(J)=XI(J,I)
12 CONTINUE
FPTT=FRET
CALL LINMIN(P,XIT,N,FRET)
IF(ABS(FPTT-FRET).GT.DEL)THEN
DEL=ABS(FPTT-FRET)
IBIG=I
ENDIF
13 CONTINUE
IF(2.*ABS(FP-FRET).LE.FTOL*(ABS(FP)+ABS(FRET)))RETURN
IF(ITER.EQ.ITMAX) PAUSE 'POWELL EXCEEDING MAXIMUM ITERATIONS.'
DO 14 J=1,N
PTT(J)=2.*P(J)-PT(J)
XIT(J)=P(J)-PT(J)
PT(J)=P(J)
14 CONTINUE
CALL FUNCZFS(PTT,FUNC,NMX,NP)
FPTT=FUNC
IF(FPTT.GE.FP)GO TO 1
T=2.*(FP-2.*FRET+FPTT)*(FP-FRET-DEL)**2-DEL*(FP-FPTT)**2
IF(T.GE.0.)GO TO 1
CALL LINMIN(P,XIT,N,FRET)
DO 15 J=1,N
XI(J,IBIG)=XIT(J)
15 CONTINUE
GO TO 1
END
SUBROUTINE POWELLINT(P,XI,N,NP,FTOL,ITER,FRET,NMX)
C PERFORMES THE INTERNAL FITTING PROCEDURE
IMPLICIT REAL*8(A-H,O-Z)
PARAMETER (NMAX=20,ITMAX=2000)
DIMENSION P(NP),XI(NP,NP),PT(NMAX),PTT(NMAX),XIT(NMAX)
COMMON /NMXCOM/ NMX1
COMMON /NVNP2/ NV
NMX1=NMX
NV=NP
CALL XISTEP2(NV,XI,P)
CALL FUNCINT(P,FUNC,NMX,NP)
FRET=FUNC
DO 11 J=1,N
PT(J)=P(J)
11 CONTINUE
ITER=0
1 ITER=ITER+1
FP=FRET
IBIG=0
DEL=0.
DO 13 I=1,N
DO 12 J=1,N
XIT(J)=XI(J,I)
12 CONTINUE
! FPTT=FRET
CALL LINMIN2(P,XIT,N,FRET)
! IF(ABS(FPTT-FRET).GT.DEL)THEN
IF(ABS(FP-FRET).GT.DEL)THEN
! DEL=ABS(FPTT-FRET)
DEL=ABS(FP-FRET)
IBIG=I
ENDIF
13 CONTINUE
IF(2.*ABS(FP-FRET).LE.FTOL*(ABS(FP)+ABS(FRET)))RETURN
IF(ITER.EQ.ITMAX) RETURN !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
IF(ITER.EQ.ITMAX) PAUSE 'POWELL EXCEEDING MAXIMUM ITERATIONS.'
DO 14 J=1,N
PTT(J)=2.*P(J)-PT(J)
XIT(J)=P(J)-PT(J)
PT(J)=P(J)
14 CONTINUE
CALL FUNCINT(PTT,FUNC,NMX,NP)
FPTT=FUNC
IF(FPTT.GE.FP)GO TO 1
T=2.*(FP-2.*FRET+FPTT)*(FP-FRET-DEL)**2-DEL*(FP-FPTT)**2
IF(T.GE.0.)GO TO 1
CALL LINMIN2(P,XIT,N,FRET)
DO 15 J=1,N
XI(J,IBIG)=XIT(J)
15 CONTINUE
GO TO 1
END
SUBROUTINE LINMIN(P,XI,N,FRET)
IMPLICIT REAL*8(A-H,O-Z)
PARAMETER (NMAX=50,TOL=1.E-4)
EXTERNAL F1DIM
DIMENSION P(N),XI(N)
COMMON /F1COM/ PCOM(NMAX),XICOM(NMAX),NCOM
NCOM=N
DO 11 J=1,N
PCOM(J)=P(J)
XICOM(J)=XI(J)
11 CONTINUE
AX=0.
XX=1
! BX=2
CALL MNBRAK(AX,XX,BX,FA,FX,FB,F1DIM)
FRET=MRENT(AX,XX,BX,F1DIM,TOL,XMIN)
DO 12 J=1,N
XI(J)=XMIN*XI(J)
P(J)=P(J)+XI(J)
12 CONTINUE
RETURN
END
SUBROUTINE LINMIN2(P,XI,N,FRET)
IMPLICIT REAL*8(A-H,O-Z)
PARAMETER (NMAX=50,TOL=1.E-4)
EXTERNAL F2DIM
DIMENSION P(N),XI(N)
COMMON /F2COM/ PCOM(NMAX),XICOM(NMAX),NCOM
NCOM=N
DO 11 J=1,N
PCOM(J)=P(J)
XICOM(J)=XI(J)
11 CONTINUE
AX=0.
XX=1
! BX=2
CALL MNBRAK2(AX,XX,BX,FA,FX,FB,F2DIM)
FRET=MRENT2(AX,XX,BX,F2DIM,TOL,XMIN)
DO 12 J=1,N
XI(J)=XMIN*XI(J)
P(J)=P(J)+XI(J)
12 CONTINUE
RETURN
END
FUNCTION F1DIM(X)
IMPLICIT REAL*8(A-H,O-Z)
PARAMETER (NMAX=50)
COMMON /F1COM/ PCOM(NMAX),XICOM(NMAX),NCOM
COMMON /NMXCOM/ NMX1
COMMON /NVNP/ NV
DIMENSION XT(NMAX)
DO 11 J=1,NCOM
XT(J)=PCOM(J)+X*XICOM(J)
11 CONTINUE
CALL FUNCZFS(XT,FUNC,NMX1,NV)
F1DIM=FUNC
RETURN
END
FUNCTION F2DIM(X)
IMPLICIT REAL*8(A-H,O-Z)
PARAMETER (NMAX=50)
COMMON /F2COM/ PCOM(NMAX),XICOM(NMAX),NCOM
COMMON /NMXCOM/ NMX1
COMMON /NVNP2/ NV
DIMENSION XT(NMAX)
DO 11 J=1,NCOM
XT(J)=PCOM(J)+X*XICOM(J)
11 CONTINUE
CALL FUNCINT(XT,FUNC,NMX1,NV)
F2DIM=FUNC
RETURN
END
SUBROUTINE MNBRAK(AX,BX,CX,FA,FB,FC,F1DIM)
IMPLICIT REAL*8(A-H,O-Z)
PARAMETER (GOLD=1.618034, GLIMIT=100., TINY=1.E-20)
FA=F1DIM(AX)
FB=F1DIM(BX)
IF(FB.GT.FA)THEN
DUM=AX
AX=BX
BX=DUM
DUM=FB
FB=FA
FA=DUM
ENDIF
CX=BX+GOLD*(BX-AX)
FC=F1DIM(CX)
1 IF(FB.GE.FC)THEN
R=(BX-AX)*(FB-FC)
Q=(BX-CX)*(FB-FA)
U=BX-((BX-CX)*Q-(BX-AX)*R)/(2.*SIGN(MAX(ABS(Q-R),TINY),Q-R))
ULIM=BX+GLIMIT*(CX-BX)
IF((BX-U)*(U-CX).GT.0.)THEN
FU=F1DIM(U)
IF(FU.LT.FC)THEN
AX=BX
FA=FB
BX=U
FB=FU
GO TO 1
ELSE IF(FU.GT.FB)THEN
CX=U
FC=FU
GO TO 1
ENDIF
U=CX+GOLD*(CX-BX)
FU=F1DIM(U)
ELSE IF((CX-U)*(U-ULIM).GT.0.)THEN
FU=F1DIM(U)
IF(FU.LT.FC)THEN
BX=CX
CX=U
U=CX+GOLD*(CX-BX)
FB=FC
FC=FU
FU=F1DIM(U)
ENDIF
ELSE IF((U-ULIM)*(ULIM-CX).GE.0.)THEN
U=ULIM
FU=F1DIM(U)
ELSE
U=CX+GOLD*(CX-BX)
FU=F1DIM(U)
ENDIF
AX=BX
BX=CX
CX=U
FA=FB
FB=FC
FC=FU
GO TO 1
ENDIF
RETURN
END
SUBROUTINE MNBRAK2(AX,BX,CX,FA,FB,FC,F2DIM)
IMPLICIT REAL*8(A-H,O-Z)
PARAMETER (GOLD=1.618034, GLIMIT=100., TINY=1.E-20)
FA=F2DIM(AX)
FB=F2DIM(BX)
IF(FB.GT.FA)THEN
DUM=AX
AX=BX
BX=DUM
DUM=FB
FB=FA
FA=DUM
ENDIF
CX=BX+GOLD*(BX-AX)
FC=F2DIM(CX)
1 IF(FB.GE.FC)THEN
R=(BX-AX)*(FB-FC)
Q=(BX-CX)*(FB-FA)
U=BX-((BX-CX)*Q-(BX-AX)*R)/(2.*SIGN(MAX(ABS(Q-R),TINY),Q-R))
ULIM=BX+GLIMIT*(CX-BX)
IF((BX-U)*(U-CX).GT.0.)THEN
FU=F2DIM(U)
IF(FU.LT.FC)THEN
AX=BX
FA=FB
BX=U
FB=FU
GO TO 1
ELSE IF(FU.GT.FB)THEN
CX=U
FC=FU
GO TO 1
ENDIF
U=CX+GOLD*(CX-BX)
FU=F2DIM(U)
ELSE IF((CX-U)*(U-ULIM).GT.0.)THEN
FU=F2DIM(U)
IF(FU.LT.FC)THEN
BX=CX
CX=U
U=CX+GOLD*(CX-BX)
FB=FC
FC=FU
FU=F2DIM(U)
ENDIF
ELSE IF((U-ULIM)*(ULIM-CX).GE.0.)THEN
U=ULIM
FU=F2DIM(U)
ELSE
U=CX+GOLD*(CX-BX)
FU=F2DIM(U)
ENDIF
AX=BX
BX=CX
CX=U
FA=FB
FB=FC
FC=FU
GO TO 1
ENDIF
RETURN
END
FUNCTION MRENT(AX,BX,CX,F1DIM,TOL,XMIN)
IMPLICIT REAL*8(A-H,O-Z)
PARAMETER (ITMAX=100,CGOLD=.3819660,ZEPS=1.0E-10)
A=MIN(AX,CX)
B=MAX(AX,CX)
V=BX
W=V
X=V
E=0.
FX=F1DIM(X)
FV=FX
FW=FX
DO 11 ITER=1,ITMAX
XM=0.5*(A+B)
TOL1=TOL*ABS(X)+ZEPS
TOL2=2.*TOL1
IF(ABS(X-XM).LE.(TOL2-.5*(B-A))) GOTO 3
IF(ABS(E).GT.TOL1) THEN
R=(X-W)*(FX-FV)
Q=(X-V)*(FX-FW)
P=(X-V)*Q-(X-W)*R
Q=2.*(Q-R)
IF(Q.GT.0.) P=-P
Q=ABS(Q)
ETEMP=E
E=D
IF(ABS(P).GE.ABS(.5*Q*ETEMP).OR.P.LE.Q*(A-X).OR.
* P.GE.Q*(B-X)) GOTO 1
D=P/Q
U=X+D
IF(U-A.LT.TOL2 .OR. B-U.LT.TOL2) D=SIGN(TOL1,XM-X)
GOTO 2
ENDIF
1 IF(X.GE.XM) THEN
E=A-X
ELSE
E=B-X
ENDIF
D=CGOLD*E
2 IF(ABS(D).GE.TOL1) THEN
U=X+D
ELSE
U=X+SIGN(TOL1,D)
ENDIF
FU=F1DIM(U)
IF(FU.LE.FX) THEN
IF(U.GE.X) THEN
A=X
ELSE
B=X
ENDIF
V=W
FV=FW
W=X
FW=FX
X=U
FX=FU
ELSE
IF(U.LT.X) THEN
A=U
ELSE
B=U
ENDIF
IF(FU.LE.FW .OR. W.EQ.X) THEN
V=W
FV=FW
W=U
FW=FU
ELSE IF(FU.LE.FV .OR. V.EQ.X .OR. V.EQ.W) THEN
V=U
FV=FU
ENDIF
ENDIF
11 CONTINUE
! PAUSE 'MRENT EXCEED MAXIMUM ITERATIONS.'
WRITE(6,*) 'MRENT EXCEED MAXIMUM ITERATIONS.'
RETURN !!!!!!!!!!!!!!!!!!!!!
3 XMIN=X
MRENT=FX
RETURN
END
FUNCTION MRENT2(AX,BX,CX,F2DIM,TOL,XMIN)
IMPLICIT REAL*8(A-H,O-Z)
PARAMETER (ITMAX=100,CGOLD=.3819660,ZEPS=1.0E-10)
A=MIN(AX,CX)
B=MAX(AX,CX)
V=BX
W=V
X=V
E=0.
FX=F2DIM(X)
FV=FX
FW=FX
DO 11 ITER=1,ITMAX
XM=0.5*(A+B)
TOL1=TOL*ABS(X)+ZEPS
TOL2=2.*TOL1
IF(ABS(X-XM).LE.(TOL2-.5*(B-A))) GOTO 3
IF(ABS(E).GT.TOL1) THEN
R=(X-W)*(FX-FV)
Q=(X-V)*(FX-FW)
P=(X-V)*Q-(X-W)*R
Q=2.*(Q-R)
IF(Q.GT.0.) P=-P
Q=ABS(Q)
ETEMP=E
E=D
IF(ABS(P).GE.ABS(.5*Q*ETEMP).OR.P.LE.Q*(A-X).OR.
* P.GE.Q*(B-X)) GOTO 1
D=P/Q
U=X+D
IF(U-A.LT.TOL2 .OR. B-U.LT.TOL2) D=SIGN(TOL1,XM-X)
GOTO 2
ENDIF
1 IF(X.GE.XM) THEN
E=A-X
ELSE
E=B-X
ENDIF
D=CGOLD*E
2 IF(ABS(D).GE.TOL1) THEN
U=X+D
ELSE
U=X+SIGN(TOL1,D)
ENDIF
FU=F2DIM(U)
IF(FU.LE.FX) THEN
IF(U.GE.X) THEN
A=X
ELSE
B=X
ENDIF
V=W
FV=FW
W=X
FW=FX
X=U
FX=FU
ELSE
IF(U.LT.X) THEN
A=U
ELSE
B=U
ENDIF
IF(FU.LE.FW .OR. W.EQ.X) THEN
V=W
FV=FW
W=U
FW=FU
ELSE IF(FU.LE.FV .OR. V.EQ.X .OR. V.EQ.W) THEN
V=U
FV=FU
ENDIF
ENDIF
11 CONTINUE
! PAUSE 'MRENT EXCEED MAXIMUM ITERATIONS.'
WRITE(6,*) 'MRENT2 EXCEED MAXIMUM ITERATIONS.'
RETURN !!!!!!!!!!!!!!!!!!!!!
3 XMIN=X
MRENT2=FX
RETURN
END
SUBROUTINE XISTEP(NV,XI,P)
DIMENSION XI(NV,NV),P(NV)
IMPLICIT REAL*8(A-H,O-Z)
COMMON /ALFASTEP/ ALFASTEP
DO 13 I=1,NV
DO 13 J=1,NV
13 XI(I,J)=ALFASTEP*P(I)
DO 14 I=1,NV
14 XI(I,I)=-ALFASTEP*P(I)
RETURN
END
SUBROUTINE XISTEP2(NV,XI,P)
DIMENSION XI(NV,NV),P(NV)
IMPLICIT REAL*8(A-H,O-Z)
COMMON /ALFASTEP/ ALFASTEP
DO 13 I=1,NV
DO 13 J=1,NV
13 XI(I,J)=ALFASTEP*P(I)
DO 14 I=1,NV
14 XI(I,I)=-ALFASTEP*P(I)
RETURN
END
SUBROUTINE F01BCF(N,TOL,Z,IZ,W,IW,D,E,C,S)
C MARK 3 RELEASE. HTTP://MEAMI.ORG/ 2009.
C MARK 4 REVISED.
C MARK 4.5 REVISED
C MARK 5C REVISED
C MARK 6B REVISED IER-125 IER-127 (OCT 2009)
C MARK 11 REVISED. VECTORISATION (OCT 2009).
C MARK 11.5(F77) REVISED. (OCTO 2009.)
C
C
C TRECX2
C F01BCF REDUCES A COMPLEX HERMITIAN MATRIX TO REAL
C TRIDIAGONAL FORM FROM WHICH THE EIGENVALUES AND EIGENVECTORS
C CAN BE FOUND USING SUBROUTINE F02AYF,(CXTQL2). THE HERMITIAN
C MATRIX A=A(1) IS REDUCED TO THE TRIDIAGONAL MATRIX A(N-1) BY
C N-2 UNITARY TRANSFORMATIONS. THE HOUSEHOLDER REDUCTION ITSELF
C DOES NOT GIVE A REAL TRIDIAGONAL MATRIX, THE OFF-DIAGONAL
C ELEMENTS ARE COMPLEX. THEY ARE SUBSEQUENTLY MADE REAL BY A
C DIAGONAL TRANSFORMATION.
C OCTOBER 21ST. 2009
C
C .. SCALAR ARGUMENTS ..
DOUBLE PRECISION TOL
INTEGER IW, IZ, N
C .. ARRAY ARGUMENTS ..
DOUBLE PRECISION C(N), D(N), E(N), S(N), W(IW,N), Z(IZ,N)
C .. LOCAL SCALARS ..
DOUBLE PRECISION CO, F, FI, FR, G, GI, GR, H, HH, R, SI
INTEGER I, II, J, K, L
C .. EXTERNAL SUBROUTINES ..
EXTERNAL F01BCY, F01BCZ
C .. INTRINSIC FUNCTIONS ..
INTRINSIC ABS, SQRT
C .. EXECUTABLE STATEMENTS ..
DO 20 I = 1, N
D(I) = Z(N,I)
E(I) = -W(N,I)
20 CONTINUE
IF (N.EQ.1) GO TO 540
DO 360 II = 2, N
I = N - II + 2
L = I - 2
G = 0.0D0
FR = D(I-1)
FI = E(I-1)
IF (L.EQ.0) GO TO 60
DO 40 K = 1, L
G = G + D(K)*D(K) + E(K)*E(K)
40 CONTINUE
60 H = G + FR*FR + FI*FI
C L IS NOW I-1
L = L + 1
IF (ABS(FR)+ABS(FI).NE.0.0D0) GO TO 80
R = 0.0D0
CO = 1.0D0
C(I) = 1.0D0
SI = 0.0D0
S(I) = 0.0D0
GO TO 140
80 IF (ABS(FR).LT.ABS(FI)) GO TO 100
R = ABS(FR)*SQRT(1.0D0+(FI/FR)**2)
GO TO 120
100 R = ABS(FI)*SQRT(1.0D0+(FR/FI)**2)
120 SI = FI/R
S(I) = -SI
CO = FR/R
C(I) = CO
140 IF (G.LE.TOL) GO TO 280
G = -SQRT(H)
E(I) = G
C E(I) HAS A REAL VALUE
H = H - R*G
C S*S + SR
D(I-1) = (R-G)*CO
E(I-1) = (R-G)*SI
DO 160 J = 1, L
Z(J,I) = D(J)
W(J,I) = E(J)
160 CONTINUE
CALL F01BCZ(Z,IZ,W,IW,L,D,E,C,S)
C FORM P
DO 180 J = 1, L
C(J) = C(J)/H
S(J) = S(J)/H
180 CONTINUE
FR = 0.0D0
DO 200 J = 1, L
FR = FR + C(J)*D(J) + S(J)*E(J)
200 CONTINUE
C FORM K
HH = FR/(H+H)
C FORM Q
DO 220 J = 1, L
C(J) = C(J) - HH*D(J)
S(J) = S(J) - HH*E(J)
220 CONTINUE
C NOW FORM REDUCED A
DO 260 J = 1, L
FR = D(J)
FI = E(J)
GR = C(J)
GI = S(J)
DO 240 K = J, L
Z(K,J) = (((Z(K,J)-GR*D(K))-GI*E(K))-FR*C(K)) - FI*S(K)
W(K,J) = (((W(K,J)-GR*E(K))+GI*D(K))-FR*S(K)) + FI*C(K)
240 CONTINUE
D(J) = Z(L,J)
Z(I,J) = 0.0D0
E(J) = -W(L,J)
W(I,J) = 0.0D0
W(J,J) = 0.0D0
260 CONTINUE
GO TO 340
280 E(I) = R
H = 0.0D0
DO 300 J = 1, L
Z(J,I) = D(J)
W(J,I) = E(J)
300 CONTINUE
DO 320 J = 1, L
Z(I,J) = 0.0D0
D(J) = Z(I-1,J)
W(I,J) = 0.0D0
E(J) = -W(I-1,J)
320 CONTINUE
340 D(I) = H
360 CONTINUE
C WE NOW FORM THE PRODUCT OF THE
C HOUSEHOLDER MATRICES, OVERWRITING
C ON Z AND W
DO 500 I = 2, N
L = I - 1
Z(N,L) = Z(L,L)
Z(L,L) = 1.0D0
W(N,L) = E(L)
W(L,L) = 0.0D0
H = D(I)
IF (H.EQ.0.0D0) GO TO 460
DO 380 K = 1, L
D(K) = 0.0D0
E(K) = 0.0D0
380 CONTINUE
CALL F01BCY(Z,IZ,W,IW,L,L,Z(1,I),W(1,I),D,E)
DO 400 K = 1, L
D(K) = D(K)/H
E(K) = -E(K)/H
400 CONTINUE
DO 440 J = 1, L
DO 420 K = 1, L
Z(K,J) = Z(K,J) - Z(K,I)*D(J) + W(K,I)*E(J)
W(K,J) = W(K,J) - Z(K,I)*E(J) - W(K,I)*D(J)
420 CONTINUE
440 CONTINUE
460 DO 480 J = 1, L
Z(J,I) = 0.0D0
W(J,I) = 0.0D0
480 CONTINUE
500 CONTINUE
W(N,N) = E(N)
DO 520 I = 1, N
D(I) = Z(N,I)
Z(N,I) = 0.0D0
E(I) = W(N,I)
W(N,I) = 0.0D0
520 CONTINUE
540 Z(N,N) = 1.0D0
W(N,N) = 0.0D0
E(1) = 0.0D0
C NOW WE MULTIPLY BY THE
C COSTHETA + I SINTHETA COLUMN
C FACTORS
CO = 1.0D0
SI = 0.0D0
IF (N.EQ.1) RETURN
DO 580 I = 2, N
F = CO*C(I) - SI*S(I)
SI = CO*S(I) + SI*C(I)
CO = F
DO 560 J = 1, N
F = Z(J,I)*CO - W(J,I)*SI
W(J,I) = Z(J,I)*SI + W(J,I)*CO
Z(J,I) = F
560 CONTINUE
580 CONTINUE
RETURN
END
SUBROUTINE F01BCY(AR,IAR,AI,IAI,M,N,BR,BI,CR,CI)
C MARK 11 RELEASE. HTTP://MEAMI.ORG/ 2009.
C MARK 11.5(F77) REVISED. (OCTO 2009.)
C
C COMPUTES C = C + (A**H)*B (COMPLEX) WHERE
C A IS RECTANGULAR M BY N.
C C MUST BE DISTINCT FROM B.
C
C
C .. SCALAR ARGUMENTS ..
INTEGER IAI, IAR, M, N
C .. ARRAY ARGUMENTS ..
DOUBLE PRECISION AI(IAI,N), AR(IAR,N), BI(M), BR(M), CI(N), CR
(N)
C .. LOCAL SCALARS ..
DOUBLE PRECISION XI, XR
INTEGER I, J
C .. EXECUTABLE STATEMENTS ..
DO 40 I = 1, N
XR = CR(I)
XI = CI(I)
DO 20 J = 1, M
XR = XR + AR(J,I)*BR(J) + AI(J,I)*BI(J)
XI = XI + AR(J,I)*BI(J) - AI(J,I)*BR(J)
20 CONTINUE
CR(I) = XR
CI(I) = XI
40 CONTINUE
RETURN
END
SUBROUTINE F01BCZ(AR,IAR,AI,IAI,N,BR,BI,CR,CI)
C MARK 11 RELEASE. HTTP://MEAMI.ORG/ 2009.
C MARK 11.5(F77) REVISED. (OCTO 2009.)
C
C COMPUTES C = A*B (COMPLEX) WHERE
C A IS A HERMITIAN N-BY-N MATRIX,
C WHOSE LOWER TRIANGLE IS STORED IN A.
C C MUST BE DISTINCT FROM B.
C
C
C .. SCALAR ARGUMENTS ..
INTEGER IAI, IAR, N
C .. ARRAY ARGUMENTS ..
DOUBLE PRECISION AI(IAI,N), AR(IAR,N), BI(N), BR(N), CI(N), CR
(N)
C .. LOCAL SCALARS ..
DOUBLE PRECISION YI, YR
INTEGER I, IP1, J, NM1
C .. EXECUTABLE STATEMENTS ..
DO 20 I = 1, N
CR(I) = 0.0D0
CI(I) = 0.0D0
20 CONTINUE
IF (N.EQ.1) GO TO 100
NM1 = N - 1
DO 80 I = 1, NM1
DO 40 J = I, N
CR(J) = CR(J) + AR(J,I)*BR(I) - AI(J,I)*BI(I)
CI(J) = CI(J) + AR(J,I)*BI(I) + AI(J,I)*BR(I)
40 CONTINUE
YR = CR(I)
YI = CI(I)
IP1 = I + 1
DO 60 J = IP1, N
YR = YR + AR(J,I)*BR(J) + AI(J,I)*BI(J)
YI = YI + AR(J,I)*BI(J) - AI(J,I)*BR(J)
60 CONTINUE
CR(I) = YR
CI(I) = YI
80 CONTINUE
100 CR(N) = CR(N) + AR(N,N)*BR(N) - AI(N,N)*BI(N)
CI(N) = CI(N) + AR(N,N)*BI(N) + AI(N,N)*BR(N)
RETURN
END
SUBROUTINE F02AXF(AR,IAR,AI,IAI,N,WR,VR,IVR,VI,IVI,WK1,WK2,WK3,
* IFAIL)
C MARK 3 RELEASE. HTTP://MEAMI.ORG/ COPYRIGHT 2009.
C MARK 4.5 REVISED
C MARK 9 REVISED. IER-327 (OCT 2009).
C MARK 11.5(F77) REVISED. (OCTO 2009.)
C MARK 13 REVISED. USE OF MARK 12 X02 FUNCTIONS (OCT 2009).
C
C EIGENVALUES AND EIGENVECTORS OF A COMPLEX HERMITIAN MATRIX
C 21ST OCTOBER 2009
C
C .. PARAMETERS ..
CHARACTER*6 SRNAME
PARAMETER (SRNAME='F02AXF')
C .. SCALAR ARGUMENTS ..
INTEGER IAI, IAR, IFAIL, IVI, IVR, N
C .. ARRAY ARGUMENTS ..
DOUBLE PRECISION AI(IAI,N), AR(IAR,N), VI(IVI,N), VR(IVR,N),
* WK1(N), WK2(N), WK3(N), WR(N)
C .. LOCAL SCALARS ..
DOUBLE PRECISION A, B, MAX, SQ, SUM, XXXX
INTEGER I, ISAVE, J
C .. LOCAL ARRAYS ..
CHARACTER*1 P01REC(1)
C .. EXTERNAL FUNCTIONS ..
DOUBLE PRECISION X02AJF, X02AKF
INTEGER P01ABF
EXTERNAL X02AJF, X02AKF, P01ABF
C .. EXTERNAL SUBROUTINES ..
EXTERNAL F01BCF, F02AYF
C .. INTRINSIC FUNCTIONS ..
INTRINSIC SQRT
C .. EXECUTABLE STATEMENTS ..
ISAVE = IFAIL
DO 40 I = 1, N
IF (AI(I,I).NE.0.0D0) GO TO 140
DO 20 J = 1, I
VR(I,J) = AR(I,J)
VI(I,J) = AI(I,J)
20 CONTINUE
40 CONTINUE
CALL F01BCF(N,X02AKF()/X02AJF(),VR,IVR,VI,IVI,WR,WK1,WK2,WK3)
IFAIL = 1
CALL F02AYF(N,X02AJF(),WR,WK1,VR,IVR,VI,IVI,IFAIL)
IF (IFAIL.EQ.0) GO TO 60
IFAIL = P01ABF(ISAVE,1,SRNAME,0,P01REC)
RETURN
C NORMALISE
60 DO 120 I = 1, N
SUM = 0.0D0
MAX = 0.0D0
DO 80 J = 1, N
SQ = VR(J,I)*VR(J,I) + VI(J,I)*VI(J,I)
SUM = SUM + SQ
IF (SQ.LE.MAX) GO TO 80
MAX = SQ
A = VR(J,I)
B = VI(J,I)
80 CONTINUE
IF (SUM.EQ.0.0D0) GO TO 120
SUM = 1.0D0/SQRT(SUM*MAX)
DO 100 J = 1, N
SQ = SUM*(VR(J,I)*A+VI(J,I)*B)
VI(J,I) = SUM*(VI(J,I)*A-VR(J,I)*B)
VR(J,I) = SQ
100 CONTINUE
120 CONTINUE
RETURN
140 IFAIL = P01ABF(ISAVE,2,SRNAME,0,P01REC)
RETURN
END
SUBROUTINE F02AYF(N,EPS,D,E,Z,IZ,W,IW,IFAIL)
C MARK 3 RELEASE. HTTP://MEAMI.ORG/ COPYRIGHT 2009.
C MARK 4 REVISED.
C MARK 4.5 REVISED
C MARK 9 REVISED. IER-326 (OCT 2009).
C MARK 11.5(F77) REVISED. (OCTO 2009.)
C
C CXTQL2
C THIS SUBROUTINE FINDS THE EIGENVALUES AND EIGENVECTORS OF A
C HERMITIAN MATRIX, WHICH HAS BEEN REDUCED TO A REAL
C TRIDIAGONAL MATRIX, T, GIVEN WITH ITS DIAGONAL ELEMENTS IN
C THE ARRAY D(N) AND ITS SUB-DIAGONAL ELEMENTS IN THE LAST N
C - 1 STORES OF THE ARRAY E(N), USING QL TRANSFORMATIONS. THE
C EIGENVALUES ARE OVERWRITTEN ON THE DIAGONAL ELEMENTS IN THE
C ARRAY D IN ASCENDING ORDER. THE REAL AND IMAGINARY PARTS OF
C THE EIGENVECTORS ARE FORMED IN THE ARRAYS Z,W(N,N)
C RESPECTIVELY, OVERWRITING THE ACCUMULATED TRANSFORMATIONS AS
C SUPPLIED BY THE SUBROUTINE F01BCF. THE SUBROUTINE WILL FAIL
C IF ALL EIGENVALUES TAKE MORE THAN 30*N ITERATIONS
C 21ST OCTOBER 2009
C
C .. PARAMETERS ..
CHARACTER*6 SRNAME
PARAMETER (SRNAME='F02AYF')
C .. SCALAR ARGUMENTS ..
DOUBLE PRECISION EPS
INTEGER IFAIL, IW, IZ, N
C .. ARRAY ARGUMENTS ..
DOUBLE PRECISION D(N), E(N), W(IW,N), Z(IZ,N)
C .. LOCAL SCALARS ..
DOUBLE PRECISION B, C, F, G, H, P, R, S
INTEGER I, I1, II, ISAVE, J, K, L, M, M1
C .. LOCAL ARRAYS ..
CHARACTER*1 P01REC(1)
C .. EXTERNAL FUNCTIONS ..
INTEGER P01ABF
EXTERNAL P01ABF
C .. INTRINSIC FUNCTIONS ..
INTRINSIC ABS, SQRT
C .. EXECUTABLE STATEMENTS ..
ISAVE = IFAIL
IF (N.EQ.1) GO TO 40
DO 20 I = 2, N
E(I-1) = E(I)
20 CONTINUE
40 E(N) = 0.0D0
B = 0.0D0
F = 0.0D0
J = 30*N
DO 300 L = 1, N
H = EPS*(ABS(D(L))+ABS(E(L)))
IF (B.LT.H) B = H
C LOOK FOR SMALL SUB-DIAG ELEMENT
DO 60 M = L, N
IF (ABS(E(M)).LE.B) GO TO 80
60 CONTINUE
80 IF (M.EQ.L) GO TO 280
100 IF (J.LE.0) GO TO 400
J = J - 1
C FORM SHIFT
G = D(L)
H = D(L+1) - G
IF (ABS(H).GE.ABS(E(L))) GO TO 120
P = H*0.5D0/E(L)
R = SQRT(P*P+1.0D0)
H = P + R
IF (P.LT.0.0D0) H = P - R
D(L) = E(L)/H
GO TO 140
120 P = 2.0D0*E(L)/H
R = SQRT(P*P+1.0D0)
D(L) = E(L)*P/(1.0D0+R)
140 H = G - D(L)
I1 = L + 1
IF (I1.GT.N) GO TO 180
DO 160 I = I1, N
D(I) = D(I) - H
160 CONTINUE
180 F = F + H
C QL TRANSFORMATION
P = D(M)
C = 1.0D0
S = 0.0D0
M1 = M - 1
DO 260 II = L, M1
I = M1 - II + L
G = C*E(I)
H = C*P
IF (ABS(P).LT.ABS(E(I))) GO TO 200
C = E(I)/P
R = SQRT(C*C+1.0D0)
E(I+1) = S*P*R
S = C/R
C = 1.0D0/R
GO TO 220
200 C = P/E(I)
R = SQRT(C*C+1.0D0)
E(I+1) = S*E(I)*R
S = 1.0D0/R
C = C/R
220 P = C*D(I) - S*G
D(I+1) = H + S*(C*G+S*D(I))
C FORM VECTOR
DO 240 K = 1, N
H = Z(K,I+1)
Z(K,I+1) = S*Z(K,I) + C*H
Z(K,I) = C*Z(K,I) - S*H
H = W(K,I+1)
W(K,I+1) = S*W(K,I) + C*H
W(K,I) = C*W(K,I) - S*H
240 CONTINUE
260 CONTINUE
E(L) = S*P
D(L) = C*P
IF (ABS(E(L)).GT.B) GO TO 100
280 D(L) = D(L) + F
300 CONTINUE
C ORDER EIGEN VALUES ANS EIGENVECTORS
DO 380 I = 1, N
K = I
P = D(I)
I1 = I + 1
IF (I1.GT.N) GO TO 340
DO 320 J = I1, N
IF (D(J).GE.P) GO TO 320
K = J
P = D(J)
320 CONTINUE
340 IF (K.EQ.I) GO TO 380
D(K) = D(I)
D(I) = P
DO 360 J = 1, N
P = Z(J,I)
Z(J,I) = Z(J,K)
Z(J,K) = P
P = W(J,I)
W(J,I) = W(J,K)
W(J,K) = P
360 CONTINUE
380 CONTINUE
IFAIL = 0
RETURN
400 IFAIL = P01ABF(ISAVE,1,SRNAME,0,P01REC)
RETURN
END
INTEGER FUNCTION P01ABF(IFAIL,IERROR,SRNAME,NREC,REC)
C MARK 11.5(F77) RELEASE. HTTP://MEAMI.ORG/ 2009.
C MARK 13 REVISED. IER-621 (OCT 1988).
C MARK 13B REVISED. IER-668 (OCT 2009).
C
C P01ABF IS THE ERROR-HANDLING ROUTINE FOR THE NAG LIBRARY.
C
C P01ABF EITHER RETURNS THE VALUE OF IERROR THROUGH THE ROUTINE
C NAME (SOFT FAILURE), OR TERMINATES EXECUTION OF THE PROGRAM
C (HARD FAILURE). DIAGNOSTIC MESSAGES MAY BE OUTPUT.
C
C IF IERROR = 0 (SUCCESSFUL EXIT FROM THE CALLING ROUTINE),
C THE VALUE 0 IS RETURNED THROUGH THE ROUTINE NAME, AND NO
C MESSAGE IS OUTPUT
C
C IF IERROR IS NON-ZERO (ABNORMAL EXIT FROM THE CALLING ROUTINE),
C THE ACTION TAKEN DEPENDS ON THE VALUE OF IFAIL.
C
C IFAIL = 1: SOFT FAILURE, SILENT EXIT (I.E. NO MESSAGES ARE
C OUTPUT)
C IFAIL = -1: SOFT FAILURE, NOISY EXIT (I.E. MESSAGES ARE OUTPUT)
C IFAIL =-13: SOFT FAILURE, NOISY EXIT BUT STANDARD MESSAGES FROM
C P01ABF ARE SUPPRESSED
C IFAIL = 0: HARD FAILURE, NOISY EXIT
C
C FOR COMPATIBILITY WITH CERTAIN ROUTINES INCLUDED BEFORE MARK 12
C P01ABF ALSO ALLOWS AN ALTERNATIVE SPECIFICATION OF IFAIL IN
WHICH
C IT IS REGARDED AS A DECIMAL INTEGER WITH LEAST SIGNIFICANT
DIGITS
C CBA. THEN
C
C A = 0: HARD FAILURE A = 1: SOFT FAILURE
C B = 0: SILENT EXIT B = 1: NOISY EXIT
C
C EXCEPT HARD FAILURE NOW ALWAYS IMPLIES A NOISY EXIT.
C
C [0]VERSION[1]HTTP://MEAMI.ORG CENTRAL OFFICE.
C
C .. SCALAR ARGUMENTS ..
INTEGER IERROR, IFAIL, NREC
CHARACTER*(*) SRNAME
C .. ARRAY ARGUMENTS ..
CHARACTER*(*) REC(*)
C .. LOCAL SCALARS ..
INTEGER I, NERR
CHARACTER*72 MESS
C .. EXTERNAL SUBROUTINES ..
EXTERNAL P01ABZ, X04AAF, X04BAF
C .. INTRINSIC FUNCTIONS ..
INTRINSIC ABS, MOD
C .. EXECUTABLE STATEMENTS ..
IF (IERROR.NE.0) THEN
C ABNORMAL EXIT FROM CALLING ROUTINE
IF (IFAIL.EQ.-1 .OR. IFAIL.EQ.0 .OR. IFAIL.EQ.-13 .OR.
* (IFAIL.GT.0 .AND. MOD(IFAIL/10,10).NE.0)) THEN
C NOISY EXIT
CALL X04AAF(0,NERR)
DO 20 I = 1, NREC
CALL X04BAF(NERR,REC(I))
20 CONTINUE
IF (IFAIL.NE.-13) THEN
WRITE (MESS,FMT=99999) SRNAME, IERROR
CALL X04BAF(NERR,MESS)
IF (ABS(MOD(IFAIL,10)).NE.1) THEN
C HARD FAILURE
CALL X04BAF(NERR,
* ' ** NAG HARD FAILURE - EXECUTION
TERMINATED'
* )
CALL P01ABZ
ELSE
C SOFT FAILURE
CALL X04BAF(NERR,
* ' ** HTTP://MEAMI.ORG/ 2009. SOFT
FAILURE - CONTROL RETURNED')
END IF
END IF
END IF
END IF
P01ABF = IERROR
RETURN
C
99999 FORMAT (' ** ABNORMAL EXIT FROM HTTP://MEAMI.ORG/ 2009. LIBRARY
ROUTINE ',A,': IFAIL',
* ' =',I6)
END
SUBROUTINE P01ABZ
C MARK 11.5(F77) RELEASE. HTTP://MEAMI.ORG/ 2009.
C
C TERMINATES EXECUTION WHEN A HARD FAILURE OCCURS.
C
C ******************** IMPLEMENTATION NOTE ********************
C THE FOLLOWING STOP STATEMENT MAY BE REPLACED BY A CALL TO AN
C IMPLEMENTATION-DEPENDENT ROUTINE TO DISPLAY A MESSAGE AND/OR
C TO ABORT THE PROGRAM.
C *************************************************************
C .. EXECUTABLE STATEMENTS ..
STOP
END
DOUBLE PRECISION FUNCTION X02AJF()
C MARK 12 RELEASE. HTTP://MEAMI.ORG/ 2009.
C
C RETURNS (1/2)*B**(1-P) IF ROUNDS IS .TRUE.
C RETURNS B**(1-P) OTHERWISE
C
C .. EXECUTABLE STATEMENTS ..
X02AJF = 0.11102230246251568000D-015
RETURN
END
DOUBLE PRECISION FUNCTION X02AKF()
C MARK 12 RELEASE. HTTP://MEAMI.ORG/ 2009.
C
C RETURNS B**(EMIN-1) (THE SMALLEST POSITIVE MODEL NUMBER)
C
C .. EXECUTABLE STATEMENTS ..
X02AKF = 0.22250738585072014000D-307
RETURN
END
SUBROUTINE X04AAF(I,NERR)
C MARK 7 RELEASE. HTTP://MEAMI.ORG/ 2009.
C MARK 7C REVISED IER-190 (OCT 2009)
C MARK 11.5(F77) REVISED. (OCTO 2009.)
C IF I = 0, SETS NERR TO CURRENT ERROR MESSAGE UNIT NUMBER
C (STORED IN NERR1).
C IF I = 1, CHANGES CURRENT ERROR MESSAGE UNIT NUMBER TO
C VALUE SPECIFIED BY NERR.
C
C .. SCALAR ARGUMENTS ..
INTEGER I, NERR
C .. LOCAL SCALARS ..
INTEGER NERR1
C .. SAVE STATEMENT ..
SAVE NERR1
C .. DATA STATEMENTS ..
DATA NERR1/0/
C .. EXECUTABLE STATEMENTS ..
IF (I.EQ.0) NERR = NERR1
IF (I.EQ.1) NERR1 = NERR
RETURN
END
SUBROUTINE X04BAF(NOUT,REC)
C MARK 11.5(F77) RELEASE. HTTP://MEAMI.ORG/ 2009.
C
C X04BAF WRITES THE CONTENTS OF REC TO THE UNIT DEFINED BY NOUT.
C
C TRAILING BLANKS ARE NOT OUTPUT, EXCEPT THAT IF REC IS ENTIRELY
C BLANK, A SINGLE BLANK CHARACTER IS OUTPUT.
C IF NOUT.LT.0, I.E. IF NOUT IS NOT A VALID FORTRAN UNIT
IDENTIFIER,
C THEN NO OUTPUT OCCURS.
C
C .. SCALAR ARGUMENTS ..
INTEGER NOUT
CHARACTER*(*) REC
C .. LOCAL SCALARS ..
INTEGER I
C .. INTRINSIC FUNCTIONS ..
INTRINSIC LEN
C .. EXECUTABLE STATEMENTS ..
IF (NOUT.GE.0) THEN
C REMOVE TRAILING BLANKS
DO 20 I = LEN(REC), 2, -1
IF (REC(I:I).NE.' ') GO TO 40
20 CONTINUE
C WRITE RECORD TO EXTERNAL FILE
40 WRITE (NOUT,FMT=99999) REC(1:I)
END IF
RETURN
C
99999 FORMAT (A)
END
SUBROUTINE GAUINT (BZ,TAUM,NMX)
C PERFORMES THE INTEGRATION ON SPACE
IMPLICIT REAL*8(A-H,O-Z)
COMMON /IPERF/AZ,AY,AX,THETA,RK,TAUC,DPARA,EPARA,PHI,S4
COMMON /TOT/ DPARATOT,EPARATOT,APARTOT,APERTOT,APERTOT2,ACONIND
COMMON /GTENSOR/ GX,GY,GZ
COMMON /STEPGAMMA/ STEPGAMMA
COMMON /RK10/ SPIN, SI
COMMON /TAU1/ TAUS0
COMMON /ECOM/ ECONT,ECROSS
COMMON /T1T2/ IREL
COMMON /TM/ TMUNO,TMUNOCONT,TMUNOCROSS
DIMENSION ARRAY(20),T(20)
DATA ARRAY/
1.,1.,.555556,.888889,.555556,.347855,.652145,.652145,.
1347855,.236927,.478629,.568889,.478629,.236927,.171324,.360762,.46
27914,.467914,.360762,.171324/
DATAT/-.577350,.577350,-.774597,0.,.774597,-.861136,-.
339981,.3399
181,.861136,-.906180,-.538469,0.,.538469,.906180,-.932470,-.
661209,
2-.238619,.238619,.661209,.932470/
TMUNO=0.
TMUNOCONT=0.
TMUNOCROSS=0.
TMUNOTOT=0.
TMUNOTOTCONT=0.
TMUNOTOTCROSS=0.
GAMMA=0.
IF(EPARA.NE.0.OR.GZ.NE.GX.OR.GZ.NE.GY.OR.AX.NE.AY)THEN
STEPGAMMA=20.
ELSE
STEPGAMMA=1.
ENDIF
IF(IREL.NE.1)STEPGAMMA=20.
IF(ACONIND.NE.0)STEPGAMMA=20.
55 A=0
H=1.5708/5.
HP=H/2.
DO 100 I=1,5
G=(2.*A+H)/2.0
A=A+H
SP=0
SPCONT=0
SPCROSS=0
DO 50 J=10,14
BETA=HP*T(J)+G
SP=SP+ARRAY(J)*E(BZ,BETA,THETA,TAUC,NMX,PHI,GAMMA)
SPCONT=SPCONT+ARRAY(J)*ECONT
SPCROSS=SPCROSS+ARRAY(J)*ECROSS
50 CONTINUE
TMUNO=TMUNO+HP*SP
TMUNOCONT=TMUNOCONT+HP*SPCONT
TMUNOCROSS=TMUNOCROSS+HP*SPCROSS
100 CONTINUE
! write(6,*)tmuno !cancellare
TMUNOTOT=TMUNOTOT+TMUNO/STEPGAMMA
TMUNOTOTCONT=TMUNOTOTCONT+TMUNOCONT/STEPGAMMA
TMUNOTOTCROSS=TMUNOTOTCROSS+TMUNOCROSS/STEPGAMMA
TMUNO=0.
TMUNOCONT=0.
TMUNOCROSS=0.
GAMMA=GAMMA+6.28/STEPGAMMA
IF(GAMMA.GE.6.28)THEN
TMUNO=TMUNOTOT
TMUNOCONT=TMUNOTOTCONT
TMUNOCROSS=TMUNOTOTCROSS
GOTO 56
ENDIF
GOTO 55
56 CONTINUE
RETURN
END
SUBROUTINE TUNO(BETA,OMI,THETA,TAUC,NMX,PHI,GAMMA)
C FOUND CONTRIBUTIONS TO T1 TO BE INTEGRATED
IMPLICIT REAL*8(A-H,O-Z)
COMPLEX*16 C(100,100,19)
COMMON /A1/ OM(1000,1000),C
COMMON /RK10/ SPIN, SI
COMMON /TAU1/ TAUS0
COMMON /TAUE/ TAUE
COMMON /MOLFRAZ/ AMOLFRA
COMMON /CONTAT/ ACONT
COMPLEX*16 F0,F1,F2,FMU,FMD
COMPLEX*16 A,B,CI,D,EI,F,G,H,AI
COMPLEX*16 FF1,FF2,FF3,FF4,FF5,FF6,FF7,FF8,FF9
COMPLEX*16 GG1,GG2,GG3,GG4,GG5,GG6,GG7,GG8,GG9
COMPLEX*16 HH1,HH2,HH3,HH4,HH5,HH6,HH7,HH8,HH9
COMPLEX*16 CCF1,CCF2,CCF3,CCF4,CCF5,CCF6,CCF7,CCF8,CCF9
COMPLEX*16 CCG1,CCG2,CCG3,CCG4,CCG5,CCG6,CCG7,CCG8,CCG9
COMPLEX*16 CCH1,CCH2,CCH3,CCH4,CCH5,CCH6,CCH7,CCH8,CCH9
COMPLEX*16 AIA,AJ,ACF,ACG,ACH
COMMON /A3/ T11,T12,T13
CT=COS(BETA)
ST=SIN(BETA)
C CONVERT DEG. ---> RAD (CA)
CONVER = ATAN(1.0)/45.0
CA=COS(THETA* CONVER)
SA=SIN(THETA* CONVER)
CF=COS(PHI* CONVER)
CF2=COS(2.*PHI* CONVER)
SF=SIN(PHI* CONVER)
SF2=SIN(2.*PHI* CONVER)
C **************** RACAH'S NORMALIZED SPHERICAL HARMONICS
F0=-(3.*CA**2-1.)/2.
F1=SQRT(6.)/2.*SA*CA*CMPLX(CF,SF)
FMU=-SQRT(6.)/2.*SA*CA*CMPLX(CF,-SF)
F2=-SQRT(6.)/4.*SA**2*CMPLX(CF2,SF2)
FMD=-SQRT(6.)/4.*SA**2*CMPLX(CF2,-SF2)
C ELEMENTS OF THE ROTATION MATRIX
A=(1.-CT**2)/4.*CMPLX(COS(2.*GAMMA),SIN(2.*GAMMA))
B=(1.+CT)*ST/(2.*SQRT(2.))*CMPLX(COS(GAMMA),SIN(GAMMA))
CI=(1.+CT)**2/4.
D=-(1.-CT)*ST/(2.*SQRT(2.))*CMPLX(COS(GAMMA),SIN(GAMMA))
EI=-ST**2/2.
F=-(1.+CT)*ST/(2.*SQRT(2.))*CMPLX(COS(GAMMA),-SIN(GAMMA))
G=(1.-CT)**2/4.
H=(1.-CT)*ST/(2.*SQRT(2.))*CMPLX(COS(GAMMA),-SIN(GAMMA))
AI=(1.-CT**2)/4.*CMPLX(COS(2.*GAMMA),-SIN(2.*GAMMA))
C ************************
GA=1/TAUC
AKL=0.
c S+S+
HH1=(A*F0*F0/20. + (B+D)*F0*FMU*(1./20.)*SQRT(3.)
& + (CI+G)*F0*FMD*(1./10.)*SQRT(1.5) + EI*FMU*FMU*
(3./20.)
& + (F+H)*FMU*FMD*(3./10.)/SQRT(2.) + AI*FMD*FMD*(3./10.))
& *C(1,1,11)
c S-S-
HH2=(AI*F0*F0/20. + (F+H)*F0*F1*(1./20.)*SQRT(3.)
& + (CI+G)*F0*F2*(1./10.)*SQRT(1.5) + EI*F1*F1*
(3./20.)
& + (B+D)*F1*F2*(3./10.)/SQRT(2.) + A*F2*F2*(3./10.))
& *C(1,1,12)
c S-S+
HH3=(-A*F0*F2*(1./10.)*SQRT(1.5) - B*F2*FMU*(3./10.)
& /SQRT(2.)-CI*F2*FMD*(3./10.) - D*F1*F0*(1./20.)*SQRT(3.)
& - EI*F1*FMU*(3./20.) - F*F1*FMD*(3./10.)/SQRT(2.)
& -G*F0*F0*(1./20.) - H*F0*FMU*(1./20.)*SQRT(3.)
& -AI*F0*FMD*(1./10.)*SQRT(1.5))
& *C(1,1,13)
c S+S-
HH4=(-A*F0*F2*(1./10.)*SQRT(1.5) - D*F2*FMU*(3./10.)
& /SQRT(2.)-G*F2*FMD*(3./10.) - B*F1*F0*(1./20.)*SQRT(3.)
& - EI*F1*FMU*(3./20.) - H*F1*FMD*(3./10.)/SQRT(2.)
& -CI*F0*F0*(1./20.) - F*F0*FMU*(1./20.)*SQRT(3.)
& -AI*F0*FMD*(1./10.)*SQRT(1.5))
& *C(1,1,14)
c S+SZ
HH5=(A*F0*F1*(1./10.)*SQRT(1.5) + B*F0*F0*SQRT(1./50.)
& + CI*F0*FMU*(1./10.)*SQRT(3./2.)+D*F1*FMU*(3./10.)/SQRT(2.)+
& EI*F0*FMU*SQRT(6./100.)+F*FMU*FMU*(3./10.)/SQRT(2.)+
& G*F1*FMD*3./10.+H*F0*FMD*SQRT(6./50.)+AI*FMU*FMD*3./10.)
& *C(1,1,15)
c SZS+
HH6=(A*F0*F1*(1./10.)*SQRT(1.5) + D*F0*F0*SQRT(1./50.)
& + G*F0*FMU*(1./10.)*SQRT(3./2.)+B*F1*FMU*(3./10.)/SQRT(2.)+
& EI*F0*FMU*SQRT(6./100.)+H*FMU*FMU*(3./10.)/SQRT(2.)+
& CI*F1*FMD*3./10.+F*F0*FMD*SQRT(6./50.)+AI*FMU*FMD*3./10.)
& *C(1,1,16)
c S-SZ
HH7=-(AI*F0*FMU*(1./10.)*SQRT(1.5)+H*F0*F0*SQRT(1./50.)
& + G*F0*F1*(1./10.)*SQRT(3./2.)+F*F1*FMU*(3./10.)/SQRT(2.)+
& EI*F0*F1*SQRT(6./100.)+D*F1*F1*(3./10.)/SQRT(2.)+
& CI*FMU*F2*3./10.+B*F0*F2*SQRT(6./50.)+A*F1*F2*3./10.)
& *C(1,1,17)
c SZS-
HH8=-(AI*F0*FMU*(1./10.)*SQRT(1.5)+F*F0*F0*SQRT(1./50.)
& + CI*F0*F1*(1./10.)*SQRT(3./2.)+H*F1*FMU*(3./10.)/SQRT(2.)+
& EI*F0*F1*SQRT(6./100.)+B*F1*F1*(3./10.)/SQRT(2.)+
& G*FMU*F2*3./10.+D*F0*F2*SQRT(6./50.)+A*F1*F2*3./10.)
& *C(1,1,18)
c SZSZ
HH9=((3./10.)*A*F1*F1+B*F1*F0*SQRT(6./50.)+(3./10.)*CI*F1
& *FMU+D*F0*F1*SQRT(6./50)+(2./5.)*EI*F0*F0+F*F0*FMU*SQRT
(6./50.)+
& (3./10.)*G*F1*FMU+H*FMU*F0*SQRT(6./50.)+(3./10.)*AI*FMU*FMU)
& *C(1,1,19)
IF (si.eq.0.and.spin.eq.1.and.acont.eq.0)then
ALL=2.*(REAL(abs(HH1)+abs(HH2)-abs(HH3)-abs(HH4)-abs(HH5)
& +abs(HH6)+abs(HH7)-abs(HH8)-abs(HH9))/FLOAT(NMX))
else
ALL=2.*(REAL(HH1+HH2+HH3+HH4+HH5+HH6+HH7+HH8+HH9)/FLOAT(NMX))
& *GA*1.E+9/(OMI**2+GA**2)+(IMAG(HH1+HH2+HH3+HH4+HH5+HH6+HH7
& +HH8+HH9)/FLOAT(NMX)*OMI-IMAG(HH1+HH2+HH3+HH4+HH5+HH6+HH7+
& HH8+HH9)/FLOAT(NMX)*OMI)
& *1.E+9/(OMI**2+GA**2)
endif
IF (si.eq.0.and.spin.eq.1.and.acont.eq.0)then
T1SC=-ALL
else
T1SC=-ALL-AKL
endif
C DIPOLAR TERM
T11=T1SC
C CONTACT TERM
IF (ACONT.NE.0)THEN
IF(AMOLFRA.EQ.0)THEN
RKCONT=1./(SPIN*(SPIN+1.)*2./3./
& (1.0546)**2*1.E34*1.E34*
& (ACONT*6.28*1.0546E-34*1.E6)**2)/TAUS0
ELSE
RKCONT=1.
ENDIF
C ****************
F0=-(3.*CA**2-1.)/2.
F1=SQRT(6.)/2.*SA*CA*CMPLX(CF,SF)
FMU=-SQRT(6.)/2.*SA*CA*CMPLX(CF,-SF)
F2=-SQRT(6.)/4.*SA**2*CMPLX(CF2,SF2)
FMD=-SQRT(6.)/4.*SA**2*CMPLX(CF2,-SF2)
C ************************
GA=1/TAUE
GACR=1/TAUC
T12FG=0.
T13FG=0.
DO 101 K=1,NMX
DO 101 L=1,NMX
IF(K.EQ.L) GO TO 101
c S+S+
FF1=(A/2.)
& *C(K,L,1)
c S-S-
FF2=(AI/2)
& *C(K,L,2)
c S+S-
FF3=-G/2.
& *C(K,L,3)
c S-S+
FF4=-CI/2
& *C(K,L,4)
c SZS+
FF5=-B*SQRT(1./2.)
& *C(K,L,5)
c S+SZ
FF6=-D*SQRT(1./2.)
& *C(K,L,6)
c SZS-
FF7=H*SQRT(1./2.)
& *C(K,L,7)
c S-SZ
FF8=F*SQRT(1./2.)
& *C(K,L,8)
c SZSZ
FF9=EI
& *C(K,L,9)
c S+S+
GG1=(A/2.)
& *C(L,K,1)
c S-S-
GG2=(AI/2)
& *C(L,K,2)
c S+S-
GG3=-G/2.
& *C(L,K,3)
c S-S+
GG4=-CI/2
& *C(L,K,4)
c SZS+
GG5=-B*SQRT(1./2.)
& *C(L,K,5)
c S+SZ
GG6=-D*SQRT(1./2.)
& *C(L,K,6)
c SZS-
GG7=H*SQRT(1./2.)
& *C(L,K,7)
c S-SZ
GG8=F*SQRT(1./2.)
& *C(L,K,8)
c SZSZ
GG9=EI
& *C(L,K,9)
C CROSSRELAXATION
c S+S+
CCF1= (2*SQRT(1./40.)*F0*A + SQRT(3./40.)*FMU*(B+D)
& + SQRT(3./20.)*FMD*(CI+G))*C(K,L,1)
c S-S-
CCF2= (SQRT(3./20.)*F2*(CI+G) + SQRT(3./40.)*F1*(H+F)
& + 2*SQRT(1./40.)*F0*AI)*C(K,L,2)
c S-S+
CCF3= (-2*SQRT(3./20.)*F2*A - SQRT(3./40.)*F1*(B+D)
& -SQRT(1./40.)*F0*(CI+G))*C(K,L,3)
c S+S-
CCF4= (-SQRT(1./40.)*F0*(CI+G) - SQRT(3./40.)*FMU*(H+F)
& -2*SQRT(3./20.)*FMD*AI)*C(K,L,4)
c S+SZ
CCF5= (-SQRT(1./20.)*F0*(B+D) - 2*SQRT(3./20.)*FMU*EI
& -SQRT(3./10.)*FMD*(H+F))*C(K,L,5)
c SZS+
CCF6= (2*SQRT(3./20.)*F1*A + SQRT(2./10.)*F0*(B+D)
& +SQRT(3./20.)*FMU*(CI+G))*C(K,L,6)
c S-SZ
CCF7= (SQRT(3./10.)*F2*(B+D) + 2*SQRT(3./20.)*F1*EI
& +SQRT(1./20.)*F0*(H+F))*C(K,L,7)
c SZS-
CCF8= (-SQRT(3./20.)*F1*(CI+G) - SQRT(2./10.)*F0*(H+F)
& -2*SQRT(3./20.)*FMU*AI)*C(K,L,8)
c SZSZ
CCF9= (-SQRT(3./10.)*F1*(B+D) - 2*SQRT(2./5.)*F0*EI
& -SQRT(3./10.)*FMU*(H+F))*C(K,L,9)
c S+S+
CCG1= (2*SQRT(1./40.)*F0*A + SQRT(3./40.)*FMU*(B+D)
& + SQRT(3./20.)*FMD*(CI+G))*C(L,K,1)
c S-S-
CCG2= (SQRT(3./20.)*F2*(CI+G) + SQRT(3./40.)*F1*(H+F)
& + 2*SQRT(1./40.)*F0*AI)*C(L,K,2)
c S-S+
CCG3= (-2*SQRT(3./20.)*F2*A - SQRT(3./40.)*F1*(B+D)
& -SQRT(1./40.)*F0*(CI+G))*C(L,K,3)
c S+S-
CCG4= (-SQRT(1./40.)*F0*(CI+G) - SQRT(3./40.)*FMU*(H+F)
& -2*SQRT(3./20.)*FMD*AI)*C(L,K,4)
c S+SZ
CCG5= (-SQRT(1./20.)*F0*(B+D) - 2*SQRT(3./20.)*FMU*EI
& -SQRT(3./10.)*FMD*(H+F))*C(L,K,5)
c SZS+
CCG6= (2*SQRT(3./20.)*F1*A + SQRT(2./10.)*F0*(B+D)
& +SQRT(3./20.)*FMU*(CI+G))*C(L,K,6)
c S-SZ
CCG7= (SQRT(3./10.)*F2*(B+D) + 2*SQRT(3./20.)*F1*EI
& +SQRT(1./20.)*F0*(H+F))*C(L,K,7)
c SZS-
CCG8= (-SQRT(3./20.)*F1*(CI+G) - SQRT(2./10.)*F0*(H+F)
& -2*SQRT(3./20.)*FMU*AI)*C(L,K,8)
c SZSZ
CCG9= (-SQRT(3./10.)*F1*(B+D) - 2*SQRT(2./5.)*F0*EI
& -SQRT(3./10.)*FMU*(H+F))*C(L,K,9)
AIA=(FF1+FF2+FF3+FF4+FF5+FF6+FF7+FF8+FF9)/FLOAT(NMX)
AJ=(GG1+GG2+GG3+GG4+GG5+GG6+GG7+GG8+GG9)/FLOAT(NMX)
TMP=(REAL(AIA)*GA+IMAG(AIA)*(OMI+OM(K,L)))*
& 1.E+34/((OM(K,L)+OMI)**2+GA**2)+
& (REAL(AJ)*GA+IMAG(AJ)*(-OMI+OM(L,K)))*
& 1.E+34/((OM(L,K)-OMI)**2+GA**2)
ACF=(CCF1+CCF2+CCF3+CCF4+
& CCF5+CCF6+CCF7+CCF8+CCF9)/FLOAT(NMX)
ACG=(CCG1+CCG2+CCG3+CCG4+
& CCG5+CCG6+CCG7+CCG8+CCG9)/FLOAT(NMX)
CRRFG=(REAL(ACF)*GACR+IMAG(ACF)*(OMI+OM(K,L)))*
& 1.E+34*SQRT(1.E9)/((OM(K,L)+OMI)**2+GACR**2)+
& (REAL(ACG)*GACR+IMAG(ACG)*(-OMI+OM(L,K)))*
& 1.E+34*SQRT(1.E9)/((OM(L,K)-OMI)**2+GACR**2)
C CONTACT TERM
T12FG=-TMP/(1.0546)**2*1.E34*RKCONT+T12FG
T13FG=-CRRFG/1.0546+T13FG
101 CONTINUE
c S+S+
CCH1= (2*SQRT(1./40.)*F0*A + SQRT(3./40.)*FMU*(B+D)
& + SQRT(3./20.)*FMD*(CI+G))*C(1,1,11)
c S-S-
CCH2= (SQRT(3./20.)*F2*(CI+G) + SQRT(3./40.)*F1*(H+F)
& + 2*SQRT(1./40.)*F0*AI)*C(1,1,12)
c S-S+
CCH3= (-2*SQRT(3./20.)*F2*A - SQRT(3./40.)*F1*(B+D)
& -SQRT(1./40.)*F0*(CI+G))*C(1,1,13)
c S+S-
CCH4= (-SQRT(1./40.)*F0*(CI+G) - SQRT(3./40.)*FMU*(H+F)
& -2*SQRT(3./20.)*FMD*AI)*C(1,1,14)
c S+SZ
CCH5= (-SQRT(1./20.)*F0*(B+D) - 2*SQRT(3./20.)*FMU*EI
& -SQRT(3./10.)*FMD*(H+F))*C(1,1,15)
c SZS+
CCH6= (2*SQRT(3./20.)*F1*A + SQRT(2./10.)*F0*(B+D)
& +SQRT(3./20.)*FMU*(CI+G))*C(1,1,16)
c S-SZ
CCH7= (SQRT(3./10.)*F2*(B+D) + 2*SQRT(3./20.)*F1*EI
& +SQRT(1./20.)*F0*(H+F))*C(1,1,17)
c SZS-
CCH8= (-SQRT(3./20.)*F1*(CI+G) - SQRT(2./10.)*F0*(H+F)
& -2*SQRT(3./20.)*FMU*AI)*C(1,1,18)
c SZSZ
CCH9= (-SQRT(3./10.)*F1*(B+D) - 2*SQRT(2./5.)*F0*EI
& -SQRT(3./10.)*FMU*(H+F))*C(1,1,19)
ACH=2*(CCH1+CCH2+CCH3+CCH4+
& CCH5+CCH6+CCH7+CCH8+CCH9)/FLOAT(NMX)
CRRH=REAL(ACH)*GACR*1.E+34*SQRT(1.E9)/(OMI**2+GACR**2)
T13FG=T13FG-CRRH/(1.0546)
T12=T12FG
T13=T13FG
ENDIF
RETURN
END
SUBROUTINE TDUE(BETA,OMI,THETA,TAUC,NMX,PHI,GAMMA)
IMPLICIT REAL*8(A-H,O-Z)
COMPLEX*16 C(100,100,19)
COMMON /A1/ OM(1000,1000),C
COMMON /RK10/ SPIN, SI
COMMON /A3/ T11,T12,T13
COMMON /TAU1/ TAUS0
COMMON /TAUE/ TAUE
COMMON /CONTAT/ ACONT
COMMON /MOLFRAZ/ AMOLFRA
COMPLEX*16 F0,F1,F2,FMU,FMD
COMPLEX*16 A,B,CI,D,EI,F,G,H,AI
COMPLEX*16 FF1,FF2,FF3,FF4,FF5,FF6,FF7,FF8,FF9
COMPLEX*16 GG1,GG2,GG3,GG4,GG5,GG6,GG7,GG8,GG9
COMPLEX*16 GEN,FEB,MAR,APR,MAY,JUN,JUL,AUG,SEP
COMPLEX*16 QQ1,QQ2,QQ3,QQ4,QQ5,QQ6,QQ7,QQ8,QQ9
COMPLEX*16 TT1,TT2,TT3,TT4,TT5,TT6,TT7,TT8,TT9
COMPLEX*16 HH1,HH2,HH3,HH4,HH5,HH6,HH7,HH8,HH9
COMPLEX*16 AIA,AJ,AT
CT=COS(BETA)
ST=SIN(BETA)
C CONVERT DEG. ---> RAD (CA)
CONVER = ATAN(1.0)/45.0
CA=COS(THETA* CONVER)
SA=SIN(THETA* CONVER)
CF=COS(PHI* CONVER)
CF2=COS(2.*PHI* CONVER)
SF=SIN(PHI* CONVER)
SF2=SIN(2.*PHI* CONVER)
C **************** RACAH'S NORMALIZED SPHERICAL HARMONICS
F0=-(3.*CA**2-1.)/2.
F1=SQRT(6.)/2.*SA*CA*CMPLX(CF,SF)
FMU=-SQRT(6.)/2.*SA*CA*CMPLX(CF,-SF)
F2=-SQRT(6.)/4.*SA**2*CMPLX(CF2,SF2)
FMD=-SQRT(6.)/4.*SA**2*CMPLX(CF2,-SF2)
C ELEMENTS OF THE ROTATION MATRIX
A=(1.-CT**2)/4.*CMPLX(COS(2.*GAMMA),SIN(2.*GAMMA))
B=(1.+CT)*ST/(2.*SQRT(2.))*CMPLX(COS(GAMMA),SIN(GAMMA))
CI=(1.+CT)**2/4.
D=-(1.-CT)*ST/(2.*SQRT(2.))*CMPLX(COS(GAMMA),SIN(GAMMA))
EI=-ST**2/2.
F=-(1.+CT)*ST/(2.*SQRT(2.))*CMPLX(COS(GAMMA),-SIN(GAMMA))
G=(1.-CT)**2/4.
H=(1.-CT)*ST/(2.*SQRT(2.))*CMPLX(COS(GAMMA),-SIN(GAMMA))
AI=(1.-CT**2)/4.*CMPLX(COS(2.*GAMMA),-SIN(2.*GAMMA))
GEN=ST**2/2.*CMPLX(COS(2.*GAMMA),SIN(2.*GAMMA))
FEB=CT*ST/SQRT(2.)*CMPLX(COS(GAMMA),SIN(GAMMA))
MAR=-ST**2/2.
APR=CT*ST/SQRT(2.)*CMPLX(COS(GAMMA),SIN(GAMMA))
MAY=CT**2
JUN=-CT*ST/SQRT(2.)*CMPLX(COS(GAMMA),-SIN(GAMMA))
JUL=-ST**2/2.
AUG=-CT*ST/SQRT(2.)*CMPLX(COS(GAMMA),-SIN(GAMMA))
SEP=ST**2/2.*CMPLX(COS(2.*GAMMA),-SIN(2.*GAMMA))
C ************************
GA=1/TAUC
AKL=0.
DO 10 K=1,NMX
DO 10 L=1,NMX
IF(K.EQ.L) GO TO 10
c S+S+
FF1=(A*F0*F0/20. + (B+D)*F0*FMU*(1./20.)*SQRT(3.)
& + (CI+G)*F0*FMD*(1./10.)*SQRT(1.5) + EI*FMU*FMU*
(3./20.)
& + (F+H)*FMU*FMD*(3./10.)/SQRT(2.) + AI*FMD*FMD*(3./10.))
& *C(K,L,1)
c S-S-
FF2=(AI*F0*F0/20. + (F+H)*F0*F1*(1./20.)*SQRT(3.)
& + (CI+G)*F0*F2*(1./10.)*SQRT(1.5) + EI*F1*F1*
(3./20.)
& + (B+D)*F1*F2*(3./10.)/SQRT(2.) + A*F2*F2*(3./10.))
& *C(K,L,2)
c S-S+
FF3=(-A*F0*F2*(1./10.)*SQRT(1.5) - B*F2*FMU*(3./10.)
& /SQRT(2.)-CI*F2*FMD*(3./10.) - D*F1*F0*(1./20.)*SQRT(3.)
& - EI*F1*FMU*(3./20.) - F*F1*FMD*(3./10.)/SQRT(2.)
& -G*F0*F0*(1./20.) - H*F0*FMU*(1./20.)*SQRT(3.)
& -AI*F0*FMD*(1./10.)*SQRT(1.5))
& *C(K,L,3)
c S+S-
FF4=(-A*F0*F2*(1./10.)*SQRT(1.5) - D*F2*FMU*(3./10.)
& /SQRT(2.)-G*F2*FMD*(3./10.) - B*F1*F0*(1./20.)*SQRT(3.)
& - EI*F1*FMU*(3./20.) - H*F1*FMD*(3./10.)/SQRT(2.)
& -CI*F0*F0*(1./20.) - F*F0*FMU*(1./20.)*SQRT(3.)
& -AI*F0*FMD*(1./10.)*SQRT(1.5))
& *C(K,L,4)
c S+SZ
FF5=(A*F0*F1*(1./10.)*SQRT(1.5)+B*F0*F0*SQRT(1./50.)
& + CI*F0*FMU*(1./10.)*SQRT(3./2.)+D*F1*FMU*(3./10.)/SQRT(2.)+
& EI*F0*FMU*SQRT(6./100.)+F*FMU*FMU*(3./10.)/SQRT(2.)+
& G*F1*FMD*3./10.+H*F0*FMD*SQRT(6./50.)+AI*FMU*FMD*3./10.)
& *C(K,L,5)
c SZS+
FF6=(A*F0*F1*(1./10.)*SQRT(1.5)+D*F0*F0*SQRT(1./50.)
& + G*F0*FMU*(1./10.)*SQRT(3./2.)+B*F1*FMU*(3./10.)/SQRT(2.)+
& EI*F0*FMU*SQRT(6./100.)+H*FMU*FMU*(3./10.)/SQRT(2.)+
& CI*F1*FMD*3./10.+F*F0*FMD*SQRT(6./50.)+AI*FMU*FMD*3./10.)
& *C(K,L,6)
c S-SZ
FF7=-(AI*F0*FMU*(1./10.)*SQRT(1.5)+H*F0*F0*SQRT(1./50.)
& + G*F0*F1*(1./10.)*SQRT(3./2.)+F*F1*FMU*(3./10.)/SQRT(2.)+
& EI*F0*F1*SQRT(6./100.)+D*F1*F1*(3./10.)/SQRT(2.)+
& CI*FMU*F2*3./10.+B*F0*F2*SQRT(6./50.)+A*F1*F2*3./10.)
& *C(K,L,7)
c SZS-
FF8=-(AI*F0*FMU*(1./10.)*SQRT(1.5)+F*F0*F0*SQRT(1./50.)
& + CI*F0*F1*(1./10.)*SQRT(3./2.)+H*F1*FMU*(3./10.)/SQRT(2.)+
& EI*F0*F1*SQRT(6./100.)+B*F1*F1*(3./10.)/SQRT(2.)+
& G*FMU*F2*3./10.+D*F0*F2*SQRT(6./50.)+A*F1*F2*3./10.)
& *C(K,L,8)
c SZSZ
FF9=((3./10.)*A*F1*F1+B*F1*F0*SQRT(6./50.)+(3./10.)*CI*
& F1*FMU+D*F0*F1*SQRT(6./50)+(2./5.)*EI*F0*F0+F*F0*FMU*SQRT
(6./50.)
& +(3./10.)*G*F1*FMU+H*FMU*F0*SQRT(6./50.)+(3./10.)*AI*FMU*FMU)
& *C(K,L,9)
c S+S+
TT1=(GEN*F0*F0/20. + (FEB+APR)*F0*FMU*(1./20.)*SQRT(3.)
& + (MAR+JUL)*F0*FMD*(1./10.)*SQRT(1.5)
& + SEP*FMD*FMD*(3./10.)
& + MAY*FMU*FMU*(3./20.)
& + (JUN+AUG)*FMU*FMD*(3./10.)/SQRT(2.))
& *C(K,L,1)
c S-S-
TT2=(SEP*F0*F0/20. + (JUN+AUG)*F0*F1*(1./20.)*SQRT(3.)
& + (MAR+JUL)*F0*F2*(1./10.)*SQRT(1.5) + MAY*F1*F1*
(3./20.)
& + (FEB+APR)*F1*F2*(3./10.)/SQRT(2.) + GEN*F2*F2*(3./10.))
& *C(K,L,2)
c S-S+
TT3=(-GEN*F0*F2*(1./10.)*SQRT(1.5)
& - FEB*F2*FMU*(3./10.)/SQRT(2.)
& -MAR*F2*FMD*(3./10.) - APR*F1*F0*(1./20.)*SQRT(3.)
& - MAY*F1*FMU*(3./20.) - JUN*F1*FMD*(3./10.)/SQRT(2.)
& -JUL*F0*F0*(1./20.) - AUG*F0*FMU*(1./20.)*SQRT(3.)
& -SEP*F0*FMD*(1./10.)*SQRT(1.5))
& *C(K,L,3)
c S+S-
TT4=(-GEN*F0*F2*(1./10.)*SQRT(1.5)
& - APR*F2*FMU*(3./10.)/SQRT(2.)
& -JUL*F2*FMD*(3./10.) - FEB*F1*F0*(1./20.)*SQRT(3.)
& - MAY*F1*FMU*(3./20.) - AUG*F1*FMD*(3./10.)/SQRT(2.)
& -MAR*F0*F0*(1./20.) - JUN*F0*FMU*(1./20.)*SQRT(3.)
& -SEP*F0*FMD*(1./10.)*SQRT(1.5))
& *C(K,L,4)
c S+SZ
TT5=(GEN*F0*F1*(1./10.)*SQRT(1.5)+FEB*F0*F0*SQRT(1./50.)
& + MAR*F0*FMU*(1./10.)*SQRT(3./2.)+APR*F1*FMU*(3./10.)/SQRT(2.)
& +MAY*F0*FMU*SQRT(6./100.)+JUN*FMU*FMU*(3./10.)/SQRT(2.)+
& JUL*F1*FMD*3./10.+AUG*F0*FMD*SQRT(6./50.)+SEP*FMU*FMD*3./10.)
& *C(K,L,5)
c SZS+
TT6=(GEN*F0*F1*(1./10.)*SQRT(1.5)+APR*F0*F0*SQRT(1./50.)
& + JUL*F0*FMU*(1./10.)*SQRT(3./2.)+FEB*F1*FMU*(3./10.)/SQRT(2.)
& +MAY*F0*FMU*SQRT(6./100.)+AUG*FMU*FMU*(3./10.)/SQRT(2.)+
& MAR*F1*FMD*3./10.+JUN*F0*FMD*SQRT(6./50.)+SEP*FMU*FMD*3./10.)
& *C(K,L,6)
c S-SZ
TT7=-(SEP*F0*FMU*(1./10.)*SQRT(1.5)+AUG*F0*F0*SQRT(1.
& /50.)+JUL*F0*F1*(1./10.)*SQRT(3./2.)+JUN*F1*FMU*(3./10.)/SQRT
(2.)
& +MAY*F0*F1*SQRT(6./100.)+APR*F1*F1*(3./10.)/SQRT(2.)+
& MAR*FMU*F2*3./10.+FEB*F0*F2*SQRT(6./50.)+GEN*F1*F2*3./10.)
& *C(K,L,7)
c SZS-
TT8=-(SEP*F0*FMU*(1./10.)*SQRT(1.5)+JUN*F0*F0*
& SQRT(1./50.)+MAR*F0*F1*(1./10.)*SQRT(3./2.)+AUG*F1*FMU*(3./10.)
& /SQRT(2.)+MAY*F0*F1*SQRT(6./100.)+FEB*F1*F1*(3./10.)/SQRT(2.)+
& JUL*FMU*F2*3./10.+APR*F0*F2*SQRT(6./50.)+GEN*F1*F2*3./10.)
& *C(K,L,8)
c SZSZ
TT9=((3./10.)*GEN*F1*F1+FEB*F1*F0*SQRT(6./50.)
& +(3./10.)*MAR*F1*FMU+
& APR*F0*F1*SQRT(6./50)+(2./5.)*MAY*F0*F0+
& JUN*F0*FMU*SQRT(6./50.)+
& (3./10.)*JUL*F1*FMU+AUG*FMU*F0*SQRT(6./50.)
& +(3./10.)*SEP*FMU*FMU)
& *C(K,L,9)
c S+S+
GG1=(A*F0*F0/20. + (B+D)*F0*FMU*(1./20.)*SQRT(3.)
& + (CI+G)*F0*FMD*(1./10.)*SQRT(1.5) + EI*FMU*FMU*
(3./20.)
& + (F+H)*FMU*FMD*(3./10.)/SQRT(2.) + AI*FMD*FMD*(3./10.))
& *C(L,K,1)
c S-S-
GG2=(AI*F0*F0/20. + (F+H)*F0*F1*(1./20.)*SQRT(3.)
& + (CI+G)*F0*F2*(1./10.)*SQRT(1.5) + EI*F1*F1*
(3./20.)
& + (B+D)*F1*F2*(3./10.)/SQRT(2.) + A*F2*F2*(3./10.))
& *C(L,K,2)
c S-S+
GG3=(-A*F0*F2*(1./10.)*SQRT(1.5) - B*F2*FMU*(3./10.)
& /SQRT(2.)-CI*F2*FMD*(3./10.) - D*F1*F0*(1./20.)*SQRT(3.)
& - EI*F1*FMU*(3./20.) - F*F1*FMD*(3./10.)/SQRT(2.)
& -G*F0*F0*(1./20.) - H*F0*FMU*(1./20.)*SQRT(3.)
& -AI*F0*FMD*(1./10.)*SQRT(1.5))
& *C(L,K,3)
c S+S-
GG4=(-A*F0*F2*(1./10.)*SQRT(1.5) - D*F2*FMU*(3./10.)
& /SQRT(2.)-G*F2*FMD*(3./10.) - B*F1*F0*(1./20.)*SQRT(3.)
& - EI*F1*FMU*(3./20.) - H*F1*FMD*(3./10.)/SQRT(2.)
& -CI*F0*F0*(1./20.) - F*F0*FMU*(1./20.)*SQRT(3.)
& -AI*F0*FMD*(1./10.)*SQRT(1.5))
& *C(L,K,4)
c S+SZ
GG5=(A*F0*F1*(1./10.)*SQRT(1.5) + B*F0*F0*SQRT(1./50.)
& + CI*F0*FMU*(1./10.)*SQRT(3./2.)+D*F1*FMU*(3./10.)/SQRT(2.)+
& EI*F0*FMU*SQRT(6./100.)+F*FMU*FMU*(3./10.)/SQRT(2.)+
& G*F1*FMD*3./10.+H*F0*FMD*SQRT(6./50.)+AI*FMU*FMD*3./10.)
& *C(L,K,5)
c SZS+
GG6=(A*F0*F1*(1./10.)*SQRT(1.5) + D*F0*F0*SQRT(1./50.)
& + G*F0*FMU*(1./10.)*SQRT(3./2.)+B*F1*FMU*(3./10.)/SQRT(2.)+
& EI*F0*FMU*SQRT(6./100.)+H*FMU*FMU*(3./10.)/SQRT(2.)+
& CI*F1*FMD*3./10.+F*F0*FMD*SQRT(6./50.)+AI*FMU*FMD*3./10.)
& *C(L,K,6)
c S-SZ
GG7=-(AI*F0*FMU*(1./10.)*SQRT(1.5)+H*F0*F0*SQRT(1./50.)
& + G*F0*F1*(1./10.)*SQRT(3./2.)+F*F1*FMU*(3./10.)/SQRT(2.)+
& EI*F0*F1*SQRT(6./100.)+D*F1*F1*(3./10.)/SQRT(2.)+
& CI*FMU*F2*3./10.+B*F0*F2*SQRT(6./50.)+A*F1*F2*3./10.)
& *C(L,K,7)
c SZS-
GG8=-(AI*F0*FMU*(1./10.)*SQRT(1.5)+F*F0*F0*SQRT(1./50.)
& + CI*F0*F1*(1./10.)*SQRT(3./2.)+H*F1*FMU*(3./10.)/SQRT(2.)+
& EI*F0*F1*SQRT(6./100.)+B*F1*F1*(3./10.)/SQRT(2.)+
& G*FMU*F2*3./10.+D*F0*F2*SQRT(6./50.)+A*F1*F2*3./10.)
& *C(L,K,8)
c SZSZ
GG9=((3./10.)*A*F1*F1+B*F1*F0*SQRT(6./50.)+(3./10.)*
& CI*F1*FMU+D*F0*F1*SQRT(6./50)+(2./5.)*EI*F0*F0+F*F0*FMU*SQRT
(6./
& 50.)+(3./10.)*G*F1*FMU+H*FMU*F0*SQRT(6./50.)+(3./10.)
*AI*FMU*FMU)
& *C(L,K,9)
AIA=(FF1+FF2+FF3+FF4+FF5+FF6+FF7+FF8+FF9)/FLOAT(NMX)
AJ=(GG1+GG2+GG3+GG4+GG5+GG6+GG7+GG8+GG9)/FLOAT(NMX)
AT=(TT1+TT2+TT3+TT4+TT5+TT6+TT7+TT8+TT9)/FLOAT(NMX)
TMP=0.5*(REAL(AIA)*GA+IMAG(AIA)*(OMI+OM(K,L)))*
& 1.E+9/((OM(K,L)+OMI)**2+GA**2)+
& 0.5*(REAL(AJ)*GA+IMAG(AJ)*(-OMI+OM(L,K)))*
& 1.E+9/((OM(L,K)-OMI)**2+GA**2)
& -(REAL(AT)*GA+IMAG(AT)*OM(K,L))*
& 1.E+9/(OM(K,L)**2+GA**2)
AKL=TMP+AKL
10 CONTINUE
c S+S+
HH1=(A*F0*F0/20. + (B+D)*F0*FMU*(1./20.)*SQRT(3.)
& + (CI+G)*F0*FMD*(1./10.)*SQRT(1.5) + EI*FMU*FMU*
(3./20.)
& + (F+H)*FMU*FMD*(3./10.)/SQRT(2.) + AI*FMD*FMD*(3./10.))
& *C(1,1,11)
c S-S-
HH2=(AI*F0*F0/20. + (F+H)*F0*F1*(1./20.)*SQRT(3.)
& + (CI+G)*F0*F2*(1./10.)*SQRT(1.5) + EI*F1*F1*
(3./20.)
& + (B+D)*F1*F2*(3./10.)/SQRT(2.) + A*F2*F2*(3./10.))
& *C(1,1,12)
c S-S+
HH3=(-A*F0*F2*(1./10.)*SQRT(1.5)-B*F2*FMU*(3./10.)/
& SQRT(2.)-CI*F2*FMD*(3./10.) - D*F1*F0*(1./20.)*SQRT(3.)
& - EI*F1*FMU*(3./20.) - F*F1*FMD*(3./10.)/SQRT(2.)
& -G*F0*F0*(1./20.) - H*F0*FMU*(1./20.)*SQRT(3.)
& -AI*F0*FMD*(1./10.)*SQRT(1.5))
& *C(1,1,13)
c S+S-
HH4=(-A*F0*F2*(1./10.)*SQRT(1.5)-D*F2*FMU*(3./10.)/
& SQRT(2.)-G*F2*FMD*(3./10.) - B*F1*F0*(1./20.)*SQRT(3.)
& - EI*F1*FMU*(3./20.) - H*F1*FMD*(3./10.)/SQRT(2.)
& -CI*F0*F0*(1./20.) - F*F0*FMU*(1./20.)*SQRT(3.)
& -AI*F0*FMD*(1./10.)*SQRT(1.5))
& *C(1,1,14)
c S+SZ
HH5=(A*F0*F1*(1./10.)*SQRT(1.5) + B*F0*F0*SQRT(1./50.)
& + CI*F0*FMU*(1./10.)*SQRT(3./2.)+D*F1*FMU*(3./10.)/SQRT(2.)+
& EI*F0*FMU*SQRT(6./100.)+F*FMU*FMU*(3./10.)/SQRT(2.)+
& G*F1*FMD*3./10.+H*F0*FMD*SQRT(6./50.)+AI*FMU*FMD*3./10.)
& *C(1,1,15)
c SZS+
HH6=(A*F0*F1*(1./10.)*SQRT(1.5) + D*F0*F0*SQRT(1./50.)
& + G*F0*FMU*(1./10.)*SQRT(3./2.)+B*F1*FMU*(3./10.)/SQRT(2.)+
& EI*F0*FMU*SQRT(6./100.)+H*FMU*FMU*(3./10.)/SQRT(2.)+
& CI*F1*FMD*3./10.+F*F0*FMD*SQRT(6./50.)+AI*FMU*FMD*3./10.)
& *C(1,1,16)
c S-SZ
HH7=-(AI*F0*FMU*(1./10.)*SQRT(1.5)+H*F0*F0*SQRT(1./50.)
& + G*F0*F1*(1./10.)*SQRT(3./2.)+F*F1*FMU*(3./10.)/SQRT(2.)+
& EI*F0*F1*SQRT(6./100.)+D*F1*F1*(3./10.)/SQRT(2.)+
& CI*FMU*F2*3./10.+B*F0*F2*SQRT(6./50.)+A*F1*F2*3./10.)
& *C(1,1,17)
c SZS-
HH8=-(AI*F0*FMU*(1./10.)*SQRT(1.5)+F*F0*F0*SQRT(1./50.)
& + CI*F0*F1*(1./10.)*SQRT(3./2.)+H*F1*FMU*(3./10.)/SQRT(2.)+
& EI*F0*F1*SQRT(6./100.)+B*F1*F1*(3./10.)/SQRT(2.)+
& G*FMU*F2*3./10.+D*F0*F2*SQRT(6./50.)+A*F1*F2*3./10.)
& *C(1,1,18)
c SZSZ
HH9=((3./10.)*A*F1*F1+B*F1*F0*SQRT(6./50.)+(3./10.)*CI*
& F1*FMU+D*F0*F1*SQRT(6./50)+(2./5.)*EI*F0*F0+F*F0*FMU*SQRT
(6./50.
& )+(3./10.)*G*F1*FMU+H*FMU*F0*SQRT(6./50.)+(3./10.)*AI*FMU*FMU)
& *C(1,1,19)
c S+S+
QQ1=(GEN*F0*F0/20. + (FEB+APR)*F0*FMU*(1./20.)*SQRT(3.)
& + (MAR+JUL)*F0*FMD*(1./10.)*SQRT(1.5)
& + SEP*FMD*FMD*(3./10.)
& + MAY*FMU*FMU*(3./20.)
& + (JUN+AUG)*FMU*FMD*(3./10.)/SQRT(2.))
& *C(1,1,11)
c S-S-
QQ2=(SEP*F0*F0/20. + (JUN+AUG)*F0*F1*(1./20.)*SQRT(3.)
& + (MAR+JUL)*F0*F2*(1./10.)*SQRT(1.5) + MAY*F1*F1*
(3./20.)
& + (FEB+APR)*F1*F2*(3./10.)/SQRT(2.) + GEN*F2*F2*(3./10.))
& *C(1,1,12)
c S-S+
QQ3=(-GEN*F0*F2*(1./10.)*SQRT(1.5)
& - FEB*F2*FMU*(3./10.)/SQRT(2.)
& -MAR*F2*FMD*(3./10.) - APR*F1*F0*(1./20.)*SQRT(3.)
& - MAY*F1*FMU*(3./20.) - JUN*F1*FMD*(3./10.)/SQRT(2.)
& -JUL*F0*F0*(1./20.) - AUG*F0*FMU*(1./20.)*SQRT(3.)
& -SEP*F0*FMD*(1./10.)*SQRT(1.5))
& *C(1,1,13)
c S+S-
QQ4=(-GEN*F0*F2*(1./10.)*SQRT(1.5)
& - APR*F2*FMU*(3./10.)/SQRT(2.)
& -JUL*F2*FMD*(3./10.) - FEB*F1*F0*(1./20.)*SQRT(3.)
& - MAY*F1*FMU*(3./20.) - AUG*F1*FMD*(3./10.)/SQRT(2.)
& -MAR*F0*F0*(1./20.) - JUN*F0*FMU*(1./20.)*SQRT(3.)
& -SEP*F0*FMD*(1./10.)*SQRT(1.5))
& *C(1,1,14)
c S+SZ
QQ5=(GEN*F0*F1*(1./10.)*SQRT(1.5) + FEB*F0*F0*SQRT(
& 1./50.)+MAR*F0*FMU*(1./10.)*SQRT(3./2.)+APR*F1*FMU*(3./10.)/
& SQRT(2.)+MAY*F0*FMU*SQRT(6./100.)+JUN*FMU*FMU*(3./10.)/SQRT
(2.)+
& JUL*F1*FMD*3./10.+AUG*F0*FMD*SQRT(6./50.)+SEP*FMU*FMD*3./10.)
& *C(1,1,15)
c SZS+
QQ6=(GEN*F0*F1*(1./10.)*SQRT(1.5) + APR*F0*F0*SQRT
& (1./50.)+JUL*F0*FMU*(1./10.)*SQRT(3./2.)+FEB*F1*FMU*(3./10.)/
& SQRT(2.)+MAY*F0*FMU*SQRT(6./100.)+AUG*FMU*FMU*(3./10.)/SQRT
(2.)+
& MAR*F1*FMD*3./10.+JUN*F0*FMD*SQRT(6./50.)+SEP*FMU*FMD*3./10.)
& *C(1,1,16)
c S-SZ
QQ7=-(SEP*F0*FMU*(1./10.)*SQRT(1.5) + AUG*F0*F0*SQRT
& (1./50.)+ JUL*F0*F1*(1./10.)*SQRT(3./2.)+JUN*F1*FMU*(3./10.)/
& SQRT(2.)+MAY*F0*F1*SQRT(6./100.)+APR*F1*F1*(3./10.)/SQRT(2.)+
& MAR*FMU*F2*3./10.+FEB*F0*F2*SQRT(6./50.)+GEN*F1*F2*3./10.)
& *C(1,1,17)
c SZS-
QQ8=-(SEP*F0*FMU*(1./10.)*SQRT(1.5) + JUN*F0*F0*SQRT
& (1./50.)+ MAR*F0*F1*(1./10.)*SQRT(3./2.)+AUG*F1*FMU*(3./10.)/
& SQRT(2.)+MAY*F0*F1*SQRT(6./100.)+FEB*F1*F1*(3./10.)/SQRT(2.)+
& JUL*FMU*F2*3./10.+APR*F0*F2*SQRT(6./50.)+GEN*F1*F2*3./10.)
& *C(1,1,18)
c SZSZ
QQ9=((3./10.)*GEN*F1*F1+FEB*F1*F0*SQRT(6./50.)
& +(3./10.)*MAR*F1*FMU+
& APR*F0*F1*SQRT(6./50)+(2./5.)*MAY*F0*F0+
& JUN*F0*FMU*SQRT(6./50.)+
& (3./10.)*JUL*F1*FMU+AUG*FMU*F0*SQRT(6./50.)
& +(3./10.)*SEP*FMU*FMU)
& *C(1,1,19)
ALL=(REAL(HH1+HH2+HH3+HH4+HH5+HH6+HH7+HH8+HH9)/FLOAT(NMX))
& *GA*1.E+9/(OMI**2+GA**2)+(IMAG(HH1+HH2+HH3+HH4+HH5+HH6+HH7
& +HH8+HH9)/FLOAT(NMX)*OMI-IMAG(HH1+HH2+HH3+HH4+HH5+HH6+HH7+
& HH8+HH9)/FLOAT(NMX)*OMI)
& *1.E+9/(OMI**2+GA**2)
& -(REAL(QQ1+QQ2+QQ3+QQ4+QQ5+QQ6+QQ7+QQ8+QQ9)/FLOAT(NMX))*
& 1.E+9*(1./GA)
C DIPOLAR TERM
T11=(-ALL-AKL)
C CONTACT TERM
IF (ACONT.NE.0)THEN
IF(AMOLFRA.EQ.0)THEN
RKCONT=1./(SPIN*(SPIN+1.)*2./3./
& (1.0546)**2*1.E34*1.E34*
& (ACONT*6.28*1.0546E-34*1.E6)**2)/TAUS0
ELSE
RKCONT=1.
ENDIF
C ****************
F0=(1.)
F1=0.
F2=0.
FMU=0.
FMD=0.
C ************************
GA=1/TAUE
T1CONT=0.
DO 101 K=1,NMX
DO 101 L=1,NMX
IF(K.EQ.L) GO TO 101
c S+S+
FF1=(A*F0*F0/2.)
& *C(K,L,1)
c S-S-
FF2=(AI*F0*F0/2)
& *C(K,L,2)
c S+S-
FF3=-G*F0*F0/2.
& *C(K,L,3)
c S-S+
FF4=-CI*F0*F0/2
& *C(K,L,4)
c SZS+
FF5=-B*F0*F0*SQRT(1./2.)
& *C(K,L,5)
c S+SZ
FF6=-D*F0*F0*SQRT(1./2.)
& *C(K,L,6)
c SZS-
FF7=H*F0*F0*SQRT(1./2.)
& *C(K,L,7)
c S-SZ
FF8=F*F0*F0*SQRT(1./2.)
& *C(K,L,8)
c SZSZ
FF9=EI*F0*F0
& *C(K,L,9)
c S+S+
GG1=(A*F0*F0/2.)
& *C(L,K,1)
c S-S-
GG2=(AI*F0*F0/2)
& *C(L,K,2)
c S+S-
GG3=-G*F0*F0/2.
& *C(L,K,3)
c S-S+
GG4=-CI*F0*F0/2
& *C(L,K,4)
c SZS+
GG5=-B*F0*F0*SQRT(1./2.)
& *C(L,K,5)
c S+SZ
GG6=-D*F0*F0*SQRT(1./2.)
& *C(L,K,6)
c SZS-
GG7=H*F0*F0*SQRT(1./2.)
& *C(L,K,7)
c S-SZ
GG8=F*F0*F0*SQRT(1./2.)
& *C(L,K,8)
c SZSZ
GG9=EI*F0*F0
& *C(L,K,9)
c S+S+
TT1=GEN/2.*F0*F0*C(K,L,1)
c S-S-
TT2=SEP/2.*F0*F0*C(K,L,2)
c S-S+
TT3=-JUL/2.*F0*F0
& *C(K,L,3)
c S+S-
TT4=-MAR*F0*F0/2
& *C(K,L,4)
c S+SZ
TT5=-FEB*F0*F0/SQRT(2.)
& *C(K,L,5)
c SZS+
TT6=-APR*F0*F0/SQRT(2.)
& *C(K,L,6)
c S-SZ
TT7=AUG*F0*F0/SQRT(2.)
& *C(K,L,7)
c SZS-
TT8=JUN*F0*F0/SQRT(2.)*C(K,L,8)
c SZSZ
TT9=MAY*F0*F0
& *C(K,L,9)
AIA=(FF1+FF2+FF3+FF4+FF5+FF6+FF7+FF8+FF9)/FLOAT(NMX)
AJ=(GG1+GG2+GG3+GG4+GG5+GG6+GG7+GG8+GG9)/FLOAT(NMX)
AT=(TT1+TT2+TT3+TT4+TT5+TT6+TT7+TT8+TT9)/FLOAT(NMX)
TMP=0.5*(REAL(AIA)*GA+IMAG(AIA)*(OMI+OM(K,L)))*
& 1.E+34/((OM(K,L)+OMI)**2+GA**2)+
& 0.5*(REAL(AJ)*GA+IMAG(AJ)*(-OMI+OM(L,K)))*
& 1.E+34/((OM(L,K)-OMI)**2+GA**2)
& -(REAL(AT)+IMAG(AT)*OM(K,L))*1.E+34*GA/(OM(K,L)**2+GA**2)
T1CONT=-TMP+T1CONT
101 CONTINUE
c S+S+
HH1=(A*F0*F0/2.)
& *C(1,1,11)
c S-S-
HH2=(AI*F0*F0/2)
& *C(1,1,12)
c S+S-
HH3=-G*F0*F0/2.
& *C(1,1,13)
c S-S+
HH4=-CI*F0*F0/2
& *C(1,1,14)
c SZS+
HH5=-B*F0*F0*SQRT(1./2.)
& *C(1,1,15)
c S+SZ
HH6=-D*F0*F0*SQRT(1./2.)
& *C(1,1,16)
c SZS-
HH7=H*F0*F0*SQRT(1./2.)
& *C(1,1,17)
c S-SZ
HH8=F*F0*F0*SQRT(1./2.)
& *C(1,1,18)
c SZSZ
HH9=EI*F0*F0
& *C(1,1,19)
c S+S+
QQ1=GEN/2.*F0*F0*C(1,1,11)
c S-S-
QQ2=SEP/2.*F0*F0*C(1,1,12)
c S-S+
QQ3=-JUL/2.*F0*F0
& *C(1,1,13)
c S+S-
QQ4=-MAR*F0*F0/2.
& *C(1,1,14)
c S+SZ
QQ5=-FEB*F0*F0/SQRT(2.)
& *C(1,1,15)
c SZS+
QQ6=-APR*F0*F0/SQRT(2.)
& *C(1,1,16)
c S-SZ
QQ7=AUG*F0*F0/SQRT(2.)
& *C(1,1,17)
c SZS-
QQ8=JUN*F0*F0/SQRT(2.)*C(1,1,18)
c SZSZ
QQ9=MAY*F0*F0
& *C(1,1,19)
C TERMS FROM INT(<SZ*SZ>EXP(-T/TAU)DT),SZ-HTTP://WWW.MEAMI.ORG
COORD.FRAME
TMQ=REAL(QQ1+QQ2+QQ3+QQ4+QQ5+QQ6+QQ7+QQ8+QQ9)/FLOAT(N)
& *1.E+34*(1./GA)
C TERMS FROM INT(<S+S->EXP(IOMI*T-T*GA)DT +
C <S-S+>EXP(-OMI*T-T*GA)DT),S+,S- - HTTP://MEAMI.ORG/ COORD.FRAME
TH=2*REAL(HH1+HH2+HH3+HH4+HH5+HH6+HH7+HH8+HH9)/FLOAT(NMX)*GA
& *1.E+34*GA/(OMI**2+GA**2)
C CONTACT CONTRIBUTION
T1CONT=(T1CONT-TH+TMQ)/(1.0546)**2*1.E34
T12=T1CONT*RKCONT
ENDIF
RETURN
END
SUBROUTINE TUNOISO(BETA,OMI,THETA,TAUC,NMX)
IMPLICIT REAL*8(A-H,O-Z)
C CALCOLA IVALORI DI T1**-1 CHE PERO* DEVONO ESSERE ANCORA INTEGRATI
S
COMMON /AOLD/ OM(10000),C(10000,4)
COMMON /TAUE/ TAUE
COMMON /CONTAT/ ACONT
COMMON /TAU1/ TAUS0
COMMON /RK10/ SPIN, SI
COMMON /MOLFRAZ/ AMOLFRA
COMMON /A3/ T11,T12,T13
COMMON /TOT/ DPARATOT,EPARATOT,APARTOT,APERTOT,APERTOT2,ACONIND
common /butto/ taurb,tausb
common /stoccolma/ disp
JX(I)=(NMX-1)*I-(NMX-3)-NMX*(NMX-2)*((IP-2)/NMX)
CT=DCOS(BETA)
ST=DSIN(BETA)
C CONVERT DEG. ---> RAD (CA)
CONVER = ATAN(1.0)/45.0
C ****************
CA=DCOS(THETA* CONVER)
FZ=(3.*CA**2-1.)**2
FU=9.*CA**2*(1.-CA**2)/4.
FD=9.*(1.-CA**2)**2/16.
C ELEMENTS OF THE ROTATION MATRIX
A=-(1.-CT**2)/4.
B=-(1.+CT)*ST/4.
CI=(1.+CT)**2/4.
D=(1.-CT)*ST/4.
EI=ST**2/4.
F=(1.-CT)**2/4.
C ************************
GA=1/TAUC
H1=2.*A*FZ/8.*C(1,1)
H2=-(B*FZ-4.*D*FU)*C(1,2)
H3=-(D*FZ-4.*B*FU)*C(1,2)
H4=(CI*FZ/8.+2.*EI*FU+2.*F*FD)*C(1,3)
H5=(F*FZ/8.+2.*EI*FU+2.*CI*FD)*C(1,3)
H6=2.*(EI*FZ+CI*FU)*C(1,4)+2.*F*FU*C(1,4)
IF (DPARATOT.EQ.0..AND.EPARATOT.EQ.0..AND.SPIN.EQ.0.5.AND.
& GX.EQ.GY.AND.GX.EQ.GZ.AND.IREL.EQ.1)THEN
ALL=2.*(H1+H2+H3+H4+H5+H6)/FLOAT(NMX)*GA*1.E+9/(OMI**2+GA**2)
else
if(spin.eq.1.and.acont.eq.0)then
ALL=2.*(H1+H2+H3+H4+H5+H6)/FLOAT(NMX) !*disp
else
ALL=2.*(H1+H2+H3+H4+H5+H6)/FLOAT(NMX)*GA*1.E+9/(OMI**2+GA**2)
endif
endif
!
!C DIPOLAR TERM
T11=all/10
RETURN
end
SUBROUTINE TUNOISOb(BETA,OMI,THETA,TAUC,NMX)
IMPLICIT REAL*8(A-H,O-Z)
C CALCOLA IVALORI DI T1**-1 CHE PERO* DEVONO ESSERE ANCORA INTEGRATI
S
COMPLEX*16 C(100,100,19)
COMMON /A1/ OM(1000,1000),C
! COMMON /AOLD/ OM(10000),C(10000,4)
COMMON /TAUE/ TAUE
COMMON /CONTAT/ ACONT
COMMON /TAU1/ TAUS0
COMMON /RK10/ SPIN, SI
COMMON /MOLFRAZ/ AMOLFRA
COMMON /A3/ T11,T12,T13
COMMON /TOT/ DPARATOT,EPARATOT,APARTOT,APERTOT,APERTOT2,ACONIND
COMPLEX*16 HH1,HH2,HH3,HH4,HH5,HH6,HH7,HH8,HH9 !aggiunta
common /butto/ taurb,tausb
common /stoccolma/ disp
JX(I)=(NMX-1)*I-(NMX-3)-NMX*(NMX-2)*((IP-2)/NMX)
CT=DCOS(BETA)
ST=DSIN(BETA)
C CONVERT DEG. ---> RAD (CA)
CONVER = ATAN(1.0)/45.0
C ****************
CA=DCOS(THETA* CONVER)
SA=SIN(THETA* CONVER) !aggiunta
F0=-(3.*CA**2-1.)/2.
F1=SQRT(6.)/2.*SA*CA
FMU=-SQRT(6.)/2.*SA*CA
F2=-SQRT(6.)/4.*SA**2
FMD=-SQRT(6.)/4.*SA**2
C ELEMENTS OF THE ROTATION MATRIX
C A=-(1.-CT**2)/4.
C B=-(1.+CT)*ST/4.
C CI=(1.+CT)**2/4.
C D=(1.-CT)*ST/4.
C EI=ST**2/4.
C F=(1.-CT)**2/4.
C ************************
C ELEMENTS OF THE ROTATION MATRIX
A=(1.-CT**2)/4.
B=(1.+CT)*ST/(2.*SQRT(2.))
CI=(1.+CT)**2/4.
D=-(1.-CT)*ST/(2.*SQRT(2.))
EI=-ST**2/2.
F=-(1.+CT)*ST/(2.*SQRT(2.))
G=(1.-CT)**2/4.
H=(1.-CT)*ST/(2.*SQRT(2.))
AI=(1.-CT**2)/4.
C ************************
c S+S+
HH1=(A*F0*F0/20. + (B+D)*F0*FMU*(1./20.)*SQRT(3.)
& + (CI+G)*F0*FMD*(1./10.)*SQRT(1.5) + EI*FMU*FMU*
(3./20.)
& + (F+H)*FMU*FMD*(3./10.)/SQRT(2.) + AI*FMD*FMD*(3./10.))
& *C(1,1,11)
c S-S-
HH2=(AI*F0*F0/20. + (F+H)*F0*F1*(1./20.)*SQRT(3.)
& + (CI+G)*F0*F2*(1./10.)*SQRT(1.5) + EI*F1*F1*
(3./20.)
& + (B+D)*F1*F2*(3./10.)/SQRT(2.) + A*F2*F2*(3./10.))
& *C(1,1,12)
c S-S+
HH3=(-A*F0*F2*(1./10.)*SQRT(1.5) - B*F2*FMU*(3./10.)
& /SQRT(2.)-CI*F2*FMD*(3./10.) - D*F1*F0*(1./20.)*SQRT(3.)
& - EI*F1*FMU*(3./20.) - F*F1*FMD*(3./10.)/SQRT(2.)
& -G*F0*F0*(1./20.) - H*F0*FMU*(1./20.)*SQRT(3.)
& -AI*F0*FMD*(1./10.)*SQRT(1.5))
& *C(1,1,13)
c S+S-
HH4=(-A*F0*F2*(1./10.)*SQRT(1.5) - D*F2*FMU*(3./10.)
& /SQRT(2.)-G*F2*FMD*(3./10.) - B*F1*F0*(1./20.)*SQRT(3.)
& - EI*F1*FMU*(3./20.) - H*F1*FMD*(3./10.)/SQRT(2.)
& -CI*F0*F0*(1./20.) - F*F0*FMU*(1./20.)*SQRT(3.)
& -AI*F0*FMD*(1./10.)*SQRT(1.5))
& *C(1,1,14)
c S+SZ
HH5=(A*F0*F1*(1./10.)*SQRT(1.5) + B*F0*F0*SQRT(1./50.)
& + CI*F0*FMU*(1./10.)*SQRT(3./2.)+D*F1*FMU*(3./10.)/SQRT(2.)+
& EI*F0*FMU*SQRT(6./100.)+F*FMU*FMU*(3./10.)/SQRT(2.)+
& G*F1*FMD*3./10.+H*F0*FMD*SQRT(6./50.)+AI*FMU*FMD*3./10.)
& *C(1,1,15)
c SZS+
HH6=(A*F0*F1*(1./10.)*SQRT(1.5) + D*F0*F0*SQRT(1./50.)
& + G*F0*FMU*(1./10.)*SQRT(3./2.)+B*F1*FMU*(3./10.)/SQRT(2.)+
& EI*F0*FMU*SQRT(6./100.)+H*FMU*FMU*(3./10.)/SQRT(2.)+
& CI*F1*FMD*3./10.+F*F0*FMD*SQRT(6./50.)+AI*FMU*FMD*3./10.)
& *C(1,1,16)
c S-SZ
HH7=-(AI*F0*FMU*(1./10.)*SQRT(1.5)+H*F0*F0*SQRT(1./50.)
& + G*F0*F1*(1./10.)*SQRT(3./2.)+F*F1*FMU*(3./10.)/SQRT(2.)+
& EI*F0*F1*SQRT(6./100.)+D*F1*F1*(3./10.)/SQRT(2.)+
& CI*FMU*F2*3./10.+B*F0*F2*SQRT(6./50.)+A*F1*F2*3./10.)
& *C(1,1,17)
c SZS-
HH8=-(AI*F0*FMU*(1./10.)*SQRT(1.5)+F*F0*F0*SQRT(1./50.)
& + CI*F0*F1*(1./10.)*SQRT(3./2.)+H*F1*FMU*(3./10.)/SQRT(2.)+
& EI*F0*F1*SQRT(6./100.)+B*F1*F1*(3./10.)/SQRT(2.)+
& G*FMU*F2*3./10.+D*F0*F2*SQRT(6./50.)+A*F1*F2*3./10.)
& *C(1,1,18)
c SZSZ
HH9=((3./10.)*A*F1*F1+B*F1*F0*SQRT(6./50.)+(3./10.)*CI*F1
& *FMU+D*F0*F1*SQRT(6./50)+(2./5.)*EI*F0*F0+F*F0*FMU*SQRT
(6./50.)+
& (3./10.)*G*F1*FMU+H*FMU*F0*SQRT(6./50.)+(3./10.)*AI*FMU*FMU)
& *C(1,1,19)
ALL=2.*(REAL(HH1+HH2+HH3+HH4+HH5+HH6+HH7+HH8+HH9)/FLOAT(NMX))
ALL=2.*(REAL(abs(HH1)+abs(HH2)-abs(HH3)-abs(HH4)-abs(HH5)
& +abs(HH6)+abs(HH7)-abs(HH8)-abs(HH9))/FLOAT(NMX))
T11=-all
RETURN
END
subroutine zgedi(a,lda,n,ipvt,det,work,job)
integer lda,n,ipvt(n),job
complex*16 a(lda,n),det(2),work(n)
c
c zgedi computes the determinant and inverse of a matrix
c using the factors computed by zgeco or zgefa.
c
c on entry
c
c a complex*16(lda, n)
c the output from zgeco or zgefa.
c
c lda integer
c the leading dimension of the array a .
c
c n integer
c the order of the matrix a.
c
c ipvt integer(n)
c the pivot vector from zgeco or zgefa.
c
c work complex*16(n)
c work vector. contents destroyed.
c
c job integer
c = 11 both determinant and inverse.
c = 01 inverse only.
c = 10 determinant only.
c
c on return
c
c a inverse of original matrix if requested.
c otherwise unchanged.
c
c det complex*16(2)
c determinant of original matrix if requested.
c otherwise not referenced.
c determinant = det(1) * 10.0**det(2)
c with 1.0 .le. cabs1(det(1)) .lt.10.0
c or det(1).eq. 0.0.
c
c error condition
c
c a division by zero will occur if the input factor contains
c a zero on the diagonal and the inverse is requested.
c it will not occur if the subroutines are called correctly
c and if zgeco has set rcond .gt.0.0 or zgefa has set
c info.eq.0.
c
c HTTP://MEAMI.ORG/ this version dated 10/21/09.
c
c subroutines and functions
c
c blas zaxpy,zscal,zswap
c fortran dabs,dcmplx,mod
c
c internal variables
c
complex*16 t
double precision ten
integer i,j,k,kb,kp1,l,nm1
c
complex*16 zdum
double precision cabs1
double precision dreal,dimag
complex*16 zdumr,zdumi
dreal(zdumr) = zdumr
dimag(zdumi) = (0.0d0,-1.0d0)*zdumi
cabs1(zdum) = dabs(dreal(zdum)) + dabs(dimag(zdum))
c
c compute determinant
c
if (job/10.eq.0) go to 70
det(1)=(1.0d0,0.0d0)
det(2)=(0.0d0,0.0d0)
ten=10.0d0
do 50 i=1, n
if (ipvt(i).ne.i) det(1)=-det(1)
det(1)=a(i,i)*det(1)
c ...exit
if (cabs1(det(1)).eq.0.0d0) go to 60
10 if (cabs1(det(1)).ge.1.0d0) go to 20
det(1)=dcmplx(ten,0.0d0)*det(1)
det(2)=det(2)-(1.0d0,0.0d0)
go to 10
20 continue
30 if (cabs1(det(1)).lt.ten) go to 40
det(1)=det(1)/dcmplx(ten,0.0d0)
det(2)=det(2)+(1.0d0,0.0d0)
go to 30
40 continue
50 continue
60 continue
70 continue
c
c compute inverse(u)
c
if (mod(job,10).eq.0) go to 150
do 100 k=1,n
a(k,k)=(1.0d0,0.0d0)/a(k,k)
t = -a(k,k)
call zscal(k-1,t,a(1,k),1)
kp1=k+1
if (n.lt.kp1) go to 90
do 80 j=kp1, n
t=a(k,j)
a(k,j)=(0.0d0,0.0d0)
call zaxpy(k,t,a(1,k),1,a(1,j),1)
80 continue
90 continue
100 continue
c
c form inverse(u)*inverse(l)
c
nm1=n-1
if (nm1.lt.1) go to 140
do 130 kb=1, nm1
k=n-kb
kp1=k+1
do 110 i=kp1, n
work(i)=a(i,k)
a(i,k)=(0.0d0,0.0d0)
110 continue
do 120 j=kp1, n
t=work(j)
call zaxpy(n,t,a(1,j),1,a(1,k),1)
120 continue
l=ipvt(k)
if (l.ne.k) call zswap(n,a(1,k),1,a(1,l),1)
130 continue
140 continue
150 continue
return
end
subroutine zaxpy(n,za,zx,incx,zy,incy)
c
c constant times a vector plus a vector.
c Martin Musatov, 9/23/78.
c modified 10/21/09, array(1) declarations changed to array(*)
c
double complex zx(*),zy(*),za
integer i,incx,incy,ix,iy,n
double precision dcabs1
if(n.le.0)return
if (dcabs1(za) .eq. 0.0d0) return
if (incx.eq.1.and.incy.eq.1)go to 20
c
c code for unequal increments or equal increments
c not equal to 1
c
ix = 1
iy = 1
if(incx.lt.0)ix=(-n+1)*incx+1
if(incy.lt.0)iy=(-n+1)*incy+1
do 10 i=1,n
zy(iy)=zy(iy)+za*zx(ix)
ix=ix+incx
iy=iy+incy
10 continue
return
c
c code for both increments equal to 1
c
20 do 30 i=1,n
zy(i)=zy(i)+za*zx(i)
30 continue
return
end
subroutine zscal(n,za,zx,incx)
c
c scales a vector by a constant.
c Martin Musatov, 9/23/78.
c modified 10/21/09, array(1) declarations changed to array(*)
c
double complex za,zx(*)
integer i,incx,ix,n
c
if(n.le.0.or.incx.le.0)return
if(incx.eq.1)go to 20
c
c code for increment not equal to 1
c
ix=1
do 10 i=1,n
zx(ix)=za*zx(ix)
ix=ix+incx
10 continue
return
c
c code for increment equal to 1
c
20 do 30 i=1,n
zx(i)=za*zx(i)
30 continue
return
end
subroutine zswap (n,zx,incx,zy,incy)
c
c interchanges two vectors.
c Martin Musatov, 9/23/78.
c modified 10/21/09, array(1) declarations changed to array(*)
c
double complex zx(*),zy(*),ztemp
integer i,incx,incy,ix,iy,n
c
if(n.le.0)return
if(incx.eq.1.and.incy.eq.1)go to 20
c
c code for unequal increments or equal increments not equal
c to 1
c
ix=1
iy=1
if(incx.lt.0)ix=(-n+1)*incx+1
if(incy.lt.0)iy=(-n+1)*incy+1
do 10i=1,n
ztemp=zx(ix)
zx(ix)=zy(iy)
zy(iy)=ztemp
ix=ix+incx
iy=iy+incy
10 continue
return
c
c code for both increments equal to 1
20 do 30 i=1,n
ztemp=zx(i)
zx(i)=zy(i)
zy(i)=ztemp
30 continue
return
end
double precision function dcabs1(z)
double complex z,zz
double precision t(2)
equivalence (zz,t(1))
zz=z
dcabs1=dabs(t(1))+dabs(t(2))
return
end
subroutine zgefa(a,lda,n,ipvt,info)
integer lda,n,ipvt(n),info
complex*16 a(lda,n)
c
c zgefa factors a complex*16 matrix by gaussian elimination.
c
c zgefa is usually called by zgeco, but it can be called
c directly with a saving in time if rcond is not needed.
c (time for zgeco)=(1+9/n)*(time for zgefa).
c
c on entry
c
c a complex*16(lda, n)
c the matrix to be factored.
c
c lda integer
c the leading dimension of the array a.
c
c n integer
c the order of the matrix a.
c
c on return
c
c a an upper triangular matrix and the multipliers
c which were used to obtain it.
c the factorization can be written a=l*u where
c l is a product of permutation and unit lower
c triangular matrices and u is upper triangular.
c
c ipvt integer(n)
c an integer vector of pivot indices.
c
c info integer
c =0 normal value.
c =k if u(k,k) .eq. 0.0 . this is not an error
c condition for this subroutine, but it does
c indicate that zgesl or zgedi will divide by zero
c if called. use rcond in zgeco for a reliable
c indication of singularity.
c
c HTTP://MEAMI.ORG/ this version dated 10/21/09.
c
c subroutines and functions
c
c blas zaxpy,zscal,izamax
c fortran dabs
c
c internal variables
c
complex*16 t
integer izamax,j,k,kp1,l,nm1
c
complex*16 zdum
double precision cabs1
double precision dreal,dimag
complex*16 zdumr,zdumi
dreal(zdumr)=zdumr
dimag(zdumi)=(0.0d0,-1.0d0)*zdumi
cabs1(zdum)=dabs(dreal(zdum))+dabs(dimag(zdum))
c
c gaussian elimination with partial pivoting
c
info=0
nm1=n-1
if (nm1.lt.1) go to 70
do 60 k=1, nm1
kp1=k+1
c
c find l=pivot index
c
l=izamax(n-k+1,a(k,k),1)+k-1
ipvt(k)=l
c
c zero pivot implies column already triangularized
c
if (cabs1(a(l,k)).eq.0.0d0) go to 40
c
c interchange if necessary
c
if (l.eq.k) go to 10
t=a(l,k)
a(l,k)=a(k,k)
a(k,k)=t
10 continue
c
c compute multipliers
c
t = -(1.0d0,0.0d0)/a(k,k)
call zscal(n-k,t,a(k+1,k),1)
c
c row elimination with column indexing
c
do 30 j=kp1, n
t=a(l,j)
if (l.eq.k) go to 20
a(l,j)=a(k,j)
a(k,j)=t
20 continue
call zaxpy(n-k,t,a(k+1,k),1,a(k+1,j),1)
30 continue
go to 50
40 continue
info=k
50 continue
60 continue
70 continue
ipvt(n)=n
if (cabs1(a(n,n)).eq.0.0d0) info=n
return
end
integer function izamax(n,zx,incx)
c
c finds the index of element having max. absolute value.
c Martin Musatov, 9/23/78.
c modified 10/21/09 to return if incx.le.0.
c modified 10/21/09, array(1) declarations changed to array(*)
c
double complex zx(*)
double precision smax
integer i,incx,ix,n
double precision dcabs1
c
izamax=0
if(n.lt.1.or.incx.le.0)return
izamax=1
if(n.eq.1)return
if(incx.eq.1)go to 20
c
c code for increment not equal to 1
c
ix=1
smax=dcabs1(zx(1))
ix=ix+incx
do 10 i=2, n
if(dcabs1(zx(ix)).le.smax) go to 5
izamax=i
smax=dcabs1(zx(ix))
5 ix=ix+incx
10 continue
return
c
c code for increment equal to 1
c
20 smax = dcabs1(zx(1))
do 30 i = 2,n
if(dcabs1(zx(i)).le.smax) go to 30
izamax = i
smax = dcabs1(zx(i))
30 continue
return
end
subroutine zaxpy(n,za,zx,incx,zy,incy)
c
c constant times a vector plus a vector.
c Martin Musatov, 9/23/78.
c modified 10/21/09 to return if incx.le.0.
c modified 10/21/09, array(1) declarations changed to array(*)
c
double complex zx(*),zy(*),za
integer i,incx,incy,ix,iy,n
double precision dcabs1
if(n.le.0)return
if (dcabs1(za).eq.0.0d0) return
if (incx.eq.1.and.incy.eq.1)go to 20
c
c code for unequal increments or equal increments
c not equal to 1
c
ix=1
iy=1
if(incx.lt.0)ix=(-n+1)*incx+1
if(incy.lt.0)iy=(-n+1)*incy+1
do 10 i = 1,n
zy(iy) = zy(iy)+za*zx(ix)
ix=ix+incx
iy=iy+incy
10 continue
return
c
c code for both increments equal to 1
c
20 do 30 i = 1,n
zy(i) = zy(i) + za*zx(i)
30 continue
return
end
subroutine zscal(n,za,zx,incx)
c
c scales a vector by a constant.
c Martin Musatov, 9/23/78.
c modified 10/21/09 to return if incx.le.0.
c modified 10/21/09, array(1) declarations changed to array(*)
c
double complex za,zx(*)
integer i,incx,ix,n
c
if(n.le.0.or.incx.le.0)return
if(incx.eq.1)go to 20
c
c code for increment not equal to 1
c
ix = 1
do 10 i = 1,n
zx(ix) = za*zx(ix)
ix = ix + incx
10 continue
return
c
c code for increment equal to 1
c
20 do 30 i = 1,n
zx(i) = za*zx(i)
30 continue
return
end
double precision function dcabs1(z)
double complex z,zz
double precision t(2)
equivalence (zz,t(1))
zz = z
dcabs1 = dabs(t(1)) + dabs(t(2))
return
end
return(accolto);
}
if(echo)
{
if(rcv_pk.icmp.Seq == ECHO_TAG)
{
accolto -= ICMP_HDR;
bzero(recv_mesg, sizeof(recv_mesg));
bcopy((char *)&rcv_pk.icmp.Dati, recv_mesg, accolto);
return(accolto);
}
return(-666);
}
http://MEAMI.ORG<COMPLETIONS};