* Fri Dec 17 13:17:02 EST 1993 /plas5/h1/ksg/src/test avecur.f C averages currents subroutine avecur( cur, jdata, ipk, kstat, jtlmod, jclk, rmark, isav, iave) integer*2 ipk(4), kstat, jtlmod, jclk, isav, iave real rmark(3), cur(512) integer*2 jdata(512) go to ( 100,200,300,400,500,600,700,800),jtlmod write (8,*)"error in avecur, jdata, jtlmod =", jtlmod stop(1) 100 call avemd1( cur, jdata, ipk, kstat, jtlmod, jclk, rmark, isav, iave) return 200 call avemd2( cur, jdata, ipk, kstat, jtlmod, jclk, rmark, isav, iave) return 300 call avemd3( cur, jdata, ipk, kstat, jtlmod, jclk, rmark, isav, iave) return 400 call avemd4( cur, jdata, ipk, kstat, jtlmod, jclk, rmark, isav, iave) return 500 call avemd5( cur, jdata, ipk, kstat, jtlmod, jclk, rmark, isav, iave) return 600 call avemd6( cur, jdata, ipk, kstat, jtlmod, jclk, rmark, isav, iave) return 700 call avemd7( cur, jdata, ipk, kstat, jtlmod, jclk, rmark, isav, iave) return 800 call avemd8( cur, jdata, ipk, kstat, jtlmod, jclk, rmark, isav, iave) return end subroutine | avemd2( cur, jdata, ipk, kstat, jtlmod, jclk, rmark, isav, iave) integer*2 ipk(4), kstat, jtlmod, jclk, isav, iave c avedm2 is for jtlmod = 2, m modes c cur is currents as measued by instrument (128,4) c ipk is estimated peak of distribution function c isav is flag to 1 average, 2 reset, 3 ignore c iave is flag to 1 return average, 2 do not return average c IMPLICIT INTEGER NUMCH, NUMCUP, PEAKCH, INCH INTEGER NUMCH, NUMCUP, PEAKCH, INCH, NUM_POINTS, NUM_CH parameter( NUMCH=60, NUMCUP=3, PEAKCH=20, NUM_CH=9 ) parameter( INCH=128, NUM_POINTS=3) c parameter( NUMCH=80, NUMCUP=4, PEAKCH=20, NUM_CH=8, INCH=128) c NUM_CH is the number of subchannels c NUM_POINTS is the number of points used in fit c numch is the number of saved channels c numsch is the number of subchannels per a normal channel == NUM_CH c numcup is number of cups considered c peakch is location of peak channel in saved array c integer*2 ipk(NUMCUP) real*4 rmarko(3) real*4 rmark(3), cur(INCH,NUMCUP) integer*2 jdata (INCH,NUMCUP) real*4 acur((NUMCH+1)*NUM_CH,NUMCUP) real*4 aicur((NUMCH+1)*NUM_CH,NUMCUP) real*4 apkoff(NUMCUP) real*4 bpkoff(NUMCUP) integer ipka(NUMCUP) c if the data is in doubt ddo not give it much weight real weight(9) logical linit integer numpt integer numsch integer ii,iii integer ip, iip, m1, m2, kk, i1, i2 integer kstato, jclko, icup, ioff c integer itlmod integer ich, ipkoff, ifine, ipkn c integer ifine0 integer max0, min0 real*4 curn, w c real*4 aweight c real*4 aq, sqrt, x1, x2, y1, y2, aweight, curn, w integer i,j,k,l integer ia( NUM_POINTS) real*8 mat8( NUM_POINTS, NUM_POINTS) real*4 mat4( NUM_POINTS, NUM_POINTS) real*4 F( 6) c real*4 xxx( NUMCUP) c real*4 xxxx( NUMCUP) c integer ipa( NUMCUP) data weight / 1.0, 8*0.0001/ c if (linit ) go to 500 linit = .true. c numppp = PEAKCH*NUM_CH numpt = NUM_POINTS numsch = NUM_CH c ii is subchannel number ii = 0 c do each equation ( real channel ) do 400 j=1, numpt do 10 i=1,numpt ia(i) = 0 10 continue do 200 k=1,numsch c iii is power of subchannel number iii = 1 do 120 l=1,numpt c ia is sum of power of subchannel number c i.e. this part otf the real current ia(l) = ia(l) + iii iii = iii*ii 120 continue ii = ii + 1 200 continue do 220 i=1,numpt mat8( i, j) = ia( i) 220 continue 400 continue call printm( mat8, numpt, numpt) call gaussj( mat8, numpt, numpt) call printm( mat8, numpt, numpt) c copy over to normal percision do 440 j=1,numpt do 420 i=1,numpt mat4( j, i) = real( mat8( i, j)) 420 continue 440 continue iip = (numpt-1)/2 m1 = iip*numsch m2 = m1 + numsch - 1 500 continue c first are we averaging this current c itlmod=jtlmod c write (8,*) "isav ",isav,", iave ", iave go to ( 1000, 2000, 3000) isav c error with isav write (8,*) "error in avecur, isave =", isav stop 1 c average this data 1000 continue kk = kstat kk = mod( kk,16) if ( .not.((kk .eq. 3) .or. (kk .eq. 15))) return if (( kstato .ne. 0) .and.( kstat .ne. kstato)) return kstato = kstat jclko = jclk rmarko(1) = rmark(1) rmarko(2) = rmark(2) rmarko(3) = rmark(3) c do each cup seperately do 1500 icup=1,NUMCUP ip = ipk(icup) c do not average if ther is no peak, in reason if ((ip + PEACKCH .ge. INCH) .or. (ip .le. 3 + PEAKCH )) | go to 1500 88 format(10f8.0) ip = ip - 1 do 1010 i=1,5 c move to peak current, not peak in distribution func if ( cur( ip, icup) .ge. cur( ip+1, icup)) go to 1020 ip = ip+1 1010 continue 1020 continue c comput peak in divided space c three points 1130 continue call multv( mat4, cur(ip -1, icup), f, NUM_POINTS, numpt) if ( f(3) .eq. 0.0 ) then pkoff = numpt * 1.5 else pkoff = -f(2)/(2.0*f(3)) end if go to 1190 c now to copy currents 1190 continue call multv( mat4, cur( ip-iip, icup), f, NUM_POINTS, numpt) ipkoff = pkoff + 0.5 dpkoff = pkoff - ipkoff ppk = (ip -iip)*numsch + pkoff if ((ipkoff .lt. 1) .or. (ipkoff .gt. numpt*numsch)) | ipkoff = numsch c aweight = 1./numch c copy data over ioff = (ip - PEAKCH)*numsch + ipkoff - numsch*2 i1 = ip -1 - PEAKCH ifine = m1 - ipkoff if ( ifine .le. 0 ) then ifine = numsch + ifine i1 = i1 + 1 if ( ifine .le. 0 ) then ifine = numsch + ifine i1 = i1 + 1 end if else if ( ifine .gt. numsch ) then ifine = ifine - numsch i1 = i1 - 1 end if i2 = min0( INCH-1, i1+NUMCH-numpt, NUMCH+ioff/numsch-2) cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c xxxx( icup) = pkoff c ipa( icup) = ip c three points 1230 continue c ifine0 = ifine do 1238 ich=i1,i2 c copy over each of the input channels curn = cur( ich, icup) if ( curn .le. 0.) then ifine = ifine + numsch go to 1238 end if call multv( mat4, cur( ich , icup), f, NUM_POINTS, numpt) call normf( cur( ich+1,icup), f, NUM_POINTS, numpt,m1,m2) w = weight( jdata( ich, icup)/256 +1) do 1234 i=m1,m2 acur( ifine, icup) = acur( ifine, icup) + | w*amax1( 0., (f(1)+i*(f(2)+i*f(3)))) aicur( ifine, icup) = aicur( ifine, icup) + w ifine = ifine+1 1234 continue 1238 continue if ((pkoff .lt. 1.) .or. (pkoff .lt. numpt*numsch)) | pkoff = numsch*numpt*0.5 apkoff(icup) = apkoff(icup) + pkoff + (ip-PEAKCH)*numsch bpkoff(icup) = bpkoff(icup) + ppk ipka( icup) = ipka( icup) + ioff c ii0 = ifine0 + (numsch*PEAKCH)/2 c ii2 = ii0 + PEAKCH*numsch c ii3 = PEAKCH*numsch go to 1400 1400 continue 1500 continue ipkn = ipkn + 1 go to 5000 c clear averages 2000 continue kstato = 0 ipkn = 0 do 2010 i=1,NUMCUP ipka(i) = 0 apkoff(i) = 0.0 bpkoff(i) = 0.0 2010 continue aicur(1,1) = 0 call mvc( aicur(1,1), aicur(2,1), numsch*4*NUMCUP*(NUMCH+1)-4 ) acur(1,1) = 0. call mvc( acur(1,1), acur(2,1), numsch*4*NUMCUP*(NUMCH+1)-4 ) c write (8,*) "arrays cleared" c icup=1 c write (8,*) "icup ",icup c write (8,88) (acur(i,icup),aicur(i,icup),i=1,NUMCH*numsch) c write (8,*) " " go to 5000 c ignore spectra 3000 continue go to 5000 c return average? 5000 continue go to ( 6000, 7000), iave c error write (8,*) "In avecur error with iave =", iave stop 1 c average 6000 continue kstat = 0 if ( ipkn .lt. 1) go to 10000 do 6500 icup=1,NUMCUP ioff = (bpkoff(icup) + 0.5)/ipkn ip = ioff/numsch ipk(icup) = ip i1 = max0( 1, ioff/numsch-PEAKCH) if ( ioff .eq. 0 ) i1 = INCH i2 = min0( i1+NUMCH, INCH) - 1 C ifine = ip*numsch - ioff ifine = ip*numsch - bpkoff(icup)/ipkn - 0.5 if ( ifine .le. 0 ) then ifine = ifine + numsch i1 = i1 + 1 if ( ifine .le. 0 ) then ifine = ifine + numsch i1 = i1 + 1 end if else if ( ifine .gt.numsch ) then ifine = ifine - numsch i1 = i1 - 1 end if cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc ifine = min0( numch, max0( 1, ifine )) c ifine0 = ifine c copy averages do 6100 ich=1,i1-1 cur( ich, icup) = -1.e-20 6100 continue do 6300 ich=i1,i2 cur( ich, icup) = 0. do 6200 i=1,numsch if ( aicur( ifine, icup) .gt. 0) then cur( ich, icup) = cur( ich, icup) + | acur( ifine, icup)/aicur( ifine, icup) end if ifine = ifine+1 6200 continue 6300 continue do 6400 ich=i2+1,INCH cur( ich, icup) = 0. 6400 continue cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 6500 continue jclk = jclko kstat = kstato rmark(1) = rmarko(1) rmark(2) = rmarko(2) rmark(3) = rmarko(3) go to 10000 c no average 7000 continue go to 10000 10000 continue return end