module OW_mseq
!		  Maximum length sequence assuming 2,3,5 distinct values
!
!       [ms]=MSEQ(baseVal,powerVal[,shift,whichSeq])
!
!       OUTPUT:
!       ms = generated maximum length sequence, of length basisVal^powerVal-1
!
!       INPUT:
!		  baseVal  -nuber of sequence levels (2,3, or 5 allowed)
!		  powerVal -power, so that sequence length is baseVal^powerVal-1
!		  shift    -cyclical shift of the sequence
!		  whichSeq -sequence istantiation to use 
!		  (numer of sequences varies with powreVal - see the code)

! (c) Giedrius T. Buracas, SNL-B, Salk Institute
! Register values are taken from: WDT Davies, System Identification
! for self-adaptive control. Wiley-Interscience, 1970
! When using mseq code for design of FMRI experiments, please, cite:
! G.T.Buracas & G.M.Boynton (2002) Efficient Design of Event-Related fMRI 
! Experiments Using M-sequences. NeuroImage, 16, 801-813.
 

subroutine owmseq(baseVal,powerVal,shift,whichSeq)

optional, integer				:: shift,whichSeq
integer 					:: baseVal,powerVal
integer					:: bitNum
integer, allocatable			:: register(:)
real					:: ms(:)

integer, dimension (1:179)	:: index=(/2, &
  4, &
  6, &
  8, &
 10, &
 12, &
 14, &
 18, &
 22, &
 26, &
 30, &
 32, &
 34, &
 38, &
 42, &
 46, &
 50, &
 52, &
 54, &
 56, &
 58, &
 62, &
 66, &
 70, &
 74, &
 78, &
 82, &
 86, &
 90, &
 94, &
 98, &
104, &
110, &
116, &
120, &
124, &
128, &
132, &
136, &
140, &
144, &
148, &
152, &
156, &
160, &
164, &
168, &
174, &
180, &
186, &
190, &
192, &
194, &
198, &
202, &
206, &
210, &
214, &
218, &
222, &
226, &
230, &
234, &
238, &
242, &
246, &
250, &
254, &
258, &
264, &
270, &
276, &
282, &
288, &
294, &
300, &
306, &
312, &
318, &
324, &
330, &
336, &
342, &
348, &
354, &
360, &
366, &
372, &
378, &
384, &
390, &
396, &
402, &
408, &
414, &
420, &
426, &
434, &
442, &
444, &
446, &
450, &
454, &
458, &
462, &
466, &
470, &
474, &
478, &
482, &
486, &
490, &
494, &
498, &
502, &
506, &
510, &
514, &
518, &
522, &
526, &
532, &
538, &
544, &
550, &
556, &
562, &
568, &
574, &
580, &
586, &
592, &
598, &
604, &
610, &
616, &
622, &
628, &
634, &
640, &
646, &
652, &
658, &
664, &
670, &
676, &
682, &
688, &
694, &
702, &
710, &
718, &
726, &
734, &
742, &
750, &
758, &
766, &
774, &
776, &
780, &
784, &
788, &
790, &
794, &
796, &
798, &
802, &
804, &
806, &
808, &
810, &
814, &
816, &
820, &
824, &
826, &
828, &
832 /)


integer, dimension (1:179)	:: seq =  (/  1,2, &
           1,3, &
 	 2,3, &
           1,4, &
 	 3,4, &
           2,5, &
	 3,5, &
	 1,2,3,5, &
	 2,3,4,5, &
	 1,2,4,5, &
	 1,3,4,5, &
           1,6, &
	 5,6, &
	 1,2,5,6, &
	 1,4,5,6, &
	 1,3,4,6, &
 	 2,3,5,6, &
           1,7, &
	 6,7, &
	 3,7, &
	 4,7, &
	 1,2,3,7, &
	 4,5,6,7, &
	 1,2,5,7, &
	 2,5,6,7, &
	 2,3,4,7, &
	 3,4,5,7, &
	 1,3,5,7, &
	 2,4,6,7, &
	 1,3,6,7, &
	 1,4,6,7, &
	 2,3,4,5,6,7, &
	 1,2,3,4,5,7, &
	 1,2,4,5,6,7, &
	 1,2,3,5,6,7, &
           1,2,7,8, &
	 1,6,7,8, &
	 1,3,5,8, &
	 3,5,7,8, &
	 2,3,4,8, &
	 4,5,6,8, &
	 2,3,5,8, &
	 3,5,6,8, &
	 2,3,6,8, &
	 2,5,6,8, &
	 2,3,7,8, &
	 1,5,6,8, &
	 1,2,3,4,6,8, &
	 2,4,5,6,7,8, &
	 1,2,3,6,7,8, &
	 1,2,5,6,7,8, &
           4,9, & 
	 5,9, &
	 3,4,6,9, &
	 3,5,6,9, &
	 4,5,8,9, &
	 1,4,5,9, &
	 1,4,8,9, &
	 1,5,8,9, &
	 2,3,5,9, &
	 4,6,7,9, &
	 5,6,8,9, &
	 1,3,4,9, &
	 2,7,8,9, &
	 1,2,7,9, &
	 2,4,7,9, &
	 2,5,7,9, &
	 2,4,8,9, &
	 1,5,7,9, &
	 1,2,4,5,6,9, &
	 3,4,5,7,8,9, &
	 1,3,4,6,7,9, &
	 2,3,5,6,8,9, &
	 3,5,6,7,8,9, &
	 1,2,3,4,6,9, &
	 1,5,6,7,8,9, &
	 1,2,3,4,8,9, &
	 1,2,3,7,8,9, &
	 1,2,6,7,8,9, &
	 1,3,5,6,8,9, &
	 1,3,4,6,8,9, &
	 1,2,3,5,6,9, &
	 3,4,6,7,8,9, &
	 2,3,6,7,8,9, &
	 1,2,3,6,7,9, &
	 1,4,5,6,8,9, &
	 1,3,4,5,8,9, &
	 1,3,6,7,8,9, &
	 1,2,3,6,8,9, &
	 2,3,4,5,6,9, &
	 3,4,5,6,7,9, &
	 2,4,6,7,8,9, &
	 1,2,3,5,7,9, &
	 2,3,4,5,7,9, &
	 2,4,5,6,7,9, &
	 1,2,4,5,7,9, &
	 2,4,5,6,7,9, &
	 1,3,4,5,6,7,8,9, &
	 1,2,3,4,5,6,8,9, &
           3,10, &
	 7,10, &
	 2,3,8,10, &
	 2,7,8,10, &
	 1,3,4,10, &
	 6,7,9,10, &
	 1,5,8,10, &
	 2,5,9,10, &
	 4,5,8,10, &
	 2,5,6,10, &
	 1,4,9,10, &
	 1,6,9,10, &
	 3,4,8,10, &
	 2,6,7,10, &
	 2,3,5,10, &
	 5,7,8,10, &
	 1,2,5,10, &
	 5,8,9,10, &
	 2,4,9,10, &
	 1,6,8,10, &
	 3,7,9,10, &
	 1,3,7,10, &
	 1,2,3,5,6,10, &
	 4,5,7,8,9,10, &
	 2,3,6,8,9,10, &
	 1,2,4,7,8,10, &
	 1,5,6,8,9,10, &
	 1,2,4,5,9,10, &
	 2,5,6,7,8,10, &
	 2,3,4,5,8,10, &
	 2,4,6,8,9,10, &
	 1,2,4,6,8,10, &
	 1,2,3,7,8,10, &
	 2,3,7,8,9,10, &
	 3,4,5,8,9,10, &
	 1,2,5,6,7,10, &
	 1,4,6,7,9,10, &
	 1,3,4,6,9,10, &
	 1,2,6,8,9,10, &
	 1,2,4,8,9,10, &
	 1,4,7,8,9,10, &
	 1,2,3,6,9,10, &
	 1,2,6,7,8,10, &
	 2,3,4,8,9,10, &
	 1,2,4,6,7,10, &
	 3,4,6,8,9,10, &
	 2,4,5,7,9,10, &
	 1,3,5,6,8,10, &
	 3,4,5,6,9,10, &
	 1,4,5,6,7,10, &
	 1,3,4,5,6,7,8,10, &
	 2,3,4,5,6,7,9,10, &
	 3,4,5,6,7,8,9,10, &
	 1,2,3,4,5,6,7,10, &
	 1,2,3,4,5,6,9,10, &
	 1,4,5,6,7,8,9,10, &
	 2,3,4,5,6,8,9,10, &
	 1,2,4,5,6,7,8,10, & 
	 1,2,3,4,6,7,9,10, &
 	 1,3,4,6,7,8,9,10, &
           9,11, &
           6,8,11,12, &
           9,10,12,13, &
           4,8,13,14, &
           14,15, &
           4,13,15,16, &
           14,17, &
           11,18, &
           14,17,18,19, &
           17,20, &
           19,21, &
           21,22, &
           18,23, &
           17,22,23,24, &
           22,25, &
           20,24,25,26, &
           22,25,26,27, &
           25,28, &
           27,29, &
           7,28,29,30/)


	if (not(present(whichSeq))) whichSeq=1
	if (not(present(shift)))    shift=1

          baseVal = 2

	bitNum=baseVal^powerVal-1
	allocate(register(powerVal)); register=1
          lengtap=1

	if (powerVal .le. 30) then
		select case (powerVal)
			case (2)
				if (whichSeq > lengtap .or. whichSeq<1) then
					write(0,*) ' WARNING: wrapping arround!' 
					whichSeq=mod(whichSeq,lengtap)+1;
			   	end if;
				select case (whichSeq)
					case
			case (3)
				lengtap=2
			case (4)
				lengtap=2
			case (5)
				lengtap=6
			case (6)
				lengtap=6
			case (7)
				lengtap=18
			case (8)
				lengtap=16
			case (9)
				lengtap=48
			case (10)
				lengtap=60
		end select

		allocate(ms(bitNum)); ms=0.
		if (not(present).whichSeq) then
			call owexec_time (rnd)
			seed(1)=int(rnd/10.)
			call owexec_time (rnd)
			seed(2)=int(rnd+2000.)
			call random_seed(put=seed)
			call random_number(rnd)
			whichSeq=ceiling(rnd*lengtap)
		else
			if (whichSeq > lengtap .or. whichSeq<1) then
				write(0,*) ' WARNING: wrapping arround!' 
				whichSeq=mod(whichSeq,lengtap)+1;
		   	end if;
		end if;

		allocate(weights(powerVal)); weights=0.
		weights(tap(whichSeq).No)=1;
  
		for i=1:bitNum
   		! calculating next digit with modulo powerVal arithmetic
   		! register, (tap(1).No)

			ms(i)=mod(weights*register+baseVal,baseVal);
  ! updating the register
  			register=[ms(i);register(1:powerVal-1)];
		end do

		ms=ms(:);
		if ~isempty(shift),
   			shift=mod(shift, length(ms));
   			ms=[ms(shift+1:end); ms(1:shift)];
		end if

		if baseVal==2,     
  			ms=ms*2-1;
		else if baseVal==3, 
  			ms(ms==2)=-1;
		else if baseVal==5, 
  			ms(ms==4)=-1;
  			ms(ms==3)=-2;
		else
			call seperr("wrong baseVal!");
		end if
	else
		call seperr ("M-sequence is not defined for powerVal > 30"))
	end if
