補正法
長くかかる場合は途中で読み込みを停止して下さい。
誤差の補正は、何を補正するかによってパラメーターが異なるが、方法は全て同じである。ズレを消すように、補正磁石の値を決めるのである。補正磁石の値を変えると、他の場所でパラメーターがどう変わるかは計算できる（応答行列）ので、補正磁石の効果は小さいとして、足し算で計算を行う。具体的にＣＯＤの場合を述べる。
ＣＯＤとは、Closed Orbit Distortion である。訳すと、閉軌道の歪み、である。中心軌道は偏向磁石だけから決まるが、４極の据え付け誤差によって偏向磁石成分が発生すると、中心軌道が理想的な場合からずれてしまう。それがＣＯＤである。ちなみに、粒子の位置は
ただし、
は観測点、
は蹴りを入れた点である。全てのステアリングの蹴り角を上手く決めて、モニターのところでのＣＯＤが消えるようにするのが補正法である。具体的には、
実例
例えば、水平方向に４極磁石の据え付け誤差を５０μｍ入れた場合、発生したCODは以下のようになる。
NUMBER OF USED EIGENVALUES =
132
E= 4850.129743 3919.300326
2073.256311 1739.623329 297.942806
266.999217
226.463567 169.618160 104.298664
89.958780
85.593049
77.900576 67.019480 57.812598
36.091776
27.143331
21.191711 18.225494
9.936966 9.525603
7.128797
6.994291 3.897269
3.878624 3.360552
3.356050
3.241471 3.233803
2.081978 2.078976
1.983052
1.980730 0.671318
0.668179 0.631202
0.627392
0.620506 0.618133
0.277303 0.273718
0.271901
0.269083 0.227722
0.226744 0.150444
0.150436
0.149604 0.149555
0.146123 0.145825
0.144599
0.144567 0.142225
0.141722 0.141521
0.141458
0.140782 0.140724
0.140678 0.140509
0.139623
0.139564 0.139537
0.139324 0.103216
0.102211
0.101520 0.099843
0.053093 0.052965
0.052416
0.052364 0.047281
0.047249 0.046082
0.046010
0.044602 0.044411
0.039973 0.039921
0.038055
0.037917 0.036815
0.036614 0.034199
0.034097
0.031691 0.031607
0.031418 0.031344
0.030636
0.030595 0.030214
0.030073 0.018561
0.018509
0.007920 0.007796
0.002951 0.002903
0.002414
0.002388 0.002302
0.002297 0.002282
0.002264
0.002191 0.002172
0.002141 0.002140
0.002055
0.002044 0.002019
0.002017 0.001946
0.001942
0.001934 0.001931
0.001876 0.001872
0.000933
0.000931 0.000926
0.000921 0.000712
0.000699
0.000667 0.000658
0.000000 0.000000
0.000000
0.000000
これらを全て用いて補正を行うと、補正後のＣＯＤは以下のようになり、スパイクが立つ。
補正前
MOD dx
dEX
rms(mm).53242 4.09162
wghted..00172 .0129
max(mm)2.05719 11.0786
at MON.66 MON.132
くちゃくちゃに補正
MOD dx
dEX
rms(mm).02746 .79119
wghted.7.039E-5 .00265
max(mm)-.05841 -2.1814
at MON.52 MON.66
蹴り角：
RMS OF ANGLE(mrad) = 0.08017 FOR132 DATA
MAX OF ANGLE(mrad) = 0.33981 AT THE 35-TH CORRECTOR
緩やかに補正
MOD dx
dEX
rms(mm).06406 .61908
wghted.1.662E-4 -7.00197
max(mm)-.17885 -1.83983
at MON.106 MON.66
蹴り角：
RMS OF ANGLE(mrad) = 0.01272 FOR132 DATA
MAX OF ANGLE(mrad) = 0.03469 AT THE109-TH CORRECTOR
くちゃくちゃに補正した方がよく直るが、蹴り角が大きく、直していない場所ではスパイク状になっている。蹴り角の最大値を比べると、緩やかに補正した方が約１０分の１になっているし、直した後のＣＯＤも緩やかになっている。一般に、発生するCODはリングのチューン程度のものであるため、その程度のものが直ってしまえば、後は残っていてもそれ程影響がない。逆に、スパイクは悪影響を及ぼすため、ダイナミックアパーチャーとしては、緩やかに直した方が稼げる。また、くちゃくちゃに補正すると、ＣＯＤは直っても分散関数が大きく発生してしまうことになり、いたちごっこになりかねない。その点も、緩やかに直せば分散関数も自然に治ってくれる方向に行くため、問題が生じない。
!+
! Dynamic Aperture Survey with ERROR
!
!-
ON LOG;
ON RFSW RADCOD RAD COD FLUC DAPERT;
FFS;
USE RING;
cell calc emit;
sextf = LINE["K2","SF"];
sextd = LINE["K2","SD"];
SAVE;
com=OpenWrite["fort.89"];
$FORM="15.10";
L = 0;
IE= 0;
Write[com,"repeat number
: ","0"];
Write[com,"horizontal tune of ring
: ",Twiss["NX","***"]/2/Pi];
Write[com,"vertical tune of ring
: ",Twiss["NY","***"]/2/Pi];
Write[com,"number of eigen value use
: ",IE];
Write[com,"number of constraint condition :
",L];
Write[com,"number of constraint BMPs below : "];
Write[com,"symbol of horizontal steering
: ","STRH "];
Write[com,"symbol of vertical steering
: ","STRV "];
Write[com,"symbol of monitor
: ","MON "];
Write[com,"direction of calculation
: ","X"];
Write[com,"output file name
: ","fort.893"];
i=0;
repeat LINE["LENGTH"];
i = i + 1;
ty=LINE["ELEMENT",i];
tty = ty // "
";
tty = tty[1,7];
$FORM="15.10";
PageWidth = 1000;
kyori[i]=Twiss["S",i];
tunex[i]=Twiss["NX",i]/2/Pi;
tuney[i]=Twiss["NY",i]/2/Pi;
betax[i]=Twiss["BX",i];
betay[i]=Twiss["BY",i];
Write[com,tty,i,kyori[i],L,L,L,tunex[i],tuney[i],betax[i],betay[i]];
until;
Close[com]
! System["cat fort.89"];
System["a.out"];
res = OpenRead["fort.893"];
IC = Read[res,Real];
Print["Number of correctors = ",IC];
i=0;
repeat IC;
i = i + 1;
{stnmpx[i], stvlpx[i]} = Read[res, {Real,Real}];
stvlx[i] = 0;
until;
Close[res];
SAVE;
USE RING;
cell calc emit fix;
optics=CalculateOptics[1,LINE["LENGTH"],Twiss['*',1],0,2];
options={Solver->Micado, SetSteer->True, Calc->True};
monitors=Monitor["MON*"];
steerings=Steer["STR*"];
correction = {'X','Y'};
iseed = 1188;
SeedRandom[iseed];
GCUT=2;
DELX 5E-5 Q*;
DELY 5E-5 Q*;
DELX 5E-5 SF*;
DELX 5E-5 SD*;
DELY 5E-5 SF*;
DELY 5E-5 SD*;
DELX 5E-5 B*
DK 5E-4 B*;
DK 5E-4 Q*;
DK 5E-4 S*;
DTHETA 2E-4 B*;
DTHETA 2E-4 Q*;
DTHETA 2E-4 S*;
GCUT = 1;
DELY 5E-5 SF*;
DELY 5E-5 SD*;
MON 0 0 0 0 MON*;
STE 0 STR*;
cell;
calc;
dump;
! disp MARKR;
codplot;
out 16 draw DX & DY MON* end;
out 12 draw DOX & DOY end;
MCOD;
TCOD;
J = 0;
repeat 15;
J = J + 1;
RENUM = J;
Print["How much sextpoles?"];
Print["Input % of sext"];
a = Read[6,Real];
IF a <=100
Print[a,"% of sextupole"];
LINE["K2","SF"] = sextf * a /100;
LINE["K2","SD"] = sextd * a /100;
sextprev = a;
ELSE
Print["not changed from previous value"];
ENDIF;
Write[$Output,"Witch direction do you want to calculate?"];
Print["X : 1 Y : 2"];
a = Read[6, Real];
Print["Input = ",a];
IF a == 1
Print["Direction is X !"];
DRCTN = "X";
ELSEIF a == 2
Print["Direction is Y !"];
DRCTN = "Y";
ELSE
Exit[];
ENDIF;
Print[DRCTN];
Print["How many eigenvalues used in calculation?"];
a = Read[6,Real];
IE = a;
Print["Number of Eigenvalues (Input) = ",a];
com=OpenWrite["fort.89"];
$FORM="15.10";
! a = 1;
! b = 64;
! c = 65;
! d = 128;
! L =4;
L = 0;
! IE=2;
Write[com,"repeat number
: ",RENUM];
Write[com,"horizontal tune of ring
: ",Twiss["NX","***"]/2/Pi];
Write[com,"vertical tune of ring
: ",Twiss["NY","***"]/2/Pi];
Write[com,"number of eigen value use
: ",IE];
Write[com,"number of constraint condition :
",L];
Write[com,"number of constraint BMPs below : "];
! Write[com,a];
! Write[com,b];
! Write[com,c];
! Write[com,d];
Write[com,"symbol of horizontal steering
: ","STRH "];
Write[com,"symbol of vertical steering
: ","STRV "];
Write[com,"symbol of monitor
: ","MON "];
Write[com,"direction of calculation
: ",DRCTN];
Write[com,"output file name
: ","fort.893"];
i=0;
repeat LINE["LENGTH"];
i = i + 1;
ty=LINE["ELEMENT",i];
tty = ty // "
";
tty = tty[1,7];
$FORM="15.10";
PageWidth = 1000;
Write[com,tty,i,kyori[i],LINE["K0",i]*1000,Twiss["DX",i]*1000,Twiss["DY",i]*1000,tunex[i],tuney[i],betax[i],betay[i]];
until;
Close[com]
System["a.out"];
res = OpenRead["fort.893"];
IC = Read[res,Real];
Print["Number of correctors = ",IC];
i=0;
repeat IC;
i = i + 1;
{stnm[i], stvl[i]} = Read[res, {Real,Real}];
LINE["K0",stnm[i]] = stvl[i]/1000;
until;
Close[res];
SAVE;
cell;
calc;
codplot;
TCOD;
MCOD;
Print["Do you want to save?"];
Print["1 : save 2 : not save"];
a = Read[6,Real];
IF a == 1;
i = 0;
repeat IC;
i = i +
1;
IF DRCTN
== "X"
stnmpx[i] = stnm[i];
stvlpx[i] = stvl[i];
ELSE
stnmpy[i] = stnm[i];
stvlpy[i] = stvl[i];
ENDIF;
until;
ELSE
i = 0;
repeat IC;
i = i +
1;
stnm[i]
= stnmpx[i];
stvl[i]
= stvlpx[i];
LINE["K0",stnm[i]]
= stvl[i]/1000;
stnm[i]
= stnmpy[i];
stvl[i]
= stvlpy[i];
LINE["K0",stnm[i]]
= stvl[i]/1000;
until;
SAVE;
cell calc;
codplot;
TCOD;
MCOD;
ENDIF;
!dump;
out 13;
draw BX & BY& EX & EY Q*
out 13 end;
out 11;
draw DX & DY MON*;
out 13 end;
out 12;
draw DOX & DOY MON*;
out 12 end;
until;
out 71;
draw DX & DY MON*;
out 71 end;
DAPWIDTH = 10;
dp = 0.01 / SIGE;
zrange = Range[-2.5 * dp, 2.5 * dp, dp /2];
range = {{0,200},{1,1},zrange};
cell;
calc;
disp MARKR;
res = DynamicApertureSurvey[range,1000,Output->6];
res;
range = {{1,1},{0,200},zrange};
cell;
calc;
res = DynamicApertureSurvey[range,1000,Output->6];
res;
stop;
stop;
!********************************************************************
!
! Diffraction & IBS Limited Ring
!
! v28_44 Tuesday,15,June,1999
!
6th plan, with BPM & STERRING
!
!********************************************************************
MOMENTUM=1.0 GEV
;
MARK
MARKI = ()
MARKR = ()
MARKN = ()
MARKNE = ()
MARKD = ()
MARKL = ()
MARKB = ()
;
CAVI CAVI=(VOLT=0.7 MV FREQ=500.1 MHZ)
;
BEND
BH =(L = 0.40
ANGLE = 7.5 DEG E1 = 0.5 E2 = 0.5 )
B =(L =
0.80 ANGLE = 15.0 DEG E1 = 0.5 E2 = 0.5 )
BB =(L = 0.40
ANGLE = 7.5 DEG E1 = 1.0 E2 = 0.0 )
;
DRIFT
D2 = (L = 0.2
)
D3 = (L = 0.25
)
D4 = (L = 0.2
)
D5 = (L = 0.53
)
DDD = (L = 1E-10 )
DA1L = (L = 0.2 )
DA2L = (L = 0.2 )
DA3L = (L = 0.53 )
DB1L =(L
=.2 )
DB2L =(L
=.3 )
DB3L =(L
=.2 )
DB4L =(L
=2.0 )
DB4L2 =(L =0.3)
DB5L =(L
=5.3)
DB6L =(L
=.35 )
DB7L =(L
=.35 )
DBLL =(L
=14.5+(0.377416292741-0.32)/4)
D2
=(L =.2 )
D3
=(L =.25 )
DJ1L =(L
=.2 )
DJ2L =(L
=.2 )
DJ3L =(L
=.2 )
DK1L =(L
=.2 )
DK2L =(L
=.2 )
DK3L =(L
=.7 )
DK4L =(L
=.2 )
DLL
=(L = 1 )
;
SEXT
SF
=(L =.2 K2 =29.8424804865968 )
SD
=(L =.2 K2 =-24.4423215023308 )
;
QUAD
QB7L
=(L =.4 K1 =.3639907665886 )
QB6L
=(L =.4 K1 =-.6472295378725 )
QB5L
=(L =.4 K1 =.3064281706477 )
QB4L
=(L =.4 K1 =-.0631552392035 )
QB3L
=(L =.4 K1 =-.7865764447298 )
QB2L
=(L =.4 K1 =1.7974790782121 )
QB1L
=(L =.3 K1 =-1.0683174753665 )
QA2L
=(L =.2 K1 =.6993394712946 )
QA1L
=(L =.2 K1 =.8922697427839 )
QD
=(L =.3 K1 =-.691457890618 )
QF
=(L =.4 K1 =1.5190240503441 )
QJ1L
=(L =.2 K1 =1.0629320220076 )
QJ2L
=(L =.3 K1 =.6013887503572 )
QK1L
=(L =.4 K1 =-.7740106904263 )
QK2L
=(L =.35 K1 =1.9234509112895 )
QK3L
=(L =.4 K1 =-1.7470706219873 )
QK4L
=(L =.35 K1 =1.4637391046 )
;
MONI
MON = ()
;
BEND
STRH = ( L = 0
)
STRV = ( L = 0
ROTATE = -90 DEG)
;
LINE
MONI = (DDD
MON DDD)
! MONI =
(MON)
RING24 = (MARKI 12*NCELL
CAVI 12*NCELL MARKR)
RING = (MARKI
-LCELLH 5*NCELL SCELLH
-SCELLH 5*NCELL LCELLH
CAVI
-LCELLH 5*NCELL SCELLH
-SCELLH 5*NCELL LCELLH
MARKR)
STER = (DDD
STRH DDD STRV DDD)
! STER =
(STRH STRV)
NCELL = (MARKN
STER SD STER D2 QD MONI D3 SF STER D4 QF MONI D5
B
D5 QF MONI D4 SF STER D3 QD MONI D2 STER SD
MARKNE
)
! NCELL = (MARKN
!
STER SD D2 QD MONI D3 SF STER D4 QF MONI D5
!
B
!
D5 QF MONI D4 SF STER D3 QD MONI D2 SD
!
MARKNE
!
)
LCELLH = (MARKN SD STER
D2 QD MONI D3 SF STER DA1L QA1L
DA2L QA2L MONI DA3L
BH
MARKD DB1L QB1L DB2L STER QB2L MONI DB3L STER QB3L MONI
DB4L DB4L2 STER QB4L MONI DB5L STER QB5L MONI
DB6L QB6L DB7L STER QB7L MONI STER
DBLL
MARKL
)
SCELLH = (MARKN SD MONI
STER D2 QD D3 MONI SF MONI STER DJ1L QJ1L
DJ2L QJ2L DJ3L
BH
MARKD DK1L QK1L DK2L QK2L MONI DK3L STER QK3L
DK4L STER QK4L MONI
DLL
MARKL
)
C
C
C PROGRAM TEVECO
C
C Simulation test of
Orbit Correction with Eigen Vector Mathod
C
AND
C
Really COD correction of VSX New Ring with SAD
C
C
Coded by N. NAKAMURA (11/11/1993)
C
Coded by M. Sato (June, 1998)
C
C
C
Definition of Global Variable
IMPLICIT REAL*8 (A-H,O-Z)
REAL*8 Z(150),ZL(150),TNX,TNY,TN,PI,PI2
REAL*8 BEM(150),PHIM(150),BEC(150),PHIC(150)
REAL*8 E(150),V(150,150),C(150,150),CC(150,150),D(150)
REAL*8 EP(150),VP(150,150)
REAL*8 X(150),Y(150),DX(150),DY(150),Z0(150),ZW(150)
INTEGER IM,IC,IE,LP(150)
INTEGER NSADM(150), NSADC(150),
RENUM
CHARACTER CBMP(150)*11,CCOR(150)*11
CHARACTER*1 DIRECT
CHARACTER NAME*12, NCORX*7, NCORY*7,
NMONI*7, COMMENT*35
CHARACTER OUTFILE*12
REAL*8 S, ANGLE, DELX, DELY, TUNEX,
TUNEY
REAL*8 BETAX, BETAY, NUMSAD, CSTNB,
CSTPT, EVNB
C
C Reading command and optics data file
from SAD
C File name may be "fort.89",
please
C
OPEN(UNIT=89,FILE='fort.89',STATUS='OLD',FORM='FORMATTED')
C
C Readinf Header and Commands
C
READ(89,404) COMMENT, A
RENUM = INT(A)
Write(6,*) "repeat", RENUM
READ(89,404) COMMENT, TUNEX
Write(6,404) COMMENT, TUNEX
READ(89,404) COMMENT, TUNEY
Write(6,404) COMMENT, TUNEY
TNX = TUNEX
TNY = TUNEY
READ(89,404) COMMENT, EVNB
Write(6,404) COMMENT, EVNB
IE = INT(EVNB)
C IF (IE. EQ. 0) GO TO 999
READ(89,404) COMMENT, CSTNB
Write(6,404) COMMENT, CSTNB
L = INT(CSTNB)
READ(89,406) COMMENT
Write(6,406) COMMENT
DO 504 i=1, L
READ(89,405) CSTPT
LP(i) = INT(CSTPT)
Write(6,*) LP(i)
504 CONTINUE
READ(89,407) COMMENT, NCORX
Write(6,407) COMMENT, NCORX
READ(89,407) COMMENT, NCORY
Write(6,407) COMMENT, NCORY
READ(89,407) COMMENT, NMONI
Write(6,407) COMMENT, NMONI
READ(89,408) COMMENT, DIRECT
Write(6,408) COMMENT, DIRECT
IF (DIRECT. EQ. "X") THEN
IDIR = 1
ELSEIF (DIRECT. EQ. "Y") THEN
IDIR = 2
ELSE
Write(6,*) "ERROR!"
GO TO 999
ENDIF
READ(89,415) COMMENT, OUTFILE
Write(6,415) COMMENT, OUTFILE
503 CONTINUE
404 FORMAT(A35,1F15.9)
401 FORMAT(A7,9F15.9)
403 FORMAT(A7,1F5.0, 1F8.2,1G12.3,7F9.4)
405 FORMAT(1F15.9)
406 FORMAT(A35)
407 FORMAT(A35,A7)
408 FORMAT(A35,A1)
409 FORMAT(A35,I10)
415 FORMAT(A35,A12)
PI = Acos(-1.0)
PI2 = PI * 2.D0
C
C Reading Optics
C
M = 0
N = 0
DO 501 I=1,10000
READ(89,401,END=502) NAME, NUMSAD,
S, ANGLE, DELX,
*DELY, TUNEX, TUNEY, BETAX, BETAY
IF (IDIR. EQ. 1. AND. NAME. EQ. NCORX)
THEN
C WRITE(6,403) NAME, NUMSAD, S, ANGLE,
DELX, DELY,
C *TUNEX, TUNEY, BETAX, BETAY
M = M + 1
TN
= TNX
CCOR(M) = NAME
NSADC(M) = INT(NUMSAD)
DX(M) = ANGLE
DY(M) = 0
PHIC(M) = TUNEX * PI2 /TN
BEC(M) = BETAX
ELSEIF (IDIR. EQ. 2. AND. NAME. EQ.
NCORY) THEN
C WRITE(6,403)
NAME, NUMSAD, S, ANGLE, DELX, DELY,
C *TUNEX, TUNEY, BETAX, BETAY
M = M + 1
TN
= TNY
CCOR(M) = NAME
NSADC(M) = INT(NUMSAD)
DX(M) = 0
DY(M) = ANGLE
PHIC(M) = TUNEY * PI2 /TN
BEC(M) = BETAY
ELSEIF (NAME. EQ. NMONI) THEN
C WRITE(6,403)
NAME, NUMSAD, S, ANGLE, DELX, DELY,
C *TUNEX, TUNEY, BETAX, BETAY
N = N + 1
CBMP(N) = NAME
NSADM(N) = INT(NUMSAD)
X(N) = DELX
Y(N) = DELY
IF(IDIR. EQ. 1) THEN
BEM(N)
= BETAX
PHIM(N)
= TUNEX * PI2 / TN
Z0(N)
= X(N)
ZW(N)
= X(N)
ELSEIF(IDIR. EQ. 2)
THEN
BEM(N)
= BETAY
PHIM(N)
= TUNEY * PI2 / TN
Z0(N)
= Y(N)
ZW(N)
= Y(N)
END IF
END IF
IM = N
IC = M
501 CONTINUE
502 CONTINUE
CLOSE(89)
C IMH=INT(IM/2)
C ICH=INT(IC/2)
C IMH2=INT((IM+1)/2)
C ICH2=INT((IC+1)/2)
C WRITE(6,*)' '
C WRITE(6,*)'** BPM OPTICS PARAMETERS
**'
C WRITE(6,*) ' NO. BETA FUNC.
PHASE NO. BETA FUNC. PHASE'
C IF(IMH.NE.IMH2)THEN
C
C DO 70 I=1,IMH
C
WRITE(6,202) NSADM(I),BEM(I),PHIM(I),
C *NSADM(I+IMH2),BEM(I+IMH2),PHIM(I+IMH2)
C 70
CONTINUE
C
C WRITE(6,212) NSADM(IMH2),BEM(IMH2),PHIM(IMH2)
C ELSE
C
C DO 71 I=1,IMH
C
WRITE(6,202) NSADM(I),BEM(I),PHIM(I)
C *,NSADM(I+IMH),BEM(I+IMH),PHIM(I+IMH)
C 71
CONTINUE
C
C ENDIF
C
C WRITE(6,*)' '
C WRITE(6,*)'** CORRECTOR OPTICS
PARAMETERS **'
C WRITE(6,*) ' NO. BETA FUNC.
PHASE NO. BETA FUNC. PHASE'
C IF(ICH.NE.ICH2)THEN
C
C DO 72 I=1,ICH
C
WRITE(6,202) NSADC(I),BEC(I),PHIC(I),
C *NSADC(I+ICH2),BEC(I+ICH2),PHIC(I+ICH2)
C 72 CONTINUE
C
C WRITE(6,212) NSADC(ICH2),BEC(ICH2),PHIC(ICH2)
C ELSE
C
C DO 73 I=1,ICH
C
WRITE(6,202) NSADC(I),BEC(I),PHIC(I),
C *NSADC(I+ICH),BEC(I+ICH),PHIC(I+ICH)
C 73 CONTINUE
C
C ENDIF
C
C 202 FORMAT(1H ,I4,X,2F10.5,3X,I4,X,2F10.5)
C 212 FORMAT(1H ,I4,X,2F10.5)
80 CONTINUE
SZ=0.D0
SZW=0.D0
ZWM1=ABS(ZW(1))
DO 86 I=1,IM
IF(ZWM1.LT.ABS(ZW(I)))THEN
ZWM1=ABS(ZW(I))
IZWM1=I
ENDIF
SZW=SZW+ZW(I)*ZW(I)
86 CONTINUE
IF(SZW. EQ. 0)THEN
Write(6,*) "RRRAAAANNNNNNDDDDDDAAAAAMMMMM"
DO 505 J=1, RENUM
DO
507 I=1, IC
D(I)
= (RAND(RENUM)*2-1)/10
IF(IDIR.
EQ. 1) DX(I) = 0
IF(IDIR.
EQ. 2) DY(I) = 0
C Write(6,*)
D(I)
507 CONTINUE
505 CONTINUE
GOTO 506
ENDIF
RMSZW1=SQRT(SZW/IM)
C WRITE(6,*) " "
C WRITE(6,*)'NUMBER OF USED EIGENVALUES
(0:RETURN TO THE 1-ST INPUT
C *MENU) ='
C write(6,*) IE
CALL EVECO(ZW,BEM,PHIM,IM,TN,BEC,PHIC,IC,D,C,E,IE,V,
*CC,L,LP,ZL,EP,VP)
WRITE(6,*)' '
WRITE(6,*)'NUMBER OF USED EIGENVALUES
=',IE
WRITE(6,203)'E=',(E(J),J=1,IE)
203 FORMAT(1H ,A2,5F14.6/,19(1H ,2X,5F14.6/))
SUMD2=0.0D0
DMAX=ABS(D(1))
DO 95 I=1,IC
SUMD2=SUMD2+D(I)*D(I)
IF(DMAX.LT.ABS(D(I)))THEN
DMAX=ABS(D(I))
IDMAX=I
ENDIF
95 CONTINUE
RMSD=SQRT(SUMD2/IC)
C
DO 104 I=1,IM
SUM=0.0D0
DO 106 J=1,IC
SUM=SUM+0.5D0*SQRT(BEM(I)*BEC(J))*D(J)*
*
COS(TN*(PI-ABS(PHIM(I)-PHIC(J))))/SIN(PI*TN)
106 CONTINUE
ZW(I)=ZW(I)+SUM
104 CONTINUE
C
SZW=0.0D0
ZWM2=ABS(ZW(1))
DO 112 I=1,IM
IF(ZWM2.LT.ABS(ZW(I)))THEN
ZWM2=ABS(ZW(I))
IZWM2=I
ENDIF
SZW=SZW+ZW(I)*ZW(I)
112 CONTINUE
RMSZW2=SQRT(SZW/IM)
C
IMH=IM/2
CCC
C WRITE(6,*)' '
C WRITE(6,*)'BEAM POSITIONS AT ALL
THE BPMS(mm)'
C WRITE(6,*)'(NOTE:* = USED BPM FOR
ORBIT CORRECTION)'
C Write(6,*) "No. is SAD element
number"
C WRITE(6,412) "No.","Uncorrected","Corrected","No.",
C *"Uncorrected","Corrected"
C 412 FORMAT(X,A3,X,A11,3X,A9,4X,A3,X,A11,2X,A9)
C
C DO 114 I=1,IMH
C Write(6,207) NSADM(I),Z0(I),ZW(I),NSADM(I+IMH)
C *,Z0(I+IMH),ZW(I+IMH)
C 114 CONTINUE
C
C 207 FORMAT(1H,X,I4,F10.5,3X,F10.5,3X,X,I4,F10.5,3X,F10.5)
C IR=MOD(IC,4)
C IF(IR.EQ.0)ICC=IC
C IF(IR.NE.0)ICC=IC-IR
C WRITE(6,*)' '
C WRITE(6,*)'DEFLECTION ANGLE(mrad)'
C Write(6,*) "No. is SAD element
number"
C WRITE(6,411) "No.","Angle","No.","Angle","No.","Angle"
C *,"No.","Angle"
C 411 FORMAT(5X,A3,3X,A5,5X,A3,3X,A5,5X,A3,3X,A5,5X,A3,3X,A5)
C
C DO 118 I=1,ICC,4
C WRITE(6,413) NSADC(I),D(I),NSADC(I+1),D(I+1),NSADC(I+2),D(I+2),
C *NSADC(I+3),D(I+3)
C 413 FORMAT(1H ,4(I5,F10.5,3X))
C 118 CONTINUE
C IF(IR.EQ.1)THEN
C WRITE(6,217) NSADC(ICC+1),D(ICC+1)
C ELSEIF(IR.EQ.2)THEN
C WRITE(6,217) NSADC(ICC+1),D(ICC+1),NSADC(ICC+2),D(ICC+2)
C ELSEIF(IR.EQ.3)THEN
C WRITE(6,217) NSADC(ICC+1),D(ICC+1),NSADC(ICC+2),D(ICC+2),
C *NSADC(ICC+3),D(ICC+3)
C ENDIF
C 227 FORMAT(1H ,I3,X,F10.5)
C 217 FORMAT(1H ,4(I4,F10.5,3X))
WRITE(6,*)' '
WRITE(6,*)'** UNCORRECTED ORBIT
**'
WRITE(6,237)'RMS OF POS. AT THE
BPMS(mm) =',RMSZW1,
*' FOR',IM,' DATA'
WRITE(6,237)'MAX OF POS. AT THE
BPMS(mm) =',ZWM1,' AT THE',
*IZWM1,'-TH BPM'
237 FORMAT(1H ,A,F10.5,A,I3,A)
247 FORMAT(1H ,A,F10.5,A,I3)
WRITE(6,*)' '
WRITE(6,*)'** CORRECTED ORBIT **'
WRITE(6,237)'RMS OF POS. AT THE
BPMS(mm) =',RMSZW2,
*' FOR',IM,' DATA'
WRITE(6,237)'MAX OF POS. AT THE
BPMS(mm) =',ZWM2,' AT THE',
*IZWM2,'-TH BPM'
RATRW=RMSZW2/RMSZW1
RATMW=ZWM2/ZWM1
WRITE(6,*)' '
WRITE(6,*)'** RATIO **'
WRITE(6,237)'RMS RATIO
=',RATRW,
*' FOR',IM,' DATA'
WRITE(6,237)'MAX RATIO
=',RATMW,
*' FOR',IM,' DATA'
WRITE(6,*)'** DEFLECTION ANGLE **'
WRITE(6,237)'RMS OF ANGLE(mrad)
=',RMSD,
*' FOR',IC,' DATA'
WRITE(6,237)'MAX OF ANGLE(mrad)
=',DMAX,' AT THE',IDMAX,
*'-TH CORRECTOR'
Write(6,*) " "
Write(6,*) " "
506 CONTINUE
C Write(6,*) "** Deflection Angle
(result) ** "
C Write(6,*) "
No. is SAD element number"
C Write(6,416) "No.","OLD value","New
value","Total value"
C 416 FORMAT(X,A3,9X,A9,8X,A9,7X,A11)
IF(IDIR.EQ.1)THEN
DO 130 I=1,IC
C
Write(6,414) "No.",NSADC(I),DX(I)
C *,"-",D(I),"=",DX(I)-D(I)
DX(I)=DX(I)-D(I)
414 FORMAT(X,A3,1I5,X,1F14.10,X,A,X,1F14.10,X,A,X,1F14.10)
130 CONTINUE
ELSE
DO 135 I=1,IC
C
Write(6,414) "No.",NSADC(I),DY(I)
C *,"-",D(I),"=",DY(I)-D(I)
DY(I)=DY(I)-D(I)
135 CONTINUE
ENDIF
OPEN(UNIT=20,FILE=OUTFILE,STATUS='UNKNOWN',FORM='FORMATTED')
Write(20,*) IC
IF(IDIR.EQ.1)THEN
DO 140 I=1,IC
WRITE(20,206) NSADC(I),DX(I)
140 CONTINUE
ELSEIF(IDIR.EQ.2) THEN
DO 150 I=1,IC
WRITE(20,206) NSADC(I),DY(I)
150 CONTINUE
END IF
206 FORMAT(I5,2X,F18.14)
CLOSE(20)
GO TO 999
999 CONTINUE
STOP
END
C *************************************************************
C *************************************************************
SUBROUTINE EVECO(Z,BEM,PHIM,IM,TN,BEC,PHIC,IC,D,C,E,IE,V,CC,L,LP,
*ZL,EP,VP)
C *************************************************************
C *************************************************************
C
C PURPOSE
C COMPUTE DEFLECTION
ANGLES BY EIGENVECTOR METHOD.
C DESCRIPTION OF PARAMETERS
C *INPUT*
C Z
- CLOSED ORBIT DISTORTION AT BPMS
C BEM - BETATRON
FUNCTION AT BPMS
C PHIM - BETATRON PHASE
AT BPMS
C IM - NUMBER
OF BPMS
C TN - BETATRON
TUNE
C BEC - BETATRON
FUNCTION AT CORRECTORS
C PHIC - BETATRON PHASE
AT CORRECTORS
C IC - NUMBER
OF CORRECTORS
C *OUTPUT*
C D
- DEFLECTION ANGLES OF CORRECTORS
C C
- CONVERSION MATRIX FROM POSITIONS(ETAS) TO ANGLES
C E
- EIGENVALUES OF MATRIX AA
C V
- EIGENVECTOR-MATRIX OF MATRIX AA
C IE - NUMBER
OF USED EIGENVALUES TO CALCULATE D & C
C EP
- EIGENVALUES OF MATRIX P
C VP
- EIGENVECTOR-MATRIX OF MATRIX P
C
C *OTHERS*
C ETA - ETA=Z/SQRT(BEM)
C A
- MATRIX EXPRESSING RELATION BETWEEN DEFLECTION ANGLES
C
AND POSITIONS (ETA=A*D)
C ATRN - TRANSPOSED MATRIX
OF A
C AA - AA=ATRN*A
C
C
CODED BY N. NAKAMURA (04/11/1993)
C
IMPLICIT REAL*8 (A-H,O-Z)
REAL*8 ETA(150),PI
REAL*8 Z(150), BEM(150),PHIM(150),BEC(150),PHIC(150)
REAL*8 D(150),C(150,150),CC(150,150),E(150),V(150,150),W(150,150)
REAL*8 EP(150),VP(150,150),P(150,150),WP(150,150)
REAL*8 QQ(150,150),Q(150,150),QCC(150,150),QCCW(150,150)
REAL*8 DM1(150,150),DM2(150,150),DM(150,150)
REAL*8 DV1(150),DV2(150),ZL(150),PWP(150,150)
REAL*8 WA(150,150),XX(150,150),XXX(150,150),XXP(150,150)
REAL*8 TN
INTEGER IM,IC,IE, L, LP(150)
REAL*8 A(150,150),ATRN(150,150),AA(150,150),EINV(150)
REAL*8 PP(150,150),EPINV(150),CCTRN(150,150)
REAL*8 OK1(150),OK2(150)
C
PI=ACOS(-1.0)
IF(IE.GT.IC)IE=IC
C SET ETA-VECTOR & A-MATRIX.
DO 10 I=1,IM
ETA(I)=Z(I)/SQRT(BEM(I))
DO 10 J=1,IC
A(I,J)=0.5D0*SQRT(BEC(J))/SIN(PI*TN)*
* COS(TN*(PI-ABS(PHIC(J)-PHIM(I))))
10 CONTINUE
C SET ATRN-MATRIX.
DO 20 I=1,IC
DO 20 J=1,IM
ATRN(I,J)=A(J,I)
20 CONTINUE
C WRITE(6,*)'A=',((A(I,J),I=1,IM),J=1,IC)
C WRITE(6,*)'AT=',((ATRN(I,J),I=1,IC),J=1,IM)
C COMPUTE AA-MATRIX.
CALL MATMUL(ATRN,A,AA,IC,IM,IC)
CALL MATSUB(CC,A,L,IC,LP)
C SET CCTRN-MATRIX. CCCCCCCCCC OK!!!
DO 21 I=1,IC
DO 21 J=1,L
CCTRN(I,J)=CC(J,I)
21 CONTINUE
Write(6,*) "Restraint Point
Number=",L
Write(6,*) "Restraint Point
=:"
Do 22 J = 1, L
Write(6,*)
" Position Monitor No.",LP(J)
22 CONTINUE
CALL VECZL(ETA,L,ZL,LP)
C WRITE(6,*)'AA=',((AA(I,J),I=1,IC),J=1,IC)
C COMPUTE EIGENVALUES AND EIGENVECTORS
OF AA-MATRIX.
CALL GIVHO(AA,E,V,IC,IC,IC)
C COMPUTE GENERALIZED INVERSE MATRIX OF AA. ( = W )
DO 30 I=1,IC
IF(E(I).LE.0.0.OR.I.GT.IE)THEN
EINV(I)=0.0
GO TO 30
ENDIF
EINV(I)=1.0D00/E(I)
30 CONTINUE
DO 40 I=1,IC
DO 40 J=1,IC
W(I,J)=0.0
DO 45 K=1,IC
W(I,J)=W(I,J)+EINV(K)*V(I,K)*V(J,K)
45 CONTINUE
40 CONTINUE
C COMPUTE PP-MATRIX.
C CALL MATRIXINV(W,WAA,IC)
CALL MATMUL(CC,W,PP,L,IC,IC)
C COMPUTE P-MATRIX.
CALL MATMUL(PP,CCTRN,P,L,IC,L)
CCCCCCCCCCCCCCCCCCCCCCCCCCCC WP=P OK!
C DO 60 I=1,L
C DO
50 J=1,L
C
C WP(I,J)=P(I,J)
C 50 CONTINUE
C 60 CONTINUE
CALL MATRIXINV(P,WP,L)
CALL MATMUL(W,CCTRN,QQ,IC,IC,L)
CALL MATMUL(QQ,WP,Q,IC,L,L)
CALL MATMUL(Q,CC,QCC,IC,L,IC)
CALL MATMUL(QCC,W,QCCW,IC,IC,IC)
CALL MATMUL(QCCW,ATRN,DM2,IC,IC,IM)
C This loop revived!
DO 50 I=1,IC
DO 50 J=1,IC
W(I,J)=-W(I,J)
50 CONTINUE
CALL MATMUL(W,ATRN,DM1,IC,IC,IM)
CALL MATADD(DM1,DM2,DM,IC,IM)
C COMPUTE KICK VECTOR(D).
CALL MATMUL(DM,ETA,DV1,IC,IM,1)
C This loop revived!
DO 60 I=1,IC
DO 60 J=1,L
Q(I,J)=-Q(I,J)
60 CONTINUE
CALL MATMUL(Q,ZL,DV2,IC,L,1)
CALL MATADD(DV1,DV2,D,IC,1)
CCCCCCCCCCC checking CCCCCCCCC
CALL MATMUL(CC,D,OK1,L,IC,1)
CALL MATADD(OK1,ZL,OK2,L,1)
RETURN
END
C *************************************************************
C *************************************************************
SUBROUTINE GIVHO(A,E,V,N,NEV,NVEC)
C *************************************************************
C *************************************************************
C PURPOSE
C COMPUTE EIGENVALUES
AND EIGENVECTORS OF A REAL NxN SYMMETRIC
C MATRIX BY THE GIVENS-HOUSEHOLDER
METHOD
C DESCRIPTION OF PARAMETERS
C A
- NxN INPUT MATRIX
C NEV - NUMBER
OF LARGEST EIGEN VALUES COMPUTED AND STORED IN E.
C NVEC - NUMBER
OF EIGEN VECTORS COMPUTED AND STORED IN V. NVEC<=NEV
C OVINV - INVERSE OF
OVERFLOW LIMIT FOR VAX COMPUTER
C IY
- INITIAL VALUE GENERATING RANDOM VECTOR BY FUNC. RAND
C COMMENT
C DIMENSION STATEMENT
VALID FOR NORDER UP TO 100
C THIS PROGRAM IS BASICALLY
VALID FOR NON-ZERO B(J).
C BUT EIGENVALUES AND
EIGENVECTORS MAY BE OBTAINED APPROXIMATELY,
C IF ZERO IS REPLACED
WITH A SMALL VALUE(1.0E-5 IN THIS CASE)
C IF YOU HAVE A TROUBLE,
YOU SHOULD USE AN ALTERNATE WAY
C (SEE J.M. ORTEGA, "ON
STURM SEQUENCES FOR TRIDIAGONAL MATRICES,"
C J. ASSOC. COMP.
MACH., 7, 260-263(1960).)
C
C
CODED BY N. NAKAMURA(25/10/1993)
IMPLICIT REAL*8 (A-H,O-Z)
REAL*8 A(150,150),E(150),V(150,150)
REAL*8 B(150),C(150),P(150),Q(150),R(150),W(150),Y(150)
DIMENSION IN(150)
INTEGER AG,IY, N, NEW, NVEC
REAL*8 KAP,NORM,LAMBDA,L,MULT
LOGICAL FIRST,IN
CNN
DATA EPS,EPS2,OVINV,IY/1.D-6,1.D-2,1.D-30,31415967D0/
CNN
NM1=N-1
NM2=N-2
C BEGIN THE REDUCTION TO TRIDIAGONAL FORM.
C THE DIAGONAL AND OFF-DIAGONAL ELEMENTS
OF THE TRIDIAGONAL MATRIX ARE
C STORED IN C AND B.
IF(N.LE.2)GO TO 9
DO 8 I=1,NM2
IP1=I+1
SS=0.
DO 1 J=IP1,N
SS=SS+A(J,I)*A(J,I)
C
Write(6,*) N, SS
1
CONTINUE
S=SQRT(SS)
IF(A(IP1,I).LT.0.)S=-S
C(I)=A(I,I)
B(I)=-S
C IF S=0 THEN ALPHA MUST BE ZERO.
ALPHA=0.
IF(S.EQ.0.0)GO TO 8
ALPHA=1./(SS+A(IP1,I)*S)
T=A(IP1,I)+S
A(IP1,I)=T
W(I+1)=T
IP2=I+2
DO 2 J=IP2,N
2
W(J)=A(J,I)
DO 4 J=IP1,N
T=0.0
DO 3 K=IP1,N
3
T=T+A(J,K)*W(K)
4
P(J)=T*ALPHA
KAP=0
DO 5 K=IP1,N
5
KAP=KAP+W(K)*P(K)
KAP=.5*KAP*ALPHA
DO 6 K=IP1,N
6
Q(K)=P(K)-KAP*W(K)
DO 7 J=IP1,N
C COMPUTE THE NEW A. BY SYMMETRY ONLY HALF
THE ELEMENTS NEED BE DONE.
DO 7 K=J,N
A(J,K)=A(J,K)-(Q(J)*W(K)+Q(K)*W(J))
7
A(K,J)=A(J,K)
8 A(I,I)=ALPHA
9 C(N-1)=A(N-1,N-1)
C(N)=A(N,N)
B(N-1)=A(N-1,N)
C WRITE(6,*)'C=',(C(I),I=1,N)
C WRITE(6,*)'B=',(B(I),I=1,N-1)
C WRITE(6,*)'A=',(A(I,I),I=1,N)
DO 100 I=1,N-1
IF(B(I).EQ.0.0)THEN
B(I)=1.E-5
C
WRITE(6,*)'B=0'
ENDIF
100 CONTINUE
C THIS COMPLETES THE REDUCTION TO TRIDIAGONAL
FORM.
C BEGIN THE EIGENVALUE COMPUTATION.
C COMPUTE NORM OF TRIDIAGONAL MATRIX, SQUARES
OF THE B(I).
NORM=ABS(C(1))+ABS(B(1))
C write(6,*) NORM
DO 10 I=2,N
T=ABS(C(I))+ABS(B(I))+ABS(B(I-1))
10 NORM = MAX(NORM,T)
C 10 NORM=AMAX1(NORM,T)
C WRITE(6,*)'NORM=',NORM
DO 11 I=1,N
11 W(I)=B(I)**2
C SET K AND UPPER AND LOWER ESTIMATES.
K=1
U=NORM
DO 12 I=1,NEV
E(I)=-NORM
12 CONTINUE
C write(6,*) NORM
C BEGIN MAIN LOOP.
13 L=E(K)
14 LAMBDA=.5*(L+U)
C write(6,*) LAMBDA
C THE CONVERGENCE TEST AS IMPLEMENTED HERE
ALLOWS THE COMPUTATION TO
C PROCEED UNTIL THE INTERVAL(L,U) CAN BE
MADE NO SMALLER.
IF((LAMBDA.EQ.L).OR.(LAMBDA.EQ.U))GO
TO 30
C BEGIN COMPUTATION OF NUMBER OF SIGN AGREEMENTS,AG
AG=0
I=1
16 S=C(I)-LAMBDA
18 IF(S.GE.0.0)AG=AG+1
IF(S.EQ.0.0)GO TO 20
I=I+1
IF(I.GT.N)GO TO 22
CNN
SNUM=S*(C(I)-LAMBDA)-W(I-1)
IF(SNUM.EQ.0.)THEN
S=0.0
GO TO 18
ENDIF
SINV=S/SNUM
IF(ABS(SINV).GE.OVINV)THEN
S=C(I)-LAMBDA-W(I-1)/S
GO TO 18
ENDIF
CNN
C IF OVERFLOW OCURRS, PROCEED AS IF S HAD
BEEN ZERO.
AG=AG+1
I=I-1
20 I=I+2
IF(I.LE.N)GO TO 16
C THE COMPUTATION OF AG IS COMPLETE. ADJUST
INTERVAL.
22 IF(AG.GE.K)GO TO 24
U=LAMBDA
GO TO 14
24 L=LAMBDA
M=MIN0(AG,NEV)
DO 26 I=K,M
26 E(I)=LAMBDA
GO TO 14
C THE K-TH EIGENVALUE IS COMPUTED. STORE
IN E(K) AND PROCEED.
30 E(K)=LAMBDA
K=K+1
IF(K.LE.NEV)GO TO 13
C THIS COMPLETES THE EIGENVALUE COMPUTATION.
C IF NVEC NOT ZERO, BEGIN THE EIGENVECTOR
COMPUTATION.
IF(NVEC.NE.0)GO TO 40
RETURN
40 DO 90 I=1,NVEC
C INITIALIZATION FOR THIS EIGENVECTOR.
DO 44 J=1,N
P(J)=0.
Q(J)=B(J)
R(J)=C(J)-E(I)
C WRITE(6,*)'J,P,Q,R=',J,P(J),Q(J),R(J)
44
Y(J)=1.
Y(N+1)=0.
Y(N+2)=0.
FIRST=.TRUE.
C BEGIN THE GAUSSIAN ELIMINATION REDUCTION
TO TRIANGULAR FORM.
DO 50 J=1,NM1
IF(ABS(R(J)).LT.ABS(B(J)))GO TO 46
C WRITE(6,*)'J,R(J)=',J,R(J)
C IF(R(J).EQ.0.)R(J)=OVINV
MULT=B(J)/R(J)
IN(J)=.FALSE.
GO TO 48
46
MULT=R(J)/B(J)
IN(J)=.TRUE.
R(J)=B(J)
T=R(J+1)
R(J+1)=Q(J)
Q(J)=T
P(J)=Q(J+1)
Q(J+1)=0.
C STORE THE MULTIPLIERS FOR LATTER USE.
48
W(J)=MULT
Q(J+1)=Q(J+1)-MULT*P(J)
R(J+1)=R(J+1)-MULT*Q(J)
IF(ABS(R(J)).EQ.0.0)R(J)=OVINV
50 CONTINUE
IF(ABS(R(N)).EQ.0.0)R(N)=OVINV
C DO 110 J=1,N
C 110 WRITE(6,*)'J,P,Q,R=',J,P(J),Q(J),R(J)
C THE COMPLETES THE REDUCTION. BEGIN THE
BACK SUBSTITUTION.
C IF E(I) IS A CLUSTER, GENERATE A RANDOM
VECTOR.
C WRITE(6,*)'NORM=',NORM
C IF(I.NE.1.AND.ABS(E(I)-E(I-1)).LT.NORM*EPS)E(I)=E(I-1)-NORM*EPS*2
IF(I.EQ.1.OR.ABS(E(I)-E(I-1)).GE.NORM*EPS)GO TO 54
51 CONTINUE
C WRITE(6,*)'RAN,I=',I
DO 52 J=1,N
52 Y(J)=RAN(IY)
C WRITE(6,*)'Y(RAN)=',(Y(J),J=1,N)
54 CONTINUE
DO 66 JI=1,N
K=N-JI+1
C SAVE Y(K) TO SCALE IF OVERFLOW.
T=Y(K)
C Y(N+1) AND Y(N+2) ARE ZERO.
62
CONTINUE
YNUM=T-Y(K+1)*Q(K)-Y(K+2)*P(K)
C WRITE(6,*)'YNUM=',YNUM
IF(ABS(YNUM).LE.OVINV)THEN
Y(K)=0.0
GO TO 66
ENDIF
YINV=R(K)/YNUM
C WRITE(6,*)'K,R,YINV=',K,R(K),YINV
IF(ABS(YINV).GE.OVINV)THEN
Y(K)=(T-Y(K+1)*Q(K)-Y(K+2)*P(K))/R(K)
GO TO 66
ENDIF
C WRITE(6,*)'Y62(OVFL)=',(Y(J),J=1,N)
DO 64 J=1,N
64 Y(J)=Y(J)*1.E-5
T=T*1.E-5
C WRITE(6,*)'Y62*1E-5=',(Y(J),J=1,N)
GO TO 62
66 CONTINUE
C THIS COMPLETES BACK SUBSTITUTION.
IF(.NOT.FIRST)GO TO 74
FIRST=.FALSE.
C DO THE REDUCTION ON THE NEW RIGHT HAND
SIDE.
DO 70 J=1,NM1
IF(IN(J))GO TO 68
Y(J+1)=Y(J+1)-W(J)*Y(J)
GO TO 70
68
T=Y(J)
Y(J)=Y(J+1)
Y(J+1)=T-W(J)*Y(J+1)
70 CONTINUE
C THIS COMPLETES REDUCTION OF RIGHT SIDE. DO ONE MORE ITERATION.
GO TO 54
C Y IS NOW AN EIGENVECTOR OF TRIDIAGONAL
MATRIX. BEGIN BACK
C TRANSFORMATION.
74 DO 78 J=1,NM2
T=0.0
K=N-J-1
M=K+1
DO 76 KK=M,N
T=T+A(KK,K)*Y(KK)
76
Continue
T=A(K,K)*T
DO 78 KK=M,N
Y(KK)=Y(KK)-T*A(KK,K)
78 Continue
C BACK TRANSFORMATION COMPLETE. Y NOW EIGENVECTOR OF A. NORMALIZE.
T=ABS(Y(1))
K=1
DO 80 J=2,N
S=ABS(Y(J))
IF(S.LE.T)GO TO 80
T=S
K=J
80 CONTINUE
T=1.0/Y(K)
DO 82 J=1,N
V(J,I)=Y(J)*T
82
CONTINUE
C WRITE(6,*)'Y=',(Y(J),J=1,N)
C I-TH COLUMN OF V IS I-TH EIGENVECTOR
NORMALIZED WITH MAX ELEMENT =1
C NORMALIZE EIGENVECTORS WITH NORM=1
TV=0.0
DO 84 J=1,N
TV=TV+V(J,I)*V(J,I)
84 CONTINUE
TV=1./SQRT(TV)
DO 86 J=1,N
V(J,I)=V(J,I)*TV
86 CONTINUE
IF(I.EQ.1)GO TO 90
IF(ABS(E(I)-E(I-1)).LT.NORM*EPS)THEN
PROD=0.0
DO 88 J=1,N
C WRITE(6,*)'V(I-1),V(I)=',V(J,I-1),V(J,I)
PROD=PROD+V(J,I-1)*V(J,I)
88 CONTINUE
C WRITE(6,*)'PROD=',PROD
IF(ABS(PROD).GE.EPS2)THEN
FIRST=.TRUE.
GO TO 51
ENDIF
ENDIF
C
90 CONTINUE
C
RETURN
END
C *************************************************************
C *************************************************************
SUBROUTINE MATRIXINV(A,B,M)
C *************************************************************
C *************************************************************
C PURPOSE
C INVERT A
( Gausse-Jordan Meathod )
C DISCRIPTION OF PARAMETERS
C A
- INPUT MATRIX
C AA
- DATA MATRIX ( M X 2M )
C
C
C
C COMMENTS
C DIMENSION STATEMENT
VALID FOR NORDER UP TO 100
C
C
CODED BY Masanori SATOH (29/5/1998)
C
IMPLICIT REAL*8 (A-H,O-Z)
INTEGER I,J,K,L,M,N
REAL*8 A(150,150),AA(150,300),B(150,150)
REAL*8 P
N=2*M
L=M+1
DO 14 I=1,M
DO
10 J=1,M
AA(I,J)=A(I,J)
10 CONTINUE
DO 12 J=L,N
IF (J.EQ.M+I) THEN
AA(I,J)=1.D0
ELSE
AA(I,J)=0.D0
END IF
12 CONTINUE
14 CONTINUE
DO 50 K=1,M
P=AA(K,K)
DO 20 J=K,N
AA(K,J)=AA(K,J)/P
20 CONTINUE
DO 40 I=1,M
IF(I-K)21,40,21
21 AIK=AA(I,K)
DO 30 J=K,N
AA(I,J)=AA(I,J)-AIK*AA(K,J)
30 CONTINUE
40 CONTINUE
50 CONTINUE
DO 70 I=1,M
DO 60 J=1,M
B(I,J)=AA(I,J+M)
60 CONTINUE
70 CONTINUE
C CLOSE(12)
130 CONTINUE
140 RETURN
END
C *************************************************************
C *************************************************************
C Matrix Calculation Series
C *************************************************************
C *************************************************************
SUBROUTINE MATADD(A,B,C,M,N)
C PURPOSE
C COMPUTE SUM OF
TWO MxN MATRICE.
C DESCRIPTION OF PARAMETER
C A,B - MxN MATRIX
C C
- C=A+B
C COMMENT
C DIMENSION STATEMENT
VALID FOR N AND M UP TO 100
C
CODED BY N. NAKAMURA (26/10/1993)
REAL*8 A(150,150),B(150,150),C(150,150)
DO 10 I=1,M
DO 10 J=1,N
C(I,J)=A(I,J)+B(I,J)
10 CONTINUE
RETURN
END
C
SUBROUTINE MATRED(A,B,C,M,N)
C PURPOSE
C COMPUTE DIFFERENCE
OF TWO MxN MATRICE.
C DESCRIPTION OF PARAMETER
C A,B - MxN MATRIX
C C
- C=A-B
C COMMENT
C DIMENSION STATEMENT
VALID FOR N AND M UP TO 100
C
CODED BY N. NAKAMURA (26/10/1993)
IMPLICIT REAL*8 (A-H,O-Z)
REAL*8 A(150,150),B(150,150),C(150,150)
DO 10 I=1,M
DO 10 J=1,N
C(I,J)=A(I,J)-B(I,J)
10 CONTINUE
RETURN
END
C
SUBROUTINE MATMUL(A,B,C,M,L,N)
C PURPOSE
C COMPUTE SUM OF
TWO MxN MATRICE.
C DESCRIPTION OF PARAMETER
C A - MxL MATRIX
C B - LxN MATRIX
C C - C=A*B,MxN
MATRIX
C COMMENT
C DIMENSION STATEMENT
VALID FOR N,L AND M UP TO 100
C
CODED BY N. NAKAMURA (26/10/1993)
IMPLICIT REAL*8 (A-H,O-Z)
REAL*8 A(150,150),B(150,150),C(150,150)
C REAL I,J,K
DO 10 I=1,M
DO 10 J=1,N
C(I,J)=0.D0
DO 10 K=1,L
C(I,J)=C(I,J)+A(I,K)*B(K,J)
10 CONTINUE
RETURN
END
C
SUBROUTINE MATSUB(CC,A,L,IC,LP)
C PURPOSE
C COMPUTE SUB-MATRIX
OF Response MATRICE.
C
C DESCRIPTION OF PARAMETER
C A - MxN Response
MATRIX
C CC - LXM
MATRIX
C
C COMMENT
C DIMENSION STATEMENT
VALID FOR N,L AND M UP TO 100
C
CODED BY Masanori SATOH (25/5/1998)
IMPLICIT REAL*8 (A-H,O-Z)
REAL*8 A(150,150),CC(150,150)
INTEGER LP(150)
DO 10 I=1,L
DO 10 J=1,IC
CC(I,J)=A(LP(I),J)
10 CONTINUE
C CC Matrix Data Writing !
C OPEN(UNIT=11,FILE='cc.dat',STATUS='OLD',FORM='FORMATTED')
C
C DO 11 I=1,L
C DO
11 J=1,IC
C
C write(11,*) CC(I,J),A(I,J)
C
C 11 CONTINUE
C CLOSE(11)
RETURN
END
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
SUBROUTINE VECZL(Z,L,ZL,LP)
IMPLICIT REAL*8 (A-H,O-Z)
REAL*8 Z(150),ZL(150)
INTEGER LP(150)
DO 10 I=1,L
ZL(I)=Z(LP(I))
10 CONTINUE
RETURN
END