macro see_oscillations kT='0.00' * APPLication comis QUIT real function EvenCos(phi) common/pawpar/par(3) EvenCos = par(1)+2.0*par(2)*cos(2.0*phi)+ + 2.0*par(3)*cos(4.0*phi) end QUIT APPLication comis QUIT real function EvenSin(phi) common/pawpar/par(3) EvenSin = par(1)+2.0*par(2)*sin(2.0*phi)+ + 2.0*par(3)*sin(4.0*phi) end QUIT * set hwid 5 set xwin 0.01 set ywin 0.01 set bwid 2 set vsiz 0.6 set asiz 0.6 set ylab 1.0 set bcol 1 set csiz 0.6 op utit op ndat vec/crea par(3) r vec/crea errpar(3) r vec/crea fourier_coefficients(5,3) vec/crea dfourier_coefficients(5,3) title 'KT = '//[kT] vec/rea phi,ro2,rs2,ros2,rl2 'Kt'//[kT]//'.dat' ! OC -/%/(2) first=1 * do ipar=2,5 case [ipar] in (1) par='lam' zone=1 ymax=[lammax] ymin=[lammin] (2) par='ro2' zone=1 temp = $SIGMA([par](8)) ymax=[temp]+3 ymin=[temp]-3 lab='R?O!^2' ylab=14 (3) par='rs2' zone=2 temp = $SIGMA([par](8)) ymax=[temp]+3 ymin=[temp]-3 lab='R?S!^2' ylab=20 (4) par='ros2' zone=3 temp = $SIGMA([par](1)) ymax=[temp]+3 ymin=[temp]-3 lab='R?OS!^2' ylab=-2 (5) par='rl2' zone=4 temp = $SIGMA([par](8)) ymax=[temp]+3 ymin=[temp]-3 lab='R?L!^2' ylab=18 endcase if ([first].eq.1) then zone 2 2 first = 0 else zone 2 2 [zone] s endif if (([ipar].eq.4).or.([ipar].eq.5)) then set yval 0.2 else set yval 200 endif if (([ipar].eq.3).or.([ipar].eq.5)) then set xval -8.9 else set xval 0.3 endif mess [ymin] [ymax] 2dhisto 100 [par] 40 0. 3.2 40 [ymin] [ymax] hi/plo 100 tic * hplo/err phi [par] ? 'd'//[par] $VDIM(phi) [markr] .4 s hplo/err phi [par] ? ? $VDIM(phi) 29 .4 s if ([icorrect].eq.2) then text 50 [ylab] [lab] 0.7 endif if ([ipar].eq.4) then atit '[f] (rad)' 'R^2! (fm^2!)' elseif ([ipar].eq.5) then atit '[f] (rad)' ' ' mess 'hi there' elseif ([ipar].eq.1) then atit ' ' '[l]' elseif ([ipar].eq.2) then atit ' ' 'R^2! (fm^2!)' endif exec see_oscillations#fit_it [ipar] 1 [par] enddo exec see_oscillations#output_coefficients return ************* macro fit_it [ipar] [icorrect] [par] case [ipar] in (1) | lambda exitm (2) | ro2 fitfunc = 'EvenCos' (3) | rs2 fitfunc = 'EvenCos' (4) | ros2 fitfunc = 'EvenSin' (5) | rl2 fitfunc = 'EvenCos' endcase * nd = $VDIM([par]) vec/crea fake_errors([nd]) r [nd]*0.01 vec/fit phi [par] fake_errors [fitfunc] 0N 3 par ! ! ! errpar fun/plo [fitfunc] 0 3.2 cs * vec/copy par(1) fourier_coefficients([ipar],1) vec/copy par(2) fourier_coefficients([ipar],2) vec/copy par(3) fourier_coefficients([ipar],3) vec/copy errpar(1) dfourier_coefficients([ipar],1) vec/copy errpar(2) dfourier_coefficients([ipar],2) vec/copy errpar(3) dfourier_coefficients([ipar],3) return ************* macro output_coefficients APPLIcation comis QUIT subroutine printout(lun) integer lun vector fourier_coefficients(5,3) vector dfourier_coefficients(5,3) character*4 name dimension name(5) * data (name(ipar),ipar=1,5)/'junk',' Ro2',' Rs2','Ros2',' Rl2'/ name(1) = 'junk' name(2) = ' Ro2' name(3) = ' Rs2' name(4) = 'Ros2' name(5) = ' Rl2' if (lun.ne.6) then open(unit=lun,file='Fourier.coefficients',status='unknown') endif c write(lun,190) c write(lun,192) * note there is nothing in index1=1 do ipar=2,5 do iorder=0,4,2 index2=(iorder)/2+1 write(lun,191)name(ipar),iorder, + fourier_coefficients(ipar,index2), + dfourier_coefficients(ipar,index2) enddo write(lun,*)' ' enddo write(lun,193) + fourier_coefficients(2,2) - + fourier_coefficients(3,2) + + 2.0*fourier_coefficients(4,2), + sqrt(dfourier_coefficients(2,2)**2 + + dfourier_coefficients(3,2)**2 + + 2.0*dfourier_coefficients(4,2)**2) write(lun,*)' ' 190 format(14x,'RP-Uncorrected',12x, + 'RP-Corrected',12x,'RP-Corrected') 192 format(14x,'StandardCoulomb',10x, + 'StandardCoulomb',10x,'BowlerCoulomb') 191 format(1x,a4,'_',i1,5x,f8.4,' +/- ',f6.4) 193 format(' SumRule',4x,f8.4,' +/- ',f6.4) write(lun,*)'SumRule is Eq.(29) in the Symmetry paper:' write(lun,*)'Ro2_2 - Rs2_2 + 2*Ros2_2 = ?' if (lun.ne.6) then close(unit=lun) endif return end QUIT call printout(6) call printout(20) return ************* macro create_vectors vec/crea lam(12) r vec/crea ro2(12) r vec/crea rs2(12) r vec/crea ros2(12) r vec/crea rl2(12) r vec/crea dlam(12) r vec/crea dro2(12) r vec/crea drs2(12) r vec/crea dros2(12) r vec/crea drl2(12) r vec/crea phi(12) r return *************** macro read_all subdir do iang=1,12 exec see_oscillations#read_one [subdir] [iang] enddo return **************** macro read_one subdir iang exec [subdir]//'/ang'//[iang].vec ***vec/cre phi(5) r 45 90 135 180 0 vec/cop angfit(1) phi([iang]) vec/cop lambda_fit(1) lam([iang]) vec/cop Rout2_fit(1) ro2([iang]) vec/cop Rside2_fit(1) rs2([iang]) vec/cop Rlong2_fit(1) rl2([iang]) vec/cop dlambda_fit(1) dlam([iang]) vec/cop dRout2_fit(1) dro2([iang]) vec/cop dRside2_fit(1) drs2([iang]) vec/cop dRlong2_fit(1) drl2([iang]) vec/cop Routside2_fit(1) ros2([iang]) vec/cop dRoutside2_fit(1) dros2([iang]) return