CHARACTER*132 JUNK CHARACTER*15 NAME,GEO real Z(405,600),R(405,600),X(405,600),Y(405,600) real VX(405,600),VY(405,600),PVX(405),PVY(405) real ZG(10,50),RG(10,50),PX(405),PY(405),BZ(405) REAL RAQ(405),RAQI(405),QQI(3,405),RQQI(3,405) REAL OA(405),OQ(405),VZ(405,600),EZ(10),ZB(405) REAL VVX(405),VVY(405),XXZ(405),YYZ(405),ZZE(10) REAL XV(405),YV(405),ZX(405),ZY(405),ALX(2),ALY(2) real tt(405),ttem(10),AP(405),API(405),APO(405) REAL OAP(405),XE(405),YE(405),TNOTX(10),TNOTY(10) real Q(405),AA(405),QI(405),AO(405),QO(405) CHARACTER*4 NUMB INTEGER ZAX(3),RAX(3),IZE(10),BZ0(3),NUB(405),LG(10) DATA BZ0/'B(Z)',' GAU','SS'/ DATA ZAX/'Z ','( mm',' )'/,RAX/'r ','( mm',' )'/ INTEGER XAX(3),YAX(3),NS(405) DATA XAX/'X ','( mm',' )'/,YAX/'Y ','( mm',' )'/ INTEGER VAX(3),VAY(3) DATA VAX(1)/'X ('/,VAX(2)/'mrad'/,VAX(3)/') '/ DATA VAY(1)/'Y ('/,VAY(2)/'mrad'/,VAY(3)/') '/ INTEGER BI(3),AQR(3) DATA BI/'I (',' euA',' )'/,AQR/' A/Q',' ',' '/ CTAA TYPE *,' NAME' CTAA READ(5,FMT='(A15)')NAME TYPE *,' GEOMETRY FILE NAME' READ(5,FMT='(A15)')GEO TYPE *, ' HOW MANY RAYS ? TYPE M ' READ(5,*) M TYPE *,' SOLENOID LOCATION (M) AND SLIT WIDTH ?' READ(5,*) SOL_LOC,SLIT_WIDTH OPEN(100,FILE=GEO,STATUS='OLD') LGG=0 DO I=1,10 READ(100,*,END=110)LG(I),(ZG(I,J),J=1,LG(I)), &(RG(I,J),J=1,LG(I)) LGG=LGG+1 ENDDO 110 DO I=1,4 ZG(4,I)=ZG(4,I)+SOL_LOC*1000. ENDDO DO I=2,3 RG(5,I)=RG(5,I)+SLIT_WIDTH*1000. ENDDO TYPE *,' EMITTANCE FITTING FRAC , TYPE NO OF LOCATIONS AND 1 THE Z_STEP ( MM ) IN THE CALCULATION?' READ(5,*) FRAC,IZ,ZSTEP TYPE *,' Z POSITIONS FOR THE EMITTANCE PLOTS?' READ(5,*) (EZ(I),I=1,IZ) ctaa ctaa emittance locations are onput in meters- conver to mm ctaa DO I=1,IZ EZ(I) = EZ(I)*1000. ENDDO TYPE *,' ray number of the calculated main species?' READ(5,*) mm,Mr OPEN(1,STATUS='OLD') OPEN(11,STATUS='NEW',form='unformatted') RMAX=0. RXX=0. DO K=1,2 READ(1,FMT='(A132)')JUNK ENDDO DO I=1,M READ(1,*) N,AA(I),Q(I),AP(I) RAQ(I)=AA(I)/Q(I) ENDDO DO I=1,M QO(I)=Q(I) AO(I)=AA(I) APO(I)=AP(I) ENDDO DO I=2,M DO J=I,M IF(RAQ(J).LT.RAQ(I-1)) THEN RAT=RAQ(I-1) QT=Q(I-1) APT=AP(I-1) RAQ(I-1)=RAQ(J) Q(I-1)=Q(J) AP(I-1)=AP(J) AP(J)=APT RAQ(J)=RAT Q(J)=QT ENDIF ENDDO ENDDO NT=0 LI=0. DO 301 I=1,M DO LN=1,NT IF(I.EQ.NS(LN)) GOTO 301 ENDDO LI=LI+1 API(LI)=0. RAQI(LI)=RAQ(I) DO 401 J=I,M DO LN=1,NT IF(J.EQ.NS(LN)) GOTO 401 ENDDO IF(RAQ(J).EQ.RAQ(I)) THEN API(LI)=API(LI)+AP(J) NT=NT+1 NS(NT)=J ENDIF 401 CONTINUE 301 CONTINUE APMX=0. DO LL=1,LI IF(API(LL).GT.APMX) APMX=API(LL) ENDDO DO K=1,12 READ(1,FMT='(A132)')JUNK ENDDO RZX=-3. KEZ=0 tte=0. DO 2 K=1,1000 do 1 I=1,M READ(1,*,END=999)N,Z(I,K),X(I,K),Y(I,K),VX(I,K),VY(I,K),VZ(I,K), 1 tt(k),tl,BZ(K) R(I,K)=(X(I,K)*X(I,K)+Y(I,K)*Y(I,K))**.5 AX=ABS(X(I,K)) AY=ABS(Y(I,K)) AZ=Z(I,K) if(tte.lt.tt(k)) tte=tt(k) IF(AX.GT.RXX) RXX=AX IF(AY.GT.RXX) RXX=AY IF(AZ.GT.RZX) RZX=AZ 1 IF(R(I,K).GT.RMAX) RMAX=R(I,K) BZ(K)=ABS(BZ(K))*10000. ZB(K)=RZX DO J=1,IZ IF((ABS(RZX-EZ(J)).LT.ZSTEP/2.).OR.(RZX-EZ(J). 1 EQ.ZSTEP/2.).OR.(RZX.EQ.EZ(J))) THEN KEZ=KEZ+1 ttem(kez)=tte IZE(KEZ)=K ZZE(KEZ)=RZX ENDIF ENDDO 2 continue 999 L=K-1 l1=l NZ=RZX/25.4+1 NR=RMAX/25.4+1 I=6 l=1 a=float(NZ) if(a.lt.8.) a=8. b=float(NR) IF(b.lt.6.) b=6. write(11)I,l,a,b,a,b,d,d i=1 l=3 b=25.4 c=0. d=-4. write(11)i,l,a,b,c,d,ZAX,ZAX a=float(NR) IF(a.lt.6.) a=6. i=2 b=25.4 d=0. write(11)i,l,a,b,c,d,RAX,RAX b=500. c=float(NZ) write(11)i,l,a,b,c,d,BZ0,BZ0 I=7 l=K-1 B=25.4 A=25.4 C=-4. D=0. do 3 n1=1,M 3 WRITE(11)i,l,a,b,c,d,(Z(n1,j),j=1,l),(R(n1,j),j=1,l) IF(RZX.LT.SOL_LOC) LGG=3 DO J=1,LGG l=LG(J) write(11)i,l,a,b,c,d,(ZG(J,k),k=1,l),(RG(J,k),k=1,l) ENDDO l=l1 b=500. write(11)i,l,a,b,c,d,(ZB(k),k=1,l),(BZ(k),k=1,l) C i=3 l=1 a=0. b=0. write(11)i,l,a,b,c,d,c,c i=6 l=1 a=8. b=8. write(11)i,l,a,b,A,B,d,d i=1 l=3 a=8. b=1.25 c=0. d=-4.*b write(11)i,l,a,b,c,d,XAX,XAX i=2 write(11)i,l,a,b,c,d,YAX,YAX i=5 l=M a=1.25 c=d DO J=1,L ZX(J)=X(J,1)/A+4. ZY(J)=Y(J,1)/B+4. NUB(J)=J ENDDO WRITE(11)i,l,a,b,c,d,(ZX(J),J=1,L),(ZY(J),J=1,L) WRITE(11)(NUB(INB),INB=1,L) C i=3 l=1 a=0. b=0. write(11)i,l,a,b,c,d,c,c i=6 l=1 a=8. b=8. write(11)i,l,a,b,A,B,d,d MQ=APMX/8.+1. MA=RAQ(M)/8.+1 i=1 l=3 a=8. b=FLOAT(MA) c=0. d=0. write(11)i,l,a,b,c,d,AQR,AQR i=2 b=FLOAT(MQ) write(11)i,l,a,b,c,d,BI,BI i=0 l=2 a=FLOAT(MA) DO NJ=1,LI RQQI(1,NJ)=RAQI(NJ) RQQI(2,NJ)=RAQI(NJ) QQI(1,NJ)=0. QQI(2,NJ)=API(NJ) ENDDO DO LJ=1,LI WRITE(11)i,l,a,b,c,d,(RQQI(j,LJ),j=1,l),(QQI(j,LJ),j=1,l) ENDDO C C 400 i=3 l=1 a=0. b=0. write(11)i,l,a,b,c,d,c,c NXY=RXX NDXY=FLOAT(NXY)/4.+2 ND=NDXY i=6 l=1 a=8. b=8. write(11)i,l,a,b,A,B,d,d i=1 l=3 a=8. b=FLOAT(ND) c=0. d=-4.*b write(11)i,l,a,b,c,d,XAX,XAX i=2 write(11)i,l,a,b,c,d,YAX,YAX i=7 l=l1 a=FLOAT(ND) c=d do 5 n1=1,M 5 WRITE(11)i,l,a,b,c,d,(X(n1,j),j=1,l),(Y(n1,j),j=1,l) C DO 1000 J=1,KEZ VMX=0. VMY=0. XMX=0. YMX=0. IA=0 II=IZE(J) DO 800 IV=mm,Mr IF(Z(IV,II).LT.ZZE(J)) GOTO 800 IA=IA+1 NUB(IA)=IV VVX(IA)=VX(IV,II)/VZ(IV,II)*1000. VVY(IA)=VY(IV,II)/VZ(IV,II)*1000. XXZ(IA)=X(IV,II) YYZ(IA)=Y(IV,II) IF(ABS(VVX(IA)).GT.VMX) VMX=ABS(VVX(IA)) IF(ABS(VVY(IA)).GT.VMY) VMY=ABS(VVY(IA)) IF(ABS(XXZ(IA)).GT.XMX) XMX=ABS(XXZ(IA)) IF(ABS(YYZ(IA)).GT.YMX) YMX=ABS(YYZ(IA)) 800 CONTINUE DVX=VMX*3./8. DX=XMX*3./8. IF(DVX.EQ.0.) GOTO 300 i=3 l=1 a=0. b=0. write(11)i,l,a,b,c,d,c,c i=6 l=1 a=8. b=8. write(11)i,l,a,b,A,B,d,d i=1 l=3 a=8. b=20. c=0. d=-4.*b write(11)i,l,a,b,c,d,XAX,XAX i=2 b=20. d=-4.*b l=3 write(11)i,l,a,b,c,d,VAX,VAX i=5 l=IA a=20. c=-4.*a DO J1=1,L ZX(J1)=XXZ(J1)/A+4. XV(J1)=VVX(J1)/B+4. ENDDO WRITE(11)i,l,a,b,c,d,(ZX(j1),j1=1,l),(XV(j1),j1=1,l) WRITE(11)(NUB(INB),INB=1,L) CALL ELLFIT(XXZ,VVX,IA,FRAC,XC,YC,AE,BE,ANG) CALL ELLIPS(AE,BE,XC,YC,ANG,IA,XE,YE,NUSE) i=7 l=IA+1 XE(L)=XE(1) YE(L)=YE(1) Ctaa OPEN(20,FILE='EFIT.DAT',STATUS='NEW') ctaa do ik=1,l ctaa write(20,*)XE(IK),YE(IK) ctaa ENDDO ctaa CLOSE(20) WRITE(11)i,l,a,b,c,d,(XE(j1),j1=1,l),(YE(j1),j1=1,l) ALX(1)=0. ALY(1)=-4.*20. ALX(2)=0. ALY(2)=4.*20. i=7 l=2 a=20. c=-4.*a WRITE(11)i,l,a,b,c,d,(ALX(j1),j1=1,l),(ALY(j1),j1=1,l) ALX(1)=-4.*20. ALY(1)=0. ALX(2)=4.*20. ALY(2)=0. i=7 l=2 a=20. c=-4.*a WRITE(11)i,l,a,b,c,d,(ALX(j1),j1=1,l),(ALY(j1),j1=1,l) CALL ELLCON(ANG,AE,BE,AREA,R12,XINT,YINT,XMAX,YMAX,YATXM,XATYM) TNOTX(1)=ZZE(J) TNOTX(2)=XINT+XC TNOTX(3)=XMAX+XC TNOTY(1)=YINT+YC TNOTY(2)=YMAX+YC TNOTY(3)=AREA TNOTY(4)=YC TNOTX(4)=XC TNOTY(5)=ANG TNOTX(5)=ANG i=9 l=5 WRITE(11)i,l,a,b,c,d,(TNOTX(j1),j1=1,l),(TNOTY(j1),j1=1,l) i=3 l=1 a=0. b=0. write(11)i,l,a,b,c,d,c,c i=6 l=1 a=8. b=8. write(11)i,l,a,b,A,B,d,d i=1 l=3 a=8. b=20. c=0. d=-4.*b write(11)i,l,a,b,c,d,XAX,XAX i=2 b=20. d=-4.*b l=3 write(11)i,l,a,b,c,d,VAX,VAX i=10 l=IA a=20. c=-4.*a WRITE(11)I,L,A,B,C,D,(XXZ(J1),J1=1,L),(VVX(J1),J1=1,L) i=7 l=IA+1 WRITE(11)i,l,a,b,c,d,(XE(j1),j1=1,l),(YE(j1),j1=1,l) ALX(1)=0. ALY(1)=-4.*20. ALX(2)=0. ALY(2)=4.*20. i=7 l=2 a=20. c=-4.*a WRITE(11)i,l,a,b,c,d,(ALX(j1),j1=1,l),(ALY(j1),j1=1,l) ALX(1)=-4.*20. ALY(1)=0. ALX(2)=4.*20. ALY(2)=0. i=7 l=2 a=20. c=-4.*a WRITE(11)i,l,a,b,c,d,(ALX(j1),j1=1,l),(ALY(j1),j1=1,l) i=9 l=5 WRITE(11)i,l,a,b,c,d,(TNOTX(j1),j1=1,l),(TNOTY(j1),j1=1,l) C 300 DVY=VMY*3./8. DY=YMX*3./8. IF(DVY.EQ.0.) GOTO 1000 NDVY=DVY+1. NDY=DY+1. i=3 l=1 a=0. b=0. write(11)i,l,a,b,c,d,c,c i=6 l=1 a=8. b=8. write(11)i,l,a,b,A,B,d,d i=1 l=3 a=8. b=20. c=0. d=-4.*b write(11)i,l,a,b,c,d,YAX,YAX i=2 b=20. d=-4.*b l=3 write(11)i,l,a,b,c,d,VAY,VAY i=5 l=IA a=20. c=-4.*a DO J1=1,L ZY(J1)=YYZ(J1)/A+4. YV(J1)=VVY(J1)/B+4. ENDDO WRITE(11)i,l,a,b,c,d,(ZY(j1),j1=1,l),(YV(j1),j1=1,l) WRITE(11)(NUB(INB),INB=1,L) CALL ELLFIT(YYZ,VVY,IA,FRAC,XC,YC,AE,BE,ANG) CALL ELLIPS(AE,BE,XC,YC,ANG,IA,XE,YE,NUSE) i=7 l=IA+1 XE(L)=XE(1) YE(L)=YE(1) WRITE(11)i,l,a,b,c,d,(XE(j1),j1=1,l),(YE(j1),j1=1,l) ALX(1)=0. ALY(1)=-4.*20. ALX(2)=0. ALY(2)=4.*20. i=7 l=2 a=20. c=-4.*a WRITE(11)i,l,a,b,c,d,(ALX(j1),j1=1,l),(ALY(j1),j1=1,l) ALX(1)=-4.*20. ALY(1)=0. ALX(2)=4.*20. ALY(2)=0. i=7 l=2 a=20. c=-4.*a WRITE(11)i,l,a,b,c,d,(ALX(j1),j1=1,l),(ALY(j1),j1=1,l) CALL ELLCON(ANG,AE,BE,AREA,R12,XINT,YINT,XMAX,YMAX,YATXM,XATYM) TNOTX(1)=ZZE(J) TNOTX(2)=XINT+XC TNOTX(3)=XMAX+XC TNOTY(1)=YINT+YC TNOTY(2)=YMAX+YC TNOTY(3)=AREA TNOTX(4)=XC TNOTY(4)=YC TNOTY(5)=ANG TNOTX(5)=ANG i=9 l=5 WRITE(11)i,l,a,b,c,d,(TNOTX(j1),j1=1,l),(TNOTY(j1),j1=1,l) i=3 l=1 a=0. b=0. write(11)i,l,a,b,c,d,c,c i=6 l=1 a=8. b=8. write(11)i,l,a,b,A,B,d,d i=1 l=3 a=8. b=20. c=0. d=-4.*b write(11)i,l,a,b,c,d,YAX,YAX i=2 b=20. d=-4.*b l=3 write(11)i,l,a,b,c,d,VAY,VAY i=10 l=IA a=20. c=-4.*a WRITE(11)I,L,A,B,C,D,(YYZ(J1),J1=1,L),(VVY(J1),J1=1,L) i=7 l=IA+1 WRITE(11)i,l,a,b,c,d,(XE(j1),j1=1,l),(YE(j1),j1=1,l) ALX(1)=0. ALY(1)=-4.*20. ALX(2)=0. ALY(2)=4.*20. i=7 l=2 a=20. c=-4.*a WRITE(11)i,l,a,b,c,d,(ALX(j1),j1=1,l),(ALY(j1),j1=1,l) ALX(1)=-4.*20. ALY(1)=0. ALX(2)=4.*20. ALY(2)=0. i=7 l=2 a=20. c=-4.*a WRITE(11)i,l,a,b,c,d,(ALX(j1),j1=1,l),(ALY(j1),j1=1,l) i=9 l=5 WRITE(11)i,l,a,b,c,d,(TNOTX(j1),j1=1,l),(TNOTY(j1),j1=1,l) 1000 CONTINUE C C i=3 l=1 a=0. b=0. write(11)i,l,a,b,c,d,c,c i=6 l=1 a=8. b=8. write(11)i,l,a,b,A,B,d,d i=1 l=3 a=8. b=FLOAT(ND) c=0. d=-4.*b write(11)i,l,a,b,c,d,XAX,XAX i=2 write(11)i,l,a,b,c,d,YAX,YAX KL=0 do 10 LI1=1,M IF(Z(LI1,L1).LT.RZX) GOTO 10 KL=KL+1 NUB(KL)=LI1 PX(KL)=X(LI1,L1) PY(KL)=Y(LI1,L1) OQ(KL)=QO(LI1) OA(KL)=AO(LI1) OAP(KL)=APO(LI1) RAQ(KL)=OA(KL)/OQ(KL) 10 CONTINUE i=5 l=KL a=FLOAT(ND) c=d DO J=1,L PX(J)=PX(J)/A+4. PY(J)=PY(J)/B+4. ENDDO WRITE(11)i,l,a,b,c,d,(PX(J),J=1,L),(PY(J),J=1,L) WRITE(11)(NUB(INB),INB=1,L) C DO I=2,KL DO J=I,KL IF(RAQ(J).LT.RAQ(I-1)) THEN RAT=RAQ(I-1) QT=OQ(I-1) APT=OAP(I-1) RAQ(I-1)=RAQ(J) OQ(I-1)=OQ(J) OAP(I-1)=OAP(J) OAP(J)=APT RAQ(J)=RAT OQ(J)=QT ENDIF ENDDO ENDDO NT=0 LI=0. DO 601 I=1,KL DO LN=1,NT IF(I.EQ.NS(LN)) GOTO 601 ENDDO LI=LI+1 API(LI)=0. RAQI(LI)=RAQ(I) DO 501 J=I,KL DO LN=1,NT IF(J.EQ.NS(LN)) GOTO 501 ENDDO IF(RAQ(J).EQ.RAQ(I)) THEN API(LI)=API(LI)+OAP(J) NT=NT+1 NS(NT)=J ENDIF 501 CONTINUE 601 CONTINUE i=3 l=1 a=0. b=0. write(11)i,l,a,b,c,d,c,c i=6 l=1 a=8. b=8. write(11)i,l,a,b,A,B,d,d i=1 l=3 a=8. b=FLOAT(MA) c=0. d=0. write(11)i,l,a,b,c,d,AQR,AQR i=2 b=FLOAT(MQ) write(11)i,l,a,b,c,d,BI,BI i=0 l=2 a=FLOAT(MA) DO NJ=1,LI RQQI(1,NJ)=RAQI(NJ) RQQI(2,NJ)=RAQI(NJ) QQI(1,NJ)=0. QQI(2,NJ)=API(NJ) ENDDO DO LJ=1,LI WRITE(11)i,l,a,b,c,d,(RQQI(j,LJ),j=1,l),(QQI(j,LJ),j=1,l) ENDDO C 600 i=3 l=1 a=0. b=0. write(11)i,l,a,b,c,d,c,c i=4 write(11)i,l,a,b,c,d,c,c STOP END SUBROUTINE ELLIPS(A,B,XC,YC,ANG,N,X,Y,NUSE) C C LIBRARY-ROUTINE C C 05/AUG/1980 C C.J. KOST SIN C C====================================================================== C====================================================================== C== == C== ELLIPS: RETURNS NUSE POINTS IN (X, Y) WHICH UNIFORMLY POPULATE == C== THE PERIMETER OF THE ELLIPSE SPECIFIED BY THE MAJOR AND == C== MINOR AXES: A & B (RESPECTIVELY), THE X AND Y COORDINATES= C== OF THE CENTER OF THE ELLIPSE: XC AND YC, AND THE ANGLE == C== OF ORIENTATION (W.R.T. THE X-AXIS): ANG. N IS THE == C== NUMBER OF POINTS WHICH ELLIPS TRYS TO PRODUCE, WHILE == C== NUSE=N/4 * 4 IS THE ACTUAL NUMBER IT RETURNS IN (X,Y). == C== == C== INPUT PARAMETERS: A,B,XC,YC,ANG,N. == C== == C== OUTPUT PARAMETERS: X,Y,NUSE. == C== == C====================================================================== C====================================================================== C== == C== SUBROUTINES CALLED: == C== == C== NAME TYPE EXPLANATION == C== == C== COS R*4 COSINE FUNCTION. == C== SIN R*4 SINE FUNCTION. == C== == C====================================================================== C====================================================================== C== == C== VARIABLES: == C== == C== NAME TYPE EXPLANATION == C== == C== A R*4 SEMI-MAJOR AXIS LENGTH. == C== ANG R*4 ANGLE OF ORIENTATION W.R.T TO THE X-AXIS. == C== MEASURED IN DEGREES. == C== B R*4 SEMI-MINOR AXIS LENGTH. == C== CANGR R*4 COSINE OF ANG. == C== I I*4 DO LOOP INDEX. == C== N I*4 NUMBER OF POINTS ELLIPS TRYS TO POPULATE THE == C== ELLIPSE WITH. ACTUAL NUMBER OF POINTS RETURN-== C== ED IN X,Y IS NUSE=N/4*4. == C== ND2M1 I*4 N/2-1. == C== ND4 I*4 N/4. == C== ND4P1 I*4 N/4+1. == C== NUSE I*4 ACTUAL NUMBER OF POINTS USED TO POPULATE THE == C== ELLIPSE. NUSE=N/4*4. == C== PI R*4 WELL KNOWN CONSTANT: 3.14159265. == C== SANGR R*4 SINE OF ANGLE ANG. == C== THINC R*4 THETA INCREMENT. == C== X R*4(N) NUSE X-COORDINATES OF THE ELLIPSE. == C== XC R*4 X-COORDINATE OF THE CENTER OF THE ELLIPSE. == C== XX R*4 X-COORDINATE OF THE ELLIPSE. == C== Y R*4(N) NUSE Y-COORDINATES OF THE ELLIPSE. == C== YC R*4 Y-COORDINATE OF THE CENTER OF THE ELLIPSE. == C== YY R*4 Y-COORDINATE OF THE ELLIPSE. == C== == C====================================================================== C====================================================================== CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC REAL X(N),Y(N) PI=3.14159265 ND4=N/4 NUSE=4*ND4 IF(ND4.EQ.0)RETURN ND4P1=ND4+1 THINC=PI/2./ND4 C====================================================================== C== == C== POPULATE 1/4 OF THE ELLIPSE WITH N/4 POINTS IN QUADRANT 1 OF == C== THE STANDARD COORDINATE SYSTEM FOR THE ELLIPSE. == C== == C====================================================================== DO 10 I=1,ND4P1 X(I)=A*COS((I-1)*THINC) Y(I)=B*SIN((I-1)*THINC) 10 CONTINUE C====================================================================== C== == C== POPULATE 1/4 OF THE ELLIPSE WITH N/4 POINTS IN QUADRANT 2 OF == C== THE STANDARD COORDINATE SYSTEM FOR THE ELLIPSE. == C== == C====================================================================== DO 20 I=1,ND4 X(I+ND4P1)=-X(ND4P1-I) Y(I+ND4P1)=Y(ND4P1-I) 20 CONTINUE C====================================================================== C== == C== POPULATE THE BOTTOM 1/2 OF THE ELLIPSE WITH N/4*2 POINTS. == C== == C====================================================================== ND2M1=2*ND4-1 DO 30 I=1,ND2M1 X(I+ND4P1+ND4)=-X(I+1) Y(I+ND4P1+ND4)=-Y(I+1) 30 CONTINUE C====================================================================== C== == C== ROTATE THE ELLIPSE SO THAT IT HAS AN ORIENTATION ANGLE OF ANG == C== DEGREES W.R.T THE X-AXIS OF THE NEW FRAME. == C== ALSO SHIFT THE ELLIPSE CENTER TO (XC,YC). == C== == C====================================================================== CANGR=COS(ANG*PI/180.) SANGR=SIN(ANG*PI/180.) DO 40 I=1,NUSE XX=X(I) YY=Y(I) X(I)=CANGR*XX-SANGR*YY+XC Y(I)=SANGR*XX+CANGR*YY+YC 40 CONTINUE RETURN END SUBROUTINE ELLCON(ANGLE,A,B,AREA,R12,XINT,YINT,XMAX,YMAX, * YATXM,XATYM) C C LIBRARY-ROUTINE C C 05/AUG/1980 C C.J. KOST SIN C C================================================================ C== == C== ELLCON: CONVERTS THE INPUT ELLIPSE COORDINATES: ANGLE,A,B== C== INTO THE OUTPUT QUANTITIES: AREA,R12,XINT,YINT, == C== XMAX,YMAX,YATXM,XATYM. THESE PARAMETERS ARE == C== DESCRIBED BELOW: == C== == C== WRITTEN BY ARTHUR HAYNES, TRIUMF U.B.C., SEPT 1978. == C== == C== INPUT PARAMETERS: ANGLE,A,B (R*4). == C== OUTPUT PARAMETERS: AREA,R12,XINT,YINT,XMAX,YMAX,YATXM, == C== XATYM (R*4). == C== == C== ANGLE : ANGLE OF THE ELLIPSE SEMI-MAJOR AXIS A RELATIVE == C== TO THE POSITIVE X-AXIS IN DEGREES. == C== A : SEMI-MAJOR AXIS OF THE ELLIPSE. == C== B : SEMI-MINOR AXIS OF THE ELLIPSE. == C== NOTE: B MAY BE GREATER THAN A. == C== AREA : AREA OF THE ELLIPSE = PI*A*B. == C== == C== FOR THE FOLLOWING PARAMETERS THE X-AXIS AND Y-AXIS ARE == C== ASSUMED TO PASS THROUGH THE CENTER OF THE ELLIPSE: == C== == C== R12 : CORRELATION COEFFICIENT OF THE ELLIPSE. == C== XINT : POSITIVE X-INTERCEPT OF THE ELLIPSE. == C== YINT : POSITIVE Y-INTERCEPT OF THE ELLIPSE. == C== NOTE: EQUATION OF ELLIPSE IS: == C== X**2/XI**2-2*R12*X*Y/XI/YI+Y**2/YI**2 = 1 == C== WHERE XI = X-INTERCEPT, YI = Y-INTERCEPT. == C== XMAX : MAXIMUM X-COORDINATE OF THE ELLIPSE (X 1/2 WIDTH)== C== YMAX : MAXIMUM Y-COORDINATE OF THE ELLIPSE (Y 1/2 WIDTH)== C== YATXM : Y COORDINATE OF THE ELLIPSE AT XMAX. == C== XATYM : X COORDINATE OF THE ELLIPSE AT YMAX. == C== == C== THE ABOVE PARAMETERS ARE CALCULATED FROM ANGLE,A,B USING == C== THE FOLLOWING FORMULAS: == C== == C== AREA = PI*A*B. == C== XMAX = SQRT((A*COS(ANGLE))**2+(B*SIN(ANGLE))**2). == C== YMAX = SQRT((A*SIN(ANGLE))**2+(B*COS(ANGLE))**2). == C== XINT = (A*B)/YMAX IF YMAX .NE. 0. == C== = 0. IF YMAX = 0. == C== YINT = (A*B)/XMAX IF XMAX .NE. 0. == C== = 0. IF XMAX = 0. == C== R12 = SIGN*SQRT(1.-(A*B/XMAX/YMAX)**2), == C== WHERE SIGN = SIGN(1.,SINCOS*(A**2-B**2)) == C== = +1 OR -1 IF SINCOS*(A**2-B**2) NE 0.== C== = 0 IF SINCOS*(A**2-B**2) = 0.== C== SINCOS = SIN(ANGLE) * COS(ANGLE). == C== == C== YATXM = R12*YMAX. == C== XATYM = R12*XMAX. == C== == C================================================================ CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC IMPLICIT REAL*8 (A-H,O-Z) REAL ANGLE,A,B,AREA,R12,XINT,YINT,XMAX,YMAX,YATXM,XATYM DATA PI/3.14159265358979323846264D0/ C================================================================ C== CALCULATE AREA. == C================================================================ AB=ABS(A*B) AREA=PI*AB C================================================================ C== CALCULATE XMAX AND YMAX. == C================================================================ ANGL=ANGLE*PI/180.D0 COSA=DCOS(ANGL) SINA=DSIN(ANGL) SINCOS=SINA*COSA ASQ=A*A BSQ=B*B COSSQ=COSA*COSA SINSQ=SINA*SINA DXMAX=DSQRT(ASQ*COSSQ+BSQ*SINSQ) DYMAX=DSQRT(ASQ*SINSQ+BSQ*COSSQ) XMAX=DXMAX YMAX=DYMAX C================================================================ C== CALCULATE XINT AND YINT. == C================================================================ XINT=0. IF(DYMAX.GT.0.D0)XINT=AB/DYMAX YINT=0. IF(DXMAX.GT.0.D0)YINT=AB/DXMAX C================================================================ C== CALCULATE R12. == C================================================================ COEFF=0.D0 IF(DXMAX*DYMAX.GT.0.D0)COEFF=1.D0-(AB/DXMAX/DYMAX)**2 IF(COEFF.LT.0.D0)COEFF=0.D0 SIGNR=SINCOS*(ASQ-BSQ) R12=0. IF(SIGNR.NE.0.D0)R12=DSQRT(COEFF)*DSIGN(1.D0,SIGNR) C================================================================ C== CALCULATE YATXM AND XATYM. == C================================================================ YATXM=R12*YMAX XATYM=R12*XMAX RETURN END SUBROUTINE ELLFIT(XIN,YIN,NPTS,FRAC,XCENTR,YCENTR,A,B,ANGLE2,*) C C LIBRARY-ROUTINE C C 05/AUG/1980 C C.J. KOST SIN C C==================================================================== C== C== ELLFIT C== C== ELLFIT FITS AN ELLIPSE TO THE NPTS INPUT DATA POINTS (XIN,YIN) C== SUCH THAT THE ELLIPSE CONTAINS "FRAC" OF THE INPUT POINTS. C== THE CENTER OF THE ELLIPSE (XCENTR,YCENTR), THE MAJOR AND MINOR C== AXIS (A & B), AND THE ANGLE IN DEGREES WITH THE X-AXIS (ANGLE, C== WHERE 0=0). C* ALSO: ROUTINE RETURNS XCM (X CENTER OF MASS) AND YCM (Y C* CENTER OF MASS). C* IN CASE OF BAD INPUT DATA, ROUTINE RETURNS THROUGH 1. C* C********************************************************************** CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC REAL XIN(NPTS),YIN(NPTS) IF(NPTS.LT.2 .OR. IFLAG.LT.1 .OR. IFLAG.GT.2) RETURN 1 C C****** CALCULATE SUMS OF: (XIN), (YIN), (XIN**2), (YIN**2), C****** AND (XIN*YIN) C SUMX=0.0 SUMY=0.0 SUMXSQ=0.0 SUMYSQ=0.0 SUMXY=0.0 DO 20 I=1,NPTS SUMX=SUMX+XIN(I) SUMY=SUMY+YIN(I) SUMXSQ=SUMXSQ + XIN(I)*XIN(I) SUMYSQ=SUMYSQ + YIN(I)*YIN(I) SUMXY =SUMXY + XIN(I)*YIN(I) 20 CONTINUE C C******** CALCULATE CENTER OF MASSES. C XCM = SUMX/NPTS YCM = SUMY/NPTS IF(IFLAG.EQ.2) GO TO 40 C C****** FIND A AND B FOR LEAST Y**2 LINE Y=AX+B C A = (NPTS*SUMXY - SUMX*SUMY) / (NPTS*SUMXSQ - SUMX*SUMX) B = (SUMY - A*SUMX) / NPTS RETURN C C******** FIND A, B, AND C FOR LEAST PERPENDICULAR**2 LINE C******** A*X + B*Y + C = 0 C 40 A=0 B=1 C=YCM D=NPTS*SUMXSQ - NPTS*SUMYSQ - SUMX*SUMX + SUMY*SUMY E=NPTS*SUMXY - SUMX*SUMY IF(ABS(D).LT.1.E-6 .AND. ABS(E).LT.1.E-6) RETURN F = E*E / (D*D + 4*E*E) ISIGN=1 IF(D.GT.0) ISIGN=-1 ASQ = (1. + ISIGN*SQRT(1.-4.*F)) / 2. A=SQRT(ASQ) B=1. IF(A.GT.1.E-6) B=SIGN(1.,-E) * SQRT(1.-ASQ) IF(A.GT.1.E-6 .AND. ABS(D).GT.1.E-6) B=((ASQ+ASQ-1.)*E) / (A*D) C= -(A*SUMX + B*SUMY)/NPTS RETURN END